MDV10

#load libraries
library(igraph)
library(curl)
library(rgdal)

## Loading required package: sp
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj

library(rgeos)

## rgeos version: 0.3-8, (SVN revision 460)
##  GEOS runtime version: 3.3.3-CAPI-1.7.4 
##  Polygon checking: TRUE

library(maptools)

## Checking rgeos availability: TRUE

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)

#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)

# 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 <- newOrbnona %over% OrbConvex

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

# Attemp to use sp::open function
# ramphsover <- data.frame(ramphshull)
# coordinates(ramphsover) <-~lon+lat