Load Libraries

require(igraph)
require(curl)
require(plyr)
require(maptools)
require(rgdal)
require(leaflet)
require(rgeos)
require(rgdal)

Load Data

orbis_edges <-read.csv(curl("https://stacks.stanford.edu/file/druid:mn425tz9757/orbis_edges_0514.csv"),head=TRUE)
orbis_nodes<-read.csv(curl("https://stacks.stanford.edu/file/druid:mn425tz9757/orbis_nodes_0514.csv"),head=TRUE)
ramphs <-read.csv(curl("https://gist.githubusercontent.com/sfsheath/5c5987269e8aad412416/raw/a12dda2b8a681f1ab01421f68ac237ab591ff0f9/roman-amphitheaters.csv"))

Convert Data to igraph Object

R will not accept the data as it was loaded because of the missing x and y values. To work around this and use all possible data, these columns will be excluded from the igraph object.

temp.df <- orbis_nodes[,c("id","label")]
o.nw  <- graph.data.frame(orbis_edges, vertices=temp.df, directed=TRUE)

Calculate Shortest Paths to Rome

In order to account for the fact that some trips are easier to make than other shorter trips, this calculation will be weighted by the “expense” column of the orbis data.

shortest.paths(o.nw,V(o.nw)[V(o.nw)$label == "Roma"], weights = E(o.nw)$expense) -> o.nw.torome
o.nw <- set.vertex.attribute(o.nw, "torome", index = V(o.nw), value = o.nw.torome)

Use igraph Calculation to Create a Dataframe

In order to plot the shortest paths calculation on a map, the data format will need to be changed to a spatial object. This process begins by extracting the columns of interest from the igraph object as a new dataframe.

data.frame(id = V(o.nw)$name, torome = V(o.nw)$torome) -> tmp.df
join(tmp.df, orbis_nodes, by="id") -> o.torome

Calculate Percentiles of Distance to Rome

The quantile function can break up the data into 10 groups. The 90th percentile of sites closest to Rome will be represented by decile #1, since these are the 10% of sites with the lowest distance to Rome.

o.torome$decile <- as.integer(with(o.torome, cut(torome, breaks=quantile(torome, probs=seq(0,1, by=0.10)), include.lowest=TRUE)))

Convert Dataframe into Spatial Object

R will not allow coordinates to be assigned unless both coordinate columns have all values filled.

torome.sp <- o.torome
torome.sp <- subset(torome.sp, x > 0)
coordinates(torome.sp) <- ~ y + x
proj4string(torome.sp) <- CRS('+proj=longlat +datum=WGS84')

Calculation of ConvexHull

Because this function will draw a polygon around all of the data points, it is necessary to create a subset so that the polygon will only be drawn around the data points nearest to Rome.

o.nearest <- o.torome[o.torome$decile=='1',]
o.nearest.sp <- o.nearest
o.nearest.sp <- subset(o.nearest.sp, x > 0)
coordinates(o.nearest.sp) <- ~ y + x
proj4string(o.nearest.sp) <- CRS('+proj=longlat +datum=WGS84')
gConvexHull(o.nearest.sp) -> convex

Determine which Amphitheatres within the Convex Hull

Using the ramphs data loaded at the start, it will be possible to use the sp::over function to select only the amphitheatres within the convex hull polygon. This first requires converting the ramphs data into a spatial object so that R can make the calculation using two objects in the same format.

coordinates(ramphs) <- ~ longitude + latitude
proj4string(ramphs) <- CRS('+proj=longlat +datum=WGS84')
over(ramphs, as(convex, "SpatialPolygons")) -> nearest.ramphs

Select Only Amphitheatres within Convex Hull and Create Spatial Object

The over function created a column (column number 36) with ‘NA’ for all amphitheatres outside of the convex hull and a number for all amphitheatres within it. Once only the amphitheatres with a complete case in column 36, these amphitheatres within the convex hull can be converted to a spatial object for mapping.

ramphs.new <- data.frame(ramphs, nearest.ramphs)
ramphs.new[complete.cases(ramphs.new[,36]),] -> ramphs.convex
coordinates(ramphs.convex) <- ~ longitude + latitude
proj4string(ramphs.convex) <- CRS('+proj=longlat +datum=WGS84')

Create a Leaflet Map

All data has been converted to spatial objects in preparation for plotting as a leaflet map. The convex hull was plotted as the first layer because it blocked the ability to click for popups on all points within it when it was mapped as one of the final layers. The points nearest to Rome were mapped in blue, while the other points were mapped in red. Amphitheatres within the convex hull were mapped in green. Each layer of points was also set to include a popup with the site name.

leaflet() %>% addPolygons(data=convex, color='blue') %>% addCircles(data= torome.sp[torome.sp$decile >1,], popup=torome.sp$label, color='red') %>% addCircles(data= torome.sp[torome.sp$decile==1,], popup=torome.sp$label, color='red') %>% addCircles(data=ramphs.convex, popup=ramphs.convex$id, color='green') %>% addTiles()