Mapping closest amphitheaters base on travel expense

MDV10.R

#load libraries
library(igraph)
library(curl)
library(rgdal)
library(rgeos)
library(maptools)
library(ggplot2)

#load data sets
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)

#subsetting two columns from orbis_nodes
temp.df <- orbis_nodes[,c("id","label")]

#setting the network
o.nw  <- graph.data.frame(orbis_edges, vertices=temp.df, directed=TRUE)

#plotting the network
# plot(o.nw)

#all of the "coastal" values in the "type" column in the edges dataset
#E(o.nw)[type=='coastal' & E(o.nw)$days >1]

#making a subgraph and putting it into an object
#subgraph.edges(o.nw, E(o.nw)[type=='coastal'], delete.vertices = TRUE) -> o.nw.coastal

#plotting the new coastal object
#vertex.label.cex:: relates to the text size
#edge.arrow.size:: excluding the arrows, so less visual mess
#vertex.size:: the size of the node
#plot(o.nw.coastal, 
# vertex.label.cex = .1, 
# edge.arrow.size = 0, 
#vertex.size = 2)

#making a subgraph and putting it into an object
#subgraph.edges(o.nw, E(o.nw)[type=='coastal' | type == 'overseas'], delete.vertices = TRUE) -> o.nw.bysea

#plotting the new object
#plot(o.nw.bysea, 
#     vertex.label.cex = .1, 
#     edge.arrow.size = 0, 
#     vertex.size = 2)

#assigning the degrees to a new object
#degree(o.nw, mode = "in") -> o.nw.degrees

#finding the betweenness by expenses, and assigning to a new object
#betweenness(o.nw, weights = E(o.nw)$expense) -> o.nw.betweenness
#o.nw <- set.vertex.attribute(o.nw, "betweenness", index= V(o.nw), value=o.nw.betweenness)

#finding the closeness by expenses, and assigning to a new object
#closeness(o.nw, mode = "in", weights = E(o.nw)$expense) -> o.nw.closeness
#o.nw <- set.vertex.attribute(o.nw, "closeness", index = V(o.nw), value = o.nw.closeness)

#finding the shortest paths referring to Roma, by expenses, and assigning to a new object
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)

# edge.betweenness.community(o.nw,weights = E(o.nw)$expense, directed = TRUE) -> o.nw.ebc
#evcent(o.nw, directed = FALSE, scale = TRUE, weights = E(o.nw)$expense) -> o.nw.evcent

# make hist of torome route
#hist(V(o.nw)$torome)

_________________

# Find shortest paths, then map the top 20% of the sites
# Merge Roma & coords & shortest path

# Step 1: add shortest.distance to a table
data.frame(label=V(o.nw)$label, torome=V(o.nw)$torome) -> df.tmp

# Step 2: merge torome table and nodes table
orbsShort <- merge(orbis_nodes, df.tmp, by="row.names")

# Step 3: find 20% percentile value: 1.2118
orbShortDate <- quantile(orbsShort$torome, probs = c(0, 0.20, 0.5, 0.95, 1))

# Step 4: subset based on every torome value < 1.2118
newOrb <- subset(orbsShort, torome < 1.2118, select=c(label.x,x,y,torome))
# subet base on evert torome value > 12.40545
farOrb <- subset(orbsShort, torome > 12.40545, select=c(label.x,x,y,torome))

# Step 5: change long lat column name
colnames(newOrb)[colnames(newOrb)=="x"] <- "lat"
colnames(newOrb)[colnames(newOrb)=="y"] <- "lon"

# Step 6: save as csv to import into CartoDB
# write.csv(newOrb, file = "newOrb.csv")

# Step 7: prepare newOrb & farOrb for plotting
# omit NAs
newOrbnona <- na.omit(newOrb)
coordinates(newOrbnona) <- ~lon+lat
proj4string(newOrbnona) <- CRS("+proj=longlat +datum=WGS84")
plot(newOrbnona)

farOrbnona <- na.omit(farOrb)
coordinates(farOrbnona) <- ~y+x
proj4string(farOrbnona) <- CRS("+proj=longlat +datum=WGS84")
plot(farOrbnona)

# Step 8: create convexhull of closest
Orb_hull <- gConvexHull(newOrbnona)
OrbConvex <- spTransform(Orb_hull, CRS("+proj=longlat +datum=WGS84"))
plot(OrbConvex)

# Step 9: read in shapefile as background
romeshape <- readShapeSpatial('/Users/mariafang/Desktop/NYU/2_AncientData/mediterranean-shapefiles/roman_empire_bc_60_extent.shp')
proj4string(romeshape) <- CRS("+proj=longlat +datum=WGS84")

# Step 10: find amphitheaters that are inside the convex full
ramphshull <- over(newOrbnona, OrbConvex)

plot(romeshape)
plot(OrbConvex, add=T)
plot(newOrbnona, cex = .5, col = "red", add=T)
plot(farOrbnona, cex = .25, col = "purple", add=T)

# Failed attempt to use result from over function
ramphsover <- data.frame(ramphshull)
coordinates(ramphsover) <-~lon+lat