Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
functions and code for plotting examples of the O point proximity metric
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment