Using Orbis to Identify Amphitheaters Close to Rome

Introduction

The basemap is Johan Åhlfeldt’s “Imperium” map as made for his Digital Atlas of the Roman Empire. The URL pattern for accessing those tiles is http://pelagios.dme.ait.ac.at/tilesets/imperium//{z}/{x}/{y}.png
The Rome boundary is courtesy of UNC’s Ancient World Mapping Center.
Sites in the Stanford Orbis model of Roman connectivity are then plotted according to their supplied latitude and longitude. That data is available at http://purl.stanford.edu/mn425tz9757 . The sites in 90th percentile of closeness to Rome as calculated using the igraph R function shortest.paths() are colored in blue. Those farther away are in red.
The convex hull of the closest sites is plotted. It was calculated using the gConvexHull() function from the rgeos package.
That convex hull is used to identify the Roman amphitheaters that lie within it. The over() function in the sp package handles that. These amphitheaters are then plotted in green.

Map of Amphitheaters “Close” to Rome as Defined by the Stanford Orbis Model of Connectivity in the Roman World

#load libraries
library(curl)
library(leaflet)
library(igraph)
library(rgeos)
library(sp)

#load data sets
library(cawd)

#setting the network
o.nw  <- graph.data.frame(orbis.edges, vertices=orbis.nodes, directed=TRUE)

#add the x and y back in. with the NA's set to 0
#set.vertex.attribute(o.nw,name = "x", index = as.character(orbis.nodes$id), value = orbis.nodes$x) -> o.nw
#set.vertex.attribute(o.nw,name = "y", index = as.character(orbis.nodes$id), value = orbis.nodes$y) -> o.nw

#set.vertex.attribute(o.nw,"x",index = V(o.nw)[is.na(x)], 0) -> o.nw
#set.vertex.attribute(o.nw,"y",index = V(o.nw)[is.na(y)], 0) -> o.nw

#finding the shortest paths referring to Roma, by expenses, and assigning to a new object
shortest.paths(o.nw,V(o.nw)[title == "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)

#new dataframe that has two columns: the label (name of place) and the distance to Rome
data.frame(label = V(o.nw)$title, torome = V(o.nw)$torome) -> tmp.df

quantile(V(o.nw)$torome, probs = 0:60/60)[5] -> q.num

V(o.nw)$title[V(o.nw)$torome >= q.num] -> o.above.q
data.frame(title = o.above.q) -> o.above.q
merge(o.above.q,orbis.nodes, by.x = "title", by.y = "title") -> o.above.q
#o.above.q <- o.above.q[!is.na(o.above.q$x),]
o.above.q.sp <- o.above.q
coordinates(o.above.q.sp) <- ~ longitude + latitude
proj4string(o.above.q.sp) <- CRS("+proj=longlat +datum=WGS84")

V(o.nw)$title[V(o.nw)$torome < q.num] -> o.below.q
data.frame(title = o.below.q) -> o.below.q
merge(o.below.q,orbis.nodes, by.x = "title", by.y = "title") -> o.below.q
#o.below.q <- o.below.q[!is.na(o.below.q$x),]
o.below.q.sp <- o.below.q[o.below.q$latitude > 0,]
coordinates(o.below.q.sp) <- ~ longitude + latitude
proj4string(o.below.q.sp)=CRS("+proj=longlat +datum=WGS84")


gConvexHull(o.below.q.sp) -> o.below.hull.sp
proj4string(o.below.hull.sp)=CRS("+proj=longlat +datum=WGS84")
sp::over(ramphs.sp,o.below.hull.sp) -> these.ramphs
ramphs.sp[!is.na(these.ramphs),] -> these.ramphs.sp
proj4string(these.ramphs.sp)=CRS("+proj=longlat +datum=WGS84")

leaflet() %>% addTiles(urlTemplate = "http://pelagios.dme.ait.ac.at/tilesets/imperium//{z}/{x}/{y}.png",
                       attribution = "Digital Atlas Roman Empire", tileOptions(opacity= .3))  %>%
  addPolygons(data = awmc.roman.empire.200.sp, color = 'black', fillOpacity = .05, stroke = F, ) %>%
  addPolygons(data = o.below.hull.sp, color = 'red', fillOpacity = .1, stroke = F) %>%
  addCircleMarkers(data = o.above.q.sp, color = "red", popup =o.above.q.sp$label, radius = 1) %>%
  addCircleMarkers(data = o.below.q.sp, popup =o.below.q.sp$label, radius = 1) %>%
  addCircleMarkers(data = these.ramphs.sp, color = "green", popup = these.ramphs.sp$label, radius = 1)