##Point-in-polygons to designate regional clusters in global Jurassic PBDB occurrence data #Library with point-in-polygon function library(SDMTools) #Object with all occurrences occurrences <- read.csv('occurrences.csv', header=TRUE, stringsAsFactors=FALSE) #Add new column to hold geographic cluster designation occurrences[,"geo_cluster"] <- rep(0, nrow(occurrences)) #Plot of global occurrences plot(occurrences[,1:2], main='Global Occurrences', xlab='Paleolongitude', ylab='Paleolatitude', las=1, type='n', axes=FALSE) axis(1, at=seq(-180, 180, by=10), labels=seq(-180, 180, by=10)) axis(2, at=seq(-180, 180, by=10), labels=seq(-180, 180, by=10), las=1) points(occurrences[,1:2]) #Read polygon coordinates that were created manually from the plot of global occurrences by using a screen shot of the plot in illustrator and then extracted using GraphClick. poly.coords <- read.csv("PolygonCoordinates.csv", header=TRUE) #Extract coordinates vector for each polygon NorthAmericapolypnts <- poly.coords[poly.coords$Polygon=="Polygon 1", 1:2] Europepolypnts <- poly.coords[poly.coords$Polygon=="Polygon 2", 1:2] SouthAmericapolypnts <- poly.coords[poly.coords$Polygon=="Polygon 4", 1:2] MiddleEastpolypnts <- poly.coords[poly.coords$Polygon=="Polygon 5", 1:2] NewZealandpolypnts <- poly.coords[poly.coords$Polygon=="Polygon 6", 1:2] #Plot polygons polygon(NorthAmericapolypnts, lwd = 2, border = "red") text(-40, 80, labels = "North America", col = "red") polygon(Europepolypnts, lwd = 2, border = "blue") text(-5, 5, labels = "Europe", col = "blue") polygon(SouthAmericapolypnts, lwd = 2, border = "green") text(-70, -50, labels = "South America", col = "green") polygon(MiddleEastpolypnts, lwd = 2, border = "orange") text(80, -20, labels = "Middle East", col = "orange") polygon(NewZealandpolypnts, lwd = 2, border = "light blue") text(20, -80, labels = "New Zealand", col = "light blue") #Create objects with the occurrences that occur in each polygon NorthAmericapolyoccs <- pnt.in.poly(occurrences[,1:2], NorthAmericapolypnts) Europepolyoccs <- pnt.in.poly(occurrences[,1:2], Europepolypnts) SouthAmericapolyoccs <- pnt.in.poly(occurrences[,1:2], SouthAmericapolypnts) MiddleEastpolyoccs <- pnt.in.poly(occurrences[,1:2], MiddleEastpolypnts) NewZealandpolyoccs <- pnt.in.poly(occurrences[,1:2], NewZealandpolypnts) #Add list of points for each polygon to the geo_cluster vector in the occurrences dataframe by modern continent/country: North America, Europe, South America, Middle East, New Zealand, Unassigned #First, get rid of magic number by naming the column number with the geoclust vector geoClustercolumn <- 41 occurrences[NorthAmericapolyoccs$pip==1, geoClustercolumn] <- "North America" occurrences[Europepolyoccs$pip==1, geoClustercolumn] <- "Europe" occurrences[SouthAmericapolyoccs$pip==1, geoClustercolumn] <- "South America" occurrences[MiddleEastpolyoccs$pip==1, geoClustercolumn] <- "Middle East" occurrences[NewZealandpolyoccs$pip==1, geoClustercolumn] <- "New Zealand" occurrences[occurrences$geo_cluster==0, geoClustercolumn] <- "Unassigned" #Make tables showing number of occurrences in each cluster for each time interval #By stage with(occurrences, table(stage, geo_cluster)) #By PBDB 10my bin with(occurrences, table(pbdb_bin, geo_cluster)) #Write new csv file that includes column with cluster designations write.table(occurrences, "clusteredOccurrences.csv", row.names=FALSE, sep=",")