functions and code for plotting examples of the O point proximity metric

pointOverlap.R

# point overlap, LDVA
# load 
library(sp)
library(dplyr)
library(ggplot2)
library(ggmap)
library(fields)
library(cowplot)
library(rgbif)
library(gridExtra)
# source the script for calculating O from Dan Warren's GH
source("https://raw.githubusercontent.com/danlwarren/arc-extensions/master/nncluster.R")
# source this function for creating overlapping polygons with their respective points

squaresOvals <- function(overlapPct){
#create 1x1 square
dat1<- matrix(c(1,1,1,2,2,2,2,1),ncol=2,byrow = T)
ch1 <- chull(dat1)
coords1 <- dat1[c(ch1, ch1[1]), ]  # closed polygon
sp_poly1 <- SpatialPolygons(list(Polygons(list(Polygon(coords1)), ID=1)))
# second square, shifted to the right of the first one by the value of overlapPct
dat2<- matrix(c(1,1,1,2,2,2,2,1),ncol=2,byrow = T)
dat2[,1] <- dat2[,1]+(overlapPct)
ch2 <- chull(dat2)
coords2 <- dat2[c(ch2, ch2[1]), ]  # closed polygon
sp_poly2 <- SpatialPolygons(list(Polygons(list(Polygon(coords2)), ID=2)))

# fortify the sp objects
sp1df <- fortify(sp_poly1)
sp2df <- fortify(sp_poly2)
#bind into single df
squares <- bind_rows(sp1df,sp2df)

# generate the random points
rpts1 <- spsample(sp_poly1, type="random", n=1000)
rpts2 <- spsample(sp_poly2, type="random", n=1000)
rpts1df <- rpts1@coords %>% as.data.frame() %>% mutate(sq="1")
rpts2df <- rpts2@coords %>% as.data.frame() %>% mutate(sq="2")
#bind the dfs
rptsB <- bind_rows(rpts1df,rpts2df)

# subset into sets of observations
rptsSubset <- list()  
npointsO <- data.frame(n=seq(50,1000,by=50),Oval=0)
for(i in seq(50,1000,by=50)) {
  rptsSubset[[i]] <- rptsB %>% group_by(sq) %>% sample_n(i)
}  
rptsSubset <-  rptsSubset[!sapply(rptsSubset, is.null)]
# run the nncluseter function to get a matrix of values
npointsO <- c()
for(j in 1:length(rptsSubset)){
  npointsO[j] <- nncluster(rptsSubset[[j]][,1:2],rptsSubset[[j]]$sq)[1,2]
}
# bind the columns into a DF
npointsOvals <- data.frame(n=seq(50,1000,by=50),Oval=npointsO)
#save a plot of the overlap
exPlot <- ggplot(squares,aes(x=long,y=lat,color=factor(id)))+geom_polygon(fill=NA)+
  geom_point(data=rptsSubset[[4]],aes(x,y,color=factor(sq)))+
  scale_color_manual(values=c("#FFA373","#A5121E"),guide=FALSE)+
  xlab("x")+ylab("y")

# mean O value for given overlap
meanOval <- mean(npointsOvals$Oval)
#save a gg obj of the O values
OvalPlot <- ggplot(npointsOvals,aes(x=n,y=Oval))+geom_point()+
  ylim(c(0,0.55))+  ylab("O")+  geom_hline(yintercept = meanOval)


return(list(OvalsGrob=OvalPlot,overlapPlot=exPlot))
}

#vector of percentage overlap between the polygons
  # actual overlap is 1 - the value
overlapVec <- c(0, 0.10, 0.25, 0.50, 0.75, 0.90, 1,1.25)

# list of results
ovalsList <- list()
# run the function for the overlap scenarios
for (i in 1:length(overlapVec)){
ovalsList[[i]] <- squaresOvals(overlapVec[i])  
}
# you can ignore the warnings
  
