wpetry/DemogDarkCorners.R

## DemogDarkCorners.R
#################################################
## 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()
	#################################################
	## 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()