Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
COM(P)ADRE blog code
#################################################
## Illuminating the dark corners of demography
## W.K. Petry
##
## README: This code reproduces the analyses shown on the COM(P)ADRE blog. The analyses were run
## at a spatial resolution of 1/12° (>9.3 million cells to cover the Earth!), and will take several hours
## to run the necessary calculations on a standard desktop/laptop. It's recommended that the code be run
## at a coarser spatial resolution first, increasing the resolution only when accuracy is needed. See the
## 'res' parameter under the 'Preliminaries' heading.
#################################################
## Preliminaries
#################################################
# load necessary packages
library(tidyverse) # data wrangling
library(sp) # general spatial classes
library(rgdal) # spatial data wrangling
library(raster) # raster/pixel based data
library(rworldxtra) # high resolution world map
library(rgeos) # spatial geometry tools
# set up folder to hold exported files
wd<-"~/COM(P)ADRE/"
setwd(wd)
# get the latest COMPADRE (comp) and COMADRE (coma) databases
comp<-"compadre.RData"
coma<-"comadre.RData"
comp_url<-"http://www.compadre-db.org/Data/CompadreDownload"
coma_url<-"http://www.compadre-db.org/Data/ComadreDownload"
if(!file.exists(comp)) download.file(comp_url,destfile=comp)
if(!file.exists(coma)) download.file(coma_url,destfile=coma)
load(comp)
load(coma)
# set focal years (for subsetting MPM databases)
focalyears<-1966:2016
# set the spatial resolution, coordinate system for analyses
res<-1 # units are degrees, use 'res<-1/12' to reproduce blog analyses
fact<-1/res # aggregation factor for producing smaller maps (defaults to produce 1°
# resolution map; use 'fact<-1' to keep analysis resolution)
p4s<-"+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# plot options
par(mar=c(0.1,0.1,0.1,0.1))
#################################################
## Georeference the MPM databases
#################################################
comp_spdf<-SpatialPointsDataFrame(coords=compadre$metadata %>%
dplyr::select(Lon,Lat) %>%
distinct(.) %>%
filter(!is.na(Lon)&!is.na(Lat)),
data=compadre$metadata %>%
distinct(Lon,Lat,.keep_all=T) %>%
filter(!is.na(Lon)&!is.na(Lat)),
proj4string=CRS(p4s))
coma_spdf<-SpatialPointsDataFrame(coords=comadre$metadata %>%
dplyr::select(Lon,Lat) %>%
distinct(.) %>%
filter(!is.na(Lon)&!is.na(Lat)),
data=comadre$metadata %>%
distinct(Lon,Lat,.keep_all=T) %>%
filter(!is.na(Lon)&!is.na(Lat)),
proj4string=CRS(p4s))
#################################################
## Generate isolation distance using raster
## approximation of a Voronoi diagram
#################################################
## Whole-dataset (i.e. all years)
# set up empty raster template
empty_rast<-raster(xmn=-180,xmx=180,ymn=-90,ymx=90,crs=CRS(p4s),
resolution=res,vals=NULL)
# calculate distance raster only once (automatically stores file in working directory)
if(!file.exists("./comp_dist.bil")){
comp_dist<-distanceFromPoints(empty_rast,comp_spdf,
filename="./comp_dist.bil",
format="EHdr",progress="text",prj=T,overwrite=T)
}else{
comp_dist<-raster("comp_dist.bil")
}
if(!file.exists("./coma_dist.bil")){
coma_dist<-distanceFromPoints(empty_rast,coma_spdf,
filename="./coma_dist.bil",
format="EHdr",progress="text",prj=T,overwrite=T)
}else{
coma_dist<-raster("coma_dist.bil",crs=p4s)
}
#################################################
## Find the most isolated points
#################################################
# (i) Any point on Earth
(comp_farthest_all<-xyFromCell(comp_dist,which.max(comp_dist),spatial=T))
comp_dist
(coma_farthest_all<-xyFromCell(coma_dist,which.max(coma_dist),spatial=T))
coma_dist
# (ii) Any point on land
data(countriesHigh)
countriesHigh<-spTransform(countriesHigh,CRS(p4s))
comp_land<-mask(comp_dist,countriesHigh)
(comp_farthest_land<-xyFromCell(comp_land,which.max(comp_land),spatial=T))
comp_land
coma_land<-mask(coma_dist,countriesHigh)
(coma_farthest_land<-xyFromCell(coma_land,which.max(coma_land),spatial=T))
coma_land
# (iii) Any point on land, excluding Antarctica
countriesHigh_noant<-countriesHigh[countriesHigh@data$ne_10m_adm!="ATA",]
comp_land_noant<-mask(comp_dist,countriesHigh_noant)
(comp_farthest_noant<-xyFromCell(comp_land_noant,
which.max(comp_land_noant),spatial=T))
comp_land_noant
coma_land_noant<-mask(coma_dist,countriesHigh_noant)
(coma_farthest_noant<-xyFromCell(coma_land_noant,
which.max(coma_land_noant),spatial=T)) # Easter Island
coma_land_noant
#################################################
## Make static maps
#################################################
# make ocean polygon to mask isolation distance rasters
# (finds the inverse of the 'countriesHigh' polygons)
worldextent<-as(extent(comp_dist),"SpatialPolygons")
proj4string(worldextent)<-p4s
ocean_mask<-gDifference(worldextent,countriesHigh)
## COMPADRE map
comp_points<-spTransform(do.call("rbind",list(comp_farthest_all[1],comp_farthest_land[1],comp_farthest_noant[1])),
CRSobj=CRS(p4s))
comp_dist_small<-aggregate(comp_dist,fact=fact)
pdf("comp_isolation_current.pdf",height=9,width=15,useDingbats=F)
plot(comp_dist_small,
col=grey(100:1/100),
legend=F,axes=F,box=F,
asp=1,bty="n",interpolate=T)
plot(ocean_mask,col=rgb(30,144,255,maxColorValue=255,alpha=0.15*255),border="transparent",asp=1,
add=T,bty="n")
points(comp_points,pch=rev(c(16,17,15)),col=c("dodgerblue2","#06ad51","#F3CA40"),cex=4)
dev.off()
## COMADRE map
coma_points<-spTransform(do.call("rbind",list(coma_farthest_all[406],coma_farthest_land[1],coma_farthest_noant[1])),
CRSobj=CRS(p4s))
coma_dist_small<-aggregate(coma_dist,fact=fact)
pdf("coma_isolation_current.pdf",height=9,width=15,useDingbats=F)
plot(coma_dist_small,
col=grey(100:1/100),
legend=F,axes=F,box=F,
asp=1,bty="n",interpolate=T)
plot(ocean_mask,col=rgb(30,144,255,maxColorValue=255,alpha=0.15*255),border="transparent",asp=1,
add=T,bty="n")
points(coma_points,pch=rev(c(16,17,15)),col=c("dodgerblue2","#06ad51","#F3CA40"),cex=4)
dev.off()
#################################################
## Calculate isolation and farthest points by year
## (= tiles for animation)
#################################################
if(!file.exists("comp_dist_brick.bil")){
comp_spdf$YearPublication[is.na(comp_spdf$YearPublication)]<-"9999"
comp_spdf_list<-lapply(as.list(focalyears),
function(x){comp_spdf[comp_spdf$YearPublication<=x,]})
# Calculate raster Voronoi for each year
comp_dist_brick<-brick(lapply(comp_spdf_list,distanceFromPoints,
object=empty_rast))
names(comp_dist_brick)<-as.character(focalyears)
comp_dist_brick<-setMinMax(comp_dist_brick)
writeRaster(comp_dist_brick,"comp_dist_brick.bil",format="EHdr")
}else{
comp_dist_brick<-brick("comp_dist_brick.bil")
names(comp_dist_brick)<-as.character(focalyears)
}
if(!file.exists("coma_dist_brick.bil")){
coma_spdf$YearPublication[is.na(coma_spdf$YearPublication)]<-"9999"
coma_spdf_list<-lapply(as.list(focalyears),
function(x){coma_spdf[coma_spdf$YearPublication<=x,]})
# Calculate raster Voronoi for each year
coma_dist_brick<-brick(lapply(coma_spdf_list,distanceFromPoints,
object=empty_rast))
names(coma_dist_brick)<-as.character(focalyears)
coma_dist_brick<-setMinMax(coma_dist_brick)
writeRaster(coma_dist_brick,"coma_dist_brick.bil",format="EHdr")
}else{
coma_dist_brick<-brick("coma_dist_brick.bil")
names(coma_dist_brick)<-as.character(focalyears)
}
## Calculate darkest points
comp_isolationpts_time<-data.frame(year=focalyears,
lat_ww=0,long_ww=0,
lat_land=0,long_land=0,
lat_noant=0,long_noant=0)
coma_isolationpts_time<-data.frame(year=focalyears,
lat_ww=0,long_ww=0,
lat_land=0,long_land=0,
lat_noant=0,long_noant=0)
# (i) whole world
comp_ww<-vector(mode="list")
for(i in 1:nlayers(comp_dist_brick)){
comp_ww[[i]]<-xyFromCell(comp_dist_brick[[i]],
which.max(comp_dist_brick[[i]]),spatial=T)[1]
}
comp_ww<-do.call("rbind",comp_ww)
coma_ww<-vector(mode="list")
for(i in 1:nlayers(coma_dist_brick)){
coma_ww[[i]]<-xyFromCell(coma_dist_brick[[i]],
which.max(coma_dist_brick[[i]]),spatial=T)[1]
}
coma_ww<-do.call("rbind",coma_ww)
# (ii) all land
comp_dist_brick_land<-mask(comp_dist_brick,countriesHigh)
comp_land<-vector(mode="list")
for(i in 1:nlayers(comp_dist_brick_land)){
comp_land[[i]]<-xyFromCell(comp_dist_brick_land[[i]],
which.max(comp_dist_brick_land[[i]]),spatial=T)[1]
}
comp_land<-do.call("rbind",comp_land)
coma_dist_brick_land<-mask(coma_dist_brick,countriesHigh)
coma_land<-vector(mode="list")
for(i in 1:nlayers(coma_dist_brick_land)){
coma_land[[i]]<-xyFromCell(coma_dist_brick_land[[i]],
which.max(coma_dist_brick_land[[i]]),spatial=T)[1]
}
coma_land<-do.call("rbind",coma_land)
# (iii) land, excluding Antarctica
comp_dist_brick_noant<-mask(comp_dist_brick,countriesHigh_noant)
comp_noant<-vector(mode="list")
for(i in 1:nlayers(comp_dist_brick_noant)){
comp_noant[[i]]<-xyFromCell(comp_dist_brick_noant[[i]],
which.max(comp_dist_brick_noant[[i]]),spatial=T)[1]
}
comp_noant<-do.call("rbind",comp_noant)
coma_dist_brick_noant<-mask(coma_dist_brick,countriesHigh_noant)
coma_noant<-vector(mode="list")
for(i in 1:nlayers(coma_dist_brick_noant)){
coma_noant[[i]]<-xyFromCell(coma_dist_brick_noant[[i]],
which.max(coma_dist_brick_noant[[i]]),spatial=T)[1]
}
coma_noant<-do.call("rbind",coma_noant)
# reduce resolution
comp_dist_brick_small<-aggregate(comp_dist_brick,fact=fact)
coma_dist_brick_small<-aggregate(coma_dist_brick,fact=fact)
## prepare map tiles for animation
for(i in 1:nlayers(comp_dist_brick_small)){
# COMPADRE
fname1<-paste0("./comp_stack/",names(comp_dist_brick_small)[i],".jpeg")
comp_points_sub<-do.call("rbind",list(comp_ww[i],comp_land[i],comp_noant[i])) %>%
spTransform(.,CRS(p4s))
jpeg(filename=fname1,height=9,width=15,units="in",res=100)
plot(comp_dist_brick_small[[i]],
col=grey(100:1/100),
legend=F,axes=F,box=F,
asp=1,bty="n",interpolate=T)
plot(ocean_mask,col=rgb(30,144,255,maxColorValue=255,alpha=0.15*255),border="transparent",asp=1,
add=T,bty="n")
points(comp_points_sub,pch=rev(c(16,17,15)),col=c("dodgerblue2","#06ad51","#F3CA40"),cex=4)
title(main=substring(names(comp_dist_brick_small)[i],2),line=0,cex.main=4)
dev.off()
# COMADRE
fname2<-paste0("./coma_stack/",names(comp_dist_brick_small)[i],".jpeg")
coma_points_sub<-do.call("rbind",list(coma_ww[i],coma_land[i],coma_noant[i])) %>%
spTransform(.,CRS(p4s))
jpeg(filename=fname2,height=9,width=15,units="in",res=100)
plot(coma_dist_brick_small[[i]],
col=grey(100:1/100),
legend=F,axes=F,box=F,
asp=1,bty="n",interpolate=T)
plot(ocean_mask,col=rgb(30,144,255,maxColorValue=255,alpha=0.15*255),border="transparent",asp=1,
add=T,bty="n")
points(coma_points_sub,pch=rev(c(16,17,15)),col=c("dodgerblue2","#06ad51","#F3CA40"),cex=4)
title(main=substring(names(coma_dist_brick_small)[i],2),line=0,cex.main=4)
dev.off()
}
#################################################
## Measure isolation changes over time
#################################################
maxchange<-data.frame(year=rep(focalyears,times=2),
dataset=ordered(rep(c("COMPADRE - plants",
"COMADRE - animals"),
each=length(focalyears)),
levels=c("COMPADRE - plants",
"COMADRE - animals")),
ww=c(maxValue(comp_dist_brick)/1000,
maxValue(coma_dist_brick)/1000),
al=c(maxValue(comp_dist_brick_land)/1000,
maxValue(coma_dist_brick_land)/1000),
noant=c(maxValue(comp_dist_brick_noant)/1000,
maxValue(coma_dist_brick_noant)/1000)) %>%
gather(series,max,3:5) %>%
mutate(series=ordered(series,levels=c("ww","al","noant")))
comp_avgdist_ww<-comp_avgdist_land<-comp_avgdist_noant<-coma_avgdist_ww<-coma_avgdist_land<-coma_avgdist_noant<-numeric(length=nlayers(comp_dist_brick))
for(i in 1:nlayers(comp_dist_brick)){
comp_avgdist_ww[i]<-cellStats(comp_dist_brick[[i]],median)
comp_avgdist_land[i]<-cellStats(comp_dist_brick_land[[i]],median)
comp_avgdist_noant[i]<-cellStats(comp_dist_brick_noant[[i]],median)
coma_avgdist_ww[i]<-cellStats(coma_dist_brick[[i]],median)
coma_avgdist_land[i]<-cellStats(coma_dist_brick_land[[i]],median)
coma_avgdist_noant[i]<-cellStats(coma_dist_brick_noant[[i]],median)
}
medchange<-data.frame(year=rep(focalyears,times=2),
dataset=ordered(rep(c("COMPADRE - plants",
"COMADRE - animals"),
each=length(focalyears)),
levels=c("COMPADRE - plants",
"COMADRE - animals")),
ww=c(comp_avgdist_ww,coma_avgdist_ww)/1000,
al=c(comp_avgdist_land,coma_avgdist_land)/1000,
noant=c(comp_avgdist_noant,coma_avgdist_noant)/1000) %>%
gather(series,median,3:5) %>%
mutate(series=ordered(series,levels=c("ww","al","noant")))
isochange<-full_join(maxchange,medchange,by=c("year","dataset","series")) %>%
gather(stat,isolation_km,4:5) %>%
filter(year!=1965)
(gg_isochange<-ggplot(isochange,aes(x=year,y=isolation_km,color=series,
linetype=stat))+
geom_line(size=1.5)+
scale_x_continuous(name="Year")+
scale_y_continuous(name="Isolation (km)",
limits=c(0,20100),expand=c(0,0))+
scale_colour_manual(name="",values=c("dodgerblue2","#06ad51","#F3CA40"),
labels=c("World","Land only",
"Land, excl. Antarctica"))+
scale_linetype_manual(name="",values=c("solid","dashed"),
labels=c("Maximum","Median"))+
facet_grid(.~dataset)+
theme_bw()+
theme(legend.position="top",
plot.margin=unit(c(1,1,1,1),"cm"),
axis.title=element_text(size=28,face="bold"),
axis.text=element_text(size=20),
legend.text=element_text(size=20),
legend.key.width=unit(15,"mm"),
strip.text=element_text(size=18,face="bold"),
panel.grid=element_blank()))
pdf("isochange.pdf",width=12,height=8,useDingbats=F)
gg_isochange
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.