#plot side by side
grid.arrange(ovalsList[[1]]$overlapPlot,ovalsList[[1]]$OvalsGrob,
             ovalsList[[2]]$overlapPlot,ovalsList[[2]]$OvalsGrob,
             ovalsList[[3]]$overlapPlot,ovalsList[[3]]$OvalsGrob,
             ovalsList[[4]]$overlapPlot,ovalsList[[4]]$OvalsGrob,
             ovalsList[[5]]$overlapPlot,ovalsList[[5]]$OvalsGrob,
             ovalsList[[6]]$overlapPlot,ovalsList[[6]]$OvalsGrob,
             ovalsList[[7]]$overlapPlot,ovalsList[[7]]$OvalsGrob,
             ovalsList[[8]]$overlapPlot,ovalsList[[8]]$OvalsGrob,ncol=2,nrow=8)

# Download and clean occurrene data for 5 species
# restricted to Mexico, no duplicate records

Owagleri <- occ_data(scientificName="Ortalis wagleri", hasCoordinate=TRUE,
                    country="MX", hasGeospatialIssue=FALSE)$data
locsDistOwag <- Owagleri %>% distinct(name,decimalLongitude,decimalLatitude) %>% 
  rename(species=name,longitude=decimalLongitude,latitude=decimalLatitude)

Opoliocephala <- occ_data(scientificName="Ortalis poliocephala", hasCoordinate=TRUE,
                     country="MX", hasGeospatialIssue=FALSE)$data
locsDistOpol <- Opoliocephala %>% distinct(name,decimalLongitude,decimalLatitude) %>% 
  rename(species=name,longitude=decimalLongitude,latitude=decimalLatitude)

Tcanescens <- occ_data(scientificName="Marmosa canescens", hasCoordinate=TRUE,
                          country="MX", hasGeospatialIssue=FALSE)$data
locsDistTcan <- Tcanescens %>% distinct(name,decimalLongitude,decimalLatitude) %>% 
  rename(species=name,longitude=decimalLongitude,latitude=decimalLatitude)

Oarenicola <- occ_data(scientificName="Onychomys arenicola", hasCoordinate=TRUE,
                       country="MX", hasGeospatialIssue=FALSE)$data
locsDistOaren <- Oarenicola %>% distinct(name,decimalLongitude,decimalLatitude) %>% 
  rename(species=name,longitude=decimalLongitude,latitude=decimalLatitude)

Nleucodon <- occ_data(scientificName="Neotoma leucodon", hasCoordinate=TRUE,
                       country="MX", hasGeospatialIssue=FALSE)$data
locsDistNleuc <- Nleucodon %>% distinct(name,decimalLongitude,decimalLatitude) %>% 
  rename(species=name,longitude=decimalLongitude,latitude=decimalLatitude)

# bind all the species
allPoints <- bind_rows(locsDistOwag,locsDistOpol,locsDistTcan,locsDistOaren,locsDistNleuc)


# need a projected coordinate system 
allPointsCopy <- allPoints
coordinates(allPointsCopy) <- 2:3
proj4string(allPointsCopy) <- CRS("+proj=longlat +datum=WGS84")
# project to the Mexico INEGI lambert conformal using Proj4 settings copied from
# http://spatialreference.org/ref/sr-org/mexico-inegi-lambert-conformal-conic/
allPointsDMX <-  spTransform(allPointsCopy, CRS("+proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +a=6378137 +b=6378136.027241431 +units=m +no_defs "))
#back to DF
allPointsDMXdf <- as.data.frame(allPointsDMX)

#calculate O
mxspeciesO <- nncluster(allPointsDMXdf[,2:3],allPointsDMXdf$species)
# clean up the table for inset plotting
mxspeciesRd <- round(mxspeciesO[sapply(mxspeciesO, is.numeric)],2) %>% 
  dplyr::select(Op=1,Tc=2,Oa=3,Nl=4)
mxspeciesRd <- mxspeciesRd[-5,]
rownames(mxspeciesRd) <- NULL
mxspeciesRd$sp <-   c("Ow","Tc","Op","Oa")
forOtable <- mxspeciesRd%>% select(sp,everything(.))

# get a base map
mxmap <- get_map(location = c(-149,16, -56,28),maptype="terrain")

#plot
  ggmap(mxmap)+
  geom_point(data=allPoints,aes(x=longitude,y=latitude,
                                fill=species,shape=species),color="black")+
  scale_shape_manual(values=c(21,22,23,24,23))+
  inset(tableGrob(forOtable,rows=NULL,theme=ttheme_minimal(base_size = 12)),
        xmin = -115,xmax=-100,ymin=10,ymax=15)