Select Git revision
figure_2 40.74 KiB
# CREATE FIGURE 2
# load libraries
library(raster)
library(maptools)
library(mapplots)
library(scales)
library(plotrix)
library(Cairo)
library(png)
library(plyr)
library(shapes)
library(shape)
setwd("C:/bio/R_Backup/cloudforests/figures/figures_17072018")
source("pcolor_scale.R")
# SET FOLDER CONNECTIONS
INPUT <-"C:/bio/Bio_Backup/TIN R Skripte/GIT/cloudforests/data/"
OUTPUT<-"C:/bio/Bio_Backup/TIN R Skripte/GIT/cloudforests/output/"
# read in the essential files
pp <- shapefile(paste0(INPUT,"cloud_forest_regions_reshaped_aspoints.shp"))
tc <- read.csv (paste0(INPUT,"cf_statistics_newid.txt"),header=T,sep="\t")
full <- read.csv (paste0(INPUT,"vertsferns_by_cfID_speclist_101518_newid.csv"),sep=";")
tcX <- read.csv (paste0(INPUT,"cf_statistics_full_newid.csv"),sep=";")
# do some data quentching on the input files
{
CFRegCount_per_spec <- ddply(full, .(taxogroup,scientificname,cf_ass,cf_end), summarise, cfIDCount = length(na.omit(cf_id)))
CFRegCountIsOne <- subset(CFRegCount_per_spec, cfIDCount ==1)
CFSingleRegEndemicsInfo <- merge(CFRegCountIsOne, subset(full,select = c("scientificname","cf_id","CFRegionName","CFRegionContinent")),
by = "scientificname", all.x= TRUE, all.y = FALSE)
agg.cfEndemics <- ddply(CFSingleRegEndemicsInfo, .(cf_id, CFRegionName, CFRegionContinent, taxogroup), summarise,
EndSR.cf.assoc = sum(cf_ass),
EndSR.cf.restr = sum(cf_end))
agg.cfEndemics <- agg.cfEndemics[order(-agg.cfEndemics$EndSR.cf.restr),]
biodiv <- ddply(full, .(cf_id, CFRegionName, CFRegionContinent, taxogroup), summarise,
cf_ass = sum(cf_ass),
cf_end = sum(cf_end))
ends<-CFSingleRegEndemicsInfo
full <-full[,-8]
#Rends<-aggregate(full$cf_id~full$scientificname,FUN=function(x){length(x)})
#Rends<-Rends[Rends$`full$cf_id`==1,]
#ends<-full[full$scientificname %in% Rends$`full$scientificname`,]
cff<-full[,1:3]
cff<-cff[duplicated(cff)==FALSE,]
cord<-as.data.frame(pp$ID)
cord$x<-coordinates(pp)[,1]
cord$y<-coordinates(pp)[,2]
colnames(cord)<-c("cf_id","x","y")
biodiv<-merge(biodiv,cord,by.x="cf_id",by.y="cf_id")
}
namematching<-read.csv("Z:/karger/revision1/cf_namematching.csv")
# read in the pa data with all pa's
# some quentching on the pa_stat file
# file with cloud forest regions tcf extent
pa_stat<- shapefile("Z:/karger/revision1/cloud_forest_regions_reshaped [ID]_v16.shp")
pa_stat<- as.data.frame(pa_stat)pa_stat<- pa_stat[,-(c(2,3))]
colnames(pa_stat)[1]<-"ID"
pa_stat<- aggregate(x=pa_stat[,2:37],by=list(pa_stat$ID),FUN=function(x){sum(x)})
colnames(pa_stat)[1]<-"ID"
pa_stat<- merge(namematching,pa_stat,by.x="ID_new",by.y="ID")
pa_stat<- pa_stat[,-c(1,3,4,5)]
# file with cloud forest regions intersection with pa's tcf extent including all pa years
pa_stat_all<- shapefile("Z:/karger/revision1/Intersect [cloud_forest_regions_reshaped [ID]]-[Intersect_cfRegions_IUCN_cf [ID]]_v16.shp")
pa_stat_all<- as.data.frame(pa_stat_all)
pa_stat_all<- pa_stat_all[,-(c(2:6))]
pa_stat_all1<- aggregate(x=pa_stat_all[,2:length(pa_stat_all[1,])],by=list(pa_stat_all$ID),FUN=function(x){sum(na.omit(x))})
pa_stat_mean<- aggregate(x=pa_stat_all[,2:length(pa_stat_all[1,])],by=list(pa_stat_all$ID),FUN=function(x){mean(na.omit(x))})
pa_stat_all<-pa_stat_all1
pa_stat_all[,20:length(pa_stat_all[1,])]<-pa_stat_mean[,20:length(pa_stat_mean[1,])]
colnames(pa_stat_all)[1]<-"ID"
pa_stat_all<- merge(namematching,pa_stat_all,by.x="ID_new",by.y="ID")
pa_stat_all<- pa_stat_all[,-c(1,3,4,5)]
pa_stat_un <-pa_stat[,2:19]-pa_stat_all[,2:19]
pa_stat_un$ID <-pa_stat_all$ID
pa_stat_un<-cbind(pa_stat_un[length(pa_stat_un[1,])],pa_stat_un[,2:(length(pa_stat_un[1,])-1)])
pa_stat_un_sd <-pa_stat_all[,20:37]
pa_stat_un_sd$ID <-pa_stat_all$ID
pa_stat_un_sd <-cbind(pa_stat_un_sd[length(pa_stat_un_sd[1,])],pa_stat_un_sd[,2:(length(pa_stat_un_sd[1,])-1)])
pa_stat_pa <-pa_stat_all[,2:20]
pa_stat_pa$ID <-pa_stat_all$ID
pa_stat_pa <-cbind(pa_stat_pa[length(pa_stat_pa[1,])],pa_stat_pa[,2:(length(pa_stat_pa[1,])-1)])
pa_stat_pa_sd <-pa_stat_all[,20:37]
pa_stat_pa_sd$ID <-pa_stat_all$ID
pa_stat_pa_sd <-cbind(pa_stat_pa_sd[length(pa_stat_pa_sd[1,])],pa_stat_pa_sd[,2:(length(pa_stat_pa_sd[1,])-1)])
# read in the pa data with all pas
# some quentching on the pa_stat file
#pa_stat20012016<- read.csv("Z:/karger/pa_cf_region_intersect_ohne2001_2016.txt",sep = "\t")
##### FIGURE 2 ##########################################################################################################################################################################################################
pdf(paste0(OUTPUT,"Fig2_C_R2.pdf"),height=8.5,width = 6.5)
par(oma=c(4,4,1.4,1),mfrow=c(4,3),tcl=-0.5)#,mar=c(4,4,2,0)
par(mai=c(0,0.2,0.2,0))
single_region_endemics<-ends[ends$cf_end==1,]
#### Data preparation a
{
specs<-full
specs<-merge(specs,single_region_endemics,by.x="scientificname",by.y="scientificname",all.x=T)
specs<-specs[,1:8]
colnames(specs)<-c("scientificname", "cf_id","CFRegionName", "CFRegionContinent","taxogroup", "cf_ass", "cf_end", "cf_srend")
specs$cf_srend<-as.numeric(specs$cf_srend)
specs$cf_srend[specs$cf_srend>=1]<-1
specs$cf_srend[is.na(specs$cf_srend)]<-0
specscf<-specs[specs$cf_ass==1,]
ends <-specs[specs$cf_end==1,]
srends <-specs[specs$cf_srend==1,]
specs_cont <- aggregate(specs$scientificname~specs$CFRegionContinent+specs$taxogroup,FUN=function(x) {length(unique(x))})
specs_end <- aggregate(srends$scientificname~srends$CFRegionContinent+srends$taxogroup,FUN=function(x) {length(unique(x))})
specs_cfend<- aggregate(ends$scientificname~ends$CFRegionContinent+ends$taxogroup,FUN=function(x) {length(unique(x))})
}
specs_stats<-cbind(specs_cont,specs_cfend[,3],specs_end[,3])
colnames(specs_stats)<-c("cont","taxogroup","cf_ass","cf_str","cf_end")write.csv(specs_stats,file=paste0(OUTPUT,"spec_stats.csv"))
continents_IDS<-c("Americas","Africa","Asia-Pacific")
continents_IDS<-sort(continents_IDS)
ids<-c(2,1,3)
for (region in ids)
{
{
contID=continents_IDS[region]
colnames(specs_cont)<-c("cont","taxogroup","scientficname")
colnames(specs_end) <-c("cont","taxogroup","scientficname")
colnames(specs_cfend) <-c("cont","taxogroup","scientficname")
mammals <-specs_cont [specs_cont$cont==contID&specs_cont$taxogroup=="Mammals",]
mammals_end <-specs_end [specs_end$cont==contID&specs_end$taxogroup=="Mammals",]
mammals_cfend<-specs_cfend[specs_cfend$cont==contID&specs_cfend$taxogroup=="Mammals",]
mammals <-cbind(mammals[1,1],mammals[1,3],mammals_cfend[1,3],mammals_end[1,3])
mammals <-as.data.frame(mammals)
colnames(mammals)=c("cf_id","z_ass","z_str","z_end")
cyathea <-specs_cont[specs_cont$cont==contID&specs_cont$taxogroup=="Cyatheaceae",]
cyathea_end <-specs_end[specs_end$cont==contID&specs_end$taxogroup=="Cyatheaceae",]
cyathea_cfend<-specs_cfend[specs_cfend$cont==contID&specs_cfend$taxogroup=="Cyatheaceae",]
cyathea <-cbind(cyathea[1,1],cyathea[1,3],cyathea_cfend[1,3],cyathea_end[1,3])
cyathea <-as.data.frame(cyathea)
colnames(cyathea)=c("cf_id","z_ass","z_str","z_end")
birds <-specs_cont[specs_cont$cont==contID&specs_cont$taxogroup=="Birds",]
birds_end <-specs_end[specs_end$cont==contID&specs_end$taxogroup=="Birds",]
birds_cfend <-specs_cfend[specs_cfend$cont==contID&specs_cfend$taxogroup=="Birds",]
birds <-cbind(birds[1,1],birds[1,3],birds_cfend[1,3],birds_end[1,3])
birds <-as.data.frame(birds)
colnames(birds)=c("cf_id","z_ass","z_str","z_end")
amphies <-specs_cont[specs_cont$cont==contID&specs_cont$taxogroup=="Amphibians",]
amphies_end <-specs_end[specs_end$cont==contID&specs_end$taxogroup=="Amphibians",]
amphies_cfend<-specs_cfend[specs_cfend$cont==contID&specs_cfend$taxogroup=="Amphibians",]
amphies <-cbind(amphies[1,1],amphies[1,3],amphies_cfend[1,3],amphies_end[1,3])
amphies <-as.data.frame(amphies)
colnames(amphies)=c("cf_id","z_ass","z_str","z_end")
#plot(0,ylim=c(-8,8),xlim=c(-8,8),col="white", main=contID,frame=F,axes=F,ylab="",xlab="")
mapx<-function (xlim, ylim, xlab = "Longitude", ylab = "Latitude",
bg = "white", ...)
{
asp <- 1/cos(sum(ylim) * pi/360)
plot(NA, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
asp = asp, ...)
}
mapx(c(-8,8),c(-8,8),bg="white",axes=F,ylab="",xlab="",col="white")
cols<-c(alpha("goldenrod1", 0.5), alpha("royalblue4", 0.5), alpha("darkgreen",0.5), alpha("firebrick",0.5))
n=which(mammals$cf_id==region)
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2),col=cols, labels="",label.dist = 0.5,border=FALSE)
cols<-c(alpha("darkorange1", 0.9), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[1],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0.9), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0)) zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_end[n],birds$z_end[n],cyathea$z_end[n],amphies$z_end[n])
zf<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[1],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0.9), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[2],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0.9), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_end[n],birds$z_end[n],cyathea$z_end[n],amphies$z_end[n])
zf<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[2],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0), alpha("darkgreen",0.9), alpha("firebrick",0))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[3],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0), alpha("black", 0.9), alpha("firebrick",0))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_end[n],birds$z_end[n],cyathea$z_end[n],amphies$z_end[n])
zf<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[3],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0.9))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[4],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0), alpha("black", 0), alpha("black", 0.9))
zs<-c(mammals$z_ass[n],birds$z_ass[n],cyathea$z_ass[n],amphies$z_ass[n])
zs[is.na(zs)]<-0.1
ze<-c(mammals$z_end[n],birds$z_end[n],cyathea$z_end[n],amphies$z_end[n])
zf<-c(mammals$z_str[n],birds$z_str[n],cyathea$z_str[n],amphies$z_str[n])
ze[is.na(ze)]<-0.1
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,0,0,radius=log(sum(zs)+2)*p[4],col=cols, labels="",border=FALSE)
}
cexX<-1
yy=log(5) text("5 sp.",x=0,y=yy,pos=3,col="black",cex=cexX,font=2)
draw.ellipse(x=0,y=0,a = yy, b = yy ,nv=400,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(50)
text("50 sp.",x=0,y=yy-0.35,pos=3,col="black",cex=cexX,font=2)
draw.ellipse(x=0,y=0,a = yy, b = yy ,nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(250)
text("250 sp.",x=0,y=yy-0.5,pos=3,col="black",cex=cexX,font=2)
draw.ellipse(x=0,y=0,a = yy, b = yy ,nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(1000)
text("1000 sp.",x=0,y=yy-0.5,pos=3,col="black",cex=cexX,font=2)
draw.ellipse(x=0,y=0,a = yy, b = yy ,nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
par(xpd=NA)
text(contID,x=0,y=yy+2,pos=3,col="black",cex=1.5,font=2)
if(contID=="Africa")
{
img<-readPNG("C:\\bio\\R_Backup\\cloudforests\\figures\\figures_9042018\\shapes\\mammals2.png")
rasterImage(img,3.8,3.8,9,8.5)
img<-readPNG("C:\\bio\\R_Backup\\cloudforests\\figures\\figures_9042018\\shapes\\birds.png")
rasterImage(img,3,-7,7,-2)
img<-readPNG("C:\\bio\\R_Backup\\cloudforests\\figures\\figures_9042018\\shapes\\treeferns.png")
rasterImage(img,-7,-7,-4,-2)
img<-readPNG("C:\\bio\\R_Backup\\cloudforests\\figures\\figures_9042018\\shapes\\amphibians.png")
rasterImage(img,-7,3,-3,7.5)
}
par(xpd=TRUE)
if (contID=="Asia-Pacific")
{
par(xpd=NA)
lines(x=c(-7,-1),y=c(-7,-1),col=alpha("royalblue4", 0.75),lwd=2)
text(x=-7,y=-7.5,"strict")
lines(x=c(-7,-3),y=c(7,3),col=alpha("firebrick", 0.25),lwd=2)
text(x=-7.2,y=7.5,"associated")
lines(x=c(7,0),y=c(7,0),col=alpha("black", 0.75),lwd=2)
text(x=6,y=7.5,"endemic")
par(xpd=F)
}
if (contID=="Americas")
{
par(xpd=F)
text(-11,9.5,"(a)",cex=2)
par(xpd=NA)
}
}
#### Data preparation b
{
tc2<-tc
tc_pr1<-cbind(tc2[,1:2],(tc2[,3:18]))
tc_un1<-cbind(tc2[,1:2],(tc2[,19:34]))
tc_pr1<-tc_pr1+tc_un1
tc_pr1[,1:2]<-tc_un1[,1:2]
tc_pr2<-cbind(tc2[,1:2],(tc2[,35:50]))
tc_un2<-cbind(tc2[,1:2],(tc2[,51:66]))
tc_pr2<-tc_pr2+tc_un2
tc_pr2[,1:2]<-tc_un2[,1:2]
tc_pr3<-cbind(tc2[,1:2],(tc2[,67:82]))
tc_un3<-cbind(tc2[,1:2],(tc2[,83:98]))
tc_pr3<-tc_pr3+tc_un3 tc_pr3[,1:2]<-tc_un3[,1:2]
tc_pr_pr1<-cbind(tc_pr1[,1:2],100-100*(tc_pr1[,3:18]/tc_pr1[,3]))
tc_pr_pr2<-cbind(tc_pr2[,1:2],100-100*(tc_pr2[,3:18]/tc_pr2[,3]))
tc_pr_pr3<-cbind(tc_pr3[,1:2],100-100*(tc_pr3[,3:18]/tc_pr3[,3]))
tc_proc1<-100-tc_pr_pr1
tc_proc2<-100-tc_pr_pr2
tc_proc3<-100-tc_pr_pr3
ccf<-full[,1:3]
ccf<-ccf[duplicated(ccf$cf_id)==FALSE,]
tc_proc1$ID<-tc_pr1$ID
tc_proc1$Name<-tc_pr1$Name
tc_proc1<-merge(tc_proc1,ccf,by.x="ID",by.y="cf_id")
tc_proc2$ID<-tc_pr2$ID
tc_proc2$Name<-tc_pr2$Name
tc_proc2<-merge(tc_proc2,ccf,by.x="ID",by.y="cf_id")
tc_proc3$ID<-tc_pr3$ID
tc_proc3$Name<-tc_pr3$Name
tc_proc3<-merge(tc_proc3,ccf,by.x="ID",by.y="cf_id")
}
{
contnames<-tc_proc3[,19:20]
id_names<-as.data.frame(pp)
id_names$Name<-gsub("_"," ",id_names$Name)
#id_names<-merge(contnames,id_names,by.x="Name",by.y="Name",all.x=T)
id_names<-id_names[,1:3]
tc_proc4<-pa_stat
#tc_proc4<-tc_proc4[,-1]
colnames(tc_proc4)[1]<-"ID"
tc_sds<-tc_proc4[,c(20:(length(tc_proc4[1,])))]
tc_proc4<-tc_proc4[,-c(20:(length(tc_proc4[1,])))]
tc_proc1<-tc_proc4
tc_proc1[,c(2:(length(tc_proc1[1,])))]<-tc_proc1[,c(2:(length(tc_proc1[1,])))]+tc_sds
tc_proc2<-tc_proc4
tc_proc2[,c(2:(length(tc_proc2[1,])))]<-tc_proc2[,c(2:(length(tc_proc2[1,])))]-tc_sds
tc_proc4<-cbind(tc_proc4[,1],100*(tc_proc4[,2:19]/tc_proc4[,2]))
tc_proc1<-cbind(tc_proc1[,1],100*(tc_proc1[,2:19]/tc_proc1[,2]))
tc_proc2<-cbind(tc_proc2[,1],100*(tc_proc2[,2:19]/tc_proc2[,2]))
colnames(tc_proc4)[1]<-"ID"
colnames(tc_proc1)[1]<-"ID"
colnames(tc_proc2)[1]<-"ID"
tc_proc4<-merge(tc_proc4,namematching,by.x="ID",by.y="ID_org")
tc_proc1<-merge(tc_proc1,namematching,by.x="ID",by.y="ID_org")
tc_proc2<-merge(tc_proc2,namematching,by.x="ID",by.y="ID_org")
tc_proc3<-tc_proc4
}
for (CONT in unique(c("Americas","Africa","Asia-Pacific")))
{
years<-seq(2001,2018,1)
tx_cont1<-tc_proc1[tc_proc1$CFRegionContinent==CONT,]
tx_cont2<-tc_proc2[tc_proc2$CFRegionContinent==CONT,]
tx_cont3<-tc_proc3[tc_proc3$CFRegionContinent==CONT,]
mini<-min(min(tx_cont1[,3:20],na.rm = T),min(tx_cont2[,3:20],na.rm = T),min(tx_cont3[,3:20],na.rm = T))
if(CONT=="Americas")
{
plot(0,ylim=c(90,100),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="")
for(n in 90:100){abline(n,0,col="grey85",lwd=0.5)}
par(xpd=NA)
title(ylab="2001 TCF area",line=3.0,cex.lab=1) text(-11,9.5,"(b)",cex=2)
axis(side=2, at=seq(90,100, by=1),labels=c("90%"," ","92%"," ","94%"," ","96%"," ","98%"," ","100%"),col="grey50",las=1)
par(xpd=F)
col=rgb(90/255, 180/255, 172/255,alpha=1)
text(x=2014,y=90.5,"51%*",col=col,cex=1)
shape::Arrows(2016.2,90,2016.2,89.9, lwd=1,arr.type="triangle",col=col,arr.width=0.15,arr.length=0.15)
}
if (CONT=="Africa"||CONT=="Asia-Pacific")
{
plot(0,ylim=c(90,100),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="")
for(n in 90:100){abline(n,0,col="grey85",lwd=0.5)}
par(xpd=NA)
axis(side=2, at=seq(90,100, by=1),col="grey50",labels=FALSE,cex.lab=1.5,las=1)
par(xpd=F)
}
axis(side=1, at=c(2001:2018),cex=0.5,col="grey50",labels=FALSE)
for (n in 1:length(tx_cont1[,1]))
{
y1<-as.numeric(as.matrix(tx_cont1[n,2:19]))
y2<-as.numeric(as.matrix(tx_cont2[n,2:19]))
y3<-as.numeric(as.matrix(tx_cont3[n,2:19]))
s1<-cbind(y1,y2,y3)
sd1<-apply(s1,1, sd, na.rm = TRUE)
sx1<-apply(s1,1, mean, na.rm = TRUE)
upper=sx1+sd1
lower=sx1-sd1
years=seq(2001,2018,1)
upper=upper[1:18]
lower=lower[1:18]
col=rgb(90/255, 180/255, 172/255,alpha=0.5)
polygon(c(rev(years), years), c(rev(lower), upper), col = col, border = NA)
col=rgb(90/255, 180/255, 172/255,alpha=0.9)
points(sx1~years,col=col,type="l",lwd=0.5)
}
}
#### Data preparation c
{
#tc2<-tc
##tc_pr1<-cbind(tc2[,1:2],(tc2[,3:18]))
##tc_un1<-cbind(tc2[,1:2],(tc2[,19:34]))
##tc_pr2<-cbind(tc2[,1:2],(tc2[,35:50]))
##tc_un2<-cbind(tc2[,1:2],(tc2[,51:66]))
##tc_pr3<-cbind(tc2[,1:2],(tc2[,67:82]))
#tc_un3<-cbind(tc2[,1:2],(tc2[,83:98])
#############################################3
tc_pr1<-pa_stat_pa[,2:18]
tc_pr2<-pa_stat_pa[,2:18]+pa_stat_pa_sd[,2:18]
tc_pr3<-pa_stat_pa[,2:18]-pa_stat_pa_sd[,2:18]
tc_pr1$ID<-pa_stat_pa$ID
tc_pr2$ID<-pa_stat_pa$ID
tc_pr3$ID<-pa_stat_pa$ID
tc_un1<-pa_stat_un[,2:18]
tc_un2<-pa_stat_un[,2:18]+pa_stat_un_sd[,2:18]
tc_un3<-pa_stat_un[,2:18]-pa_stat_un_sd[,2:18]
tc_un1$ID<-pa_stat_un$ID
tc_un2$ID<-pa_stat_un$ID
tc_un3$ID<-pa_stat_un$ID
tc_un_pr1<-cbind(100-100*(tc_un1[,1:17]/tc_un1[,1]))
tc_pr_pr1<-cbind(100-100*(tc_pr1[,1:17]/tc_pr1[,1]))
tc_un_pr2<-cbind(100-100*(tc_un2[,1:17]/tc_un2[,1]))
tc_pr_pr2<-cbind(100-100*(tc_pr2[,1:17]/tc_pr2[,1]))
tc_un_pr3<-cbind(100-100*(tc_un3[,1:17]/tc_un3[,1])) tc_pr_pr3<-cbind(100-100*(tc_pr3[,1:17]/tc_pr3[,1]))
tc_proc1<-tc_pr_pr1/(tc_un_pr1+tc_pr_pr1)
tc_proc2<-tc_pr_pr2/(tc_un_pr2+tc_pr_pr2)
tc_proc3<-tc_pr_pr3/(tc_un_pr3+tc_pr_pr3)
tc_proc1$ID<-pa_stat$ID
tc_proc2$ID<-pa_stat$ID
tc_proc3$ID<-pa_stat$ID
ccf<-full[,1:3]
ccf<-ccf[duplicated(ccf$cf_id)==FALSE,]
#tc_proc1$ID<-tc_pr1$ID
#tc_proc1$Name<-tc_pr1$Name
tc_proc1<-merge(tc_proc1,ccf,by.x="ID",by.y="cf_id",all.x=T)
tc_proc1[,2]<-0.5
#tc_proc2$ID<-tc_proc1$ID
#tc_proc2$Name<-tc_pr2$Name
tc_proc2<-merge(tc_proc2,ccf,by.x="ID",by.y="cf_id",all.x=T)
tc_proc2[,2]<-0.5
#tc_proc3$ID<-tc_proc1$ID
#tc_proc3$Name<-tc_pr3$Name
tc_proc3<-merge(tc_proc3,ccf,by.x="ID",by.y="cf_id",all.x=T)
tc_proc3[,2]<-0.5
}
for (CONT in unique(c("Americas","Africa","Asia-Pacific")))
{
years<-seq(2001,2018,1)
tx_cont1<-tc_proc1[tc_proc1$CFRegionContinent==CONT,]
tx_cont2<-tc_proc2[tc_proc2$CFRegionContinent==CONT,]
tx_cont3<-tc_proc3[tc_proc3$CFRegionContinent==CONT,]
tx_cont2<-tx_cont2[complete.cases(tx_cont1),]
tx_cont3<-tx_cont3[complete.cases(tx_cont1),]
tx_cont1<-tx_cont1[complete.cases(tx_cont1),]
if (CONT=="Americas")
{
plot(0,ylim=c(0,1),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="")
for(n in seq(0,1,0.1)){abline(n,0,col="grey85",lwd=0.5)}
axis(side=2, at=seq(0,1, by=0.1),labels=c("0","","0.2","","0.4","","0.6","","0.8","","1"),col="grey50",las=1)
par(xpd=NA)
title(ylab="Rate ratio",line=3.0,cex.lab=1)
par(xpd=F)
}
if (CONT=="Africa"||CONT=="Asia-Pacific")
{
plot(0,ylim=c(0,1),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="",yaxt="n")
for(n in seq(0,1,0.1)){abline(n,0,col="grey85",lwd=0.5)}
axis(side=2, at=seq(0,1, by=0.1),col="grey50",labels=FALSE)
abline(0.5,0,col="red")
}
axis(side=1, at=c(2001:2018),cex=0.5,col="grey50",labels=FALSE)
n_blue<-0
n_red <-0
for (n in 1:length(tx_cont1[,1]))
{ y1<-as.numeric(as.matrix(tx_cont1[n,2:18]))
y2<-as.numeric(as.matrix(tx_cont2[n,2:18]))
y3<-as.numeric(as.matrix(tx_cont3[n,2:18]))
s1<-cbind(y1,y2,y3)
sd1<-apply(s1,1, sd, na.rm = TRUE)
sx1<-apply(s1,1, mean, na.rm = TRUE)
upper=sx1+sd1#y2
lower=sx1-sd1#y3
years=seq(2002,2018,1)
years=years-0.5
upper=upper[2:18]
lower=lower[2:18]
lower[lower<=0]<-0
col=rgb(0, 0, 0,alpha=0.25)
try(lm1<-lm(((y1+y2+y3)/3)[1:length(y1)]~years))
try(s1<-summary(lm1))
try(assign(paste0("lm_",CONT,"_",n),lm1))
p1 <- s1$coefficients[8]
slp1<- s1$coefficients[2]
if (p1 <= 0.05 && slp1<=0.0)
{
n_red<-n_red+1
col=rgb(0, 0.25, 1,alpha=0.5)
}
if (p1 <= 0.05 && slp1>0.0)
{
n_blue<-n_blue+1
col=rgb(1, 0, 0,alpha=0.5)
}
polygon(c(rev(years), years), c(rev(lower), upper), col = col, border = NA)
if (p1<=0.05)
{
# col=rgb(216/255, 179/255, 101/255,alpha=0.9)
}
points(sx1[2:18]~years,col=col,type="l",lwd=1)
slp1<-0
p1<-1
}
cola="black"
if (CONT=="Americas")
{
abline(0.5,0,col="green")
shape::Arrows(x0=2001,y0=0.35,x1=2001,y1=0.05,col=cola,lwd=1,arr.type="triangle", arr.width=0.1)
shape::Arrows(x0=2001,y0=0.65,x1=2001,y1=0.95,col=cola,lwd=1,arr.type="triangle", arr.width=0.1)
text(x=2001,y=0.95,labels="rate higher inside PA",pos=4,col=cola)
text(x=2001,y=0.05,labels="rate higher outside PA",pos=4,col=cola)
}
if (CONT=="Africa"||CONT=="Asia-Pacific")
{abline(0.5,0,col="green")}
cols<-c("blue","red","grey75")
n_grey<-length(tx_cont1[,1])-(n_red+n_blue)
zs<-c(n_red,n_blue,n_grey)
add.pie(z=zs,2014,0.9,radius=0.1,col=cols, labels="",label.dist = 0.5,border=FALSE)
if (CONT=="Americas")
{
cexll=1.5
text(x=2013.25,y=0.925,"n.s.",cex=0.75)
text(x=2014.75,y=0.925,"-",cex=cexll)
text(x=2014,y=0.85,"+",col="white",cex=cexll)
text(x=2014,y=0.96,"trends",cex=0.75,pos=3)
}
}
nN=0
#### Data preparation d
{
tcX<-merge(tcX,cff,by.x="ID",by.y="cf_id")
tc_all <-aggregate(tcX$X2001~tcX$CFRegionContinent,FUN=function(x){sum(x)})
colnames(tc_all)<-c("CFRegionContinent","cf2001")
tc_all$cf2002 <-aggregate(tcX$X2002~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2003 <-aggregate(tcX$X2003~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2004 <-aggregate(tcX$X2004~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2005 <-aggregate(tcX$X2005~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2006 <-aggregate(tcX$X2006~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2007 <-aggregate(tcX$X2007~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2008 <-aggregate(tcX$X2008~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2009 <-aggregate(tcX$X2009~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2010 <-aggregate(tcX$X2010~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2011 <-aggregate(tcX$X2011~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2012 <-aggregate(tcX$X2012~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2013 <-aggregate(tcX$X2013~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2014 <-aggregate(tcX$X2014~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2015 <-aggregate(tcX$X2015~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc_all$cf2016 <-aggregate(tcX$X2016~tcX$CFRegionContinent,FUN=function(x){sum(x)})[,2]
tc2<-tc
tc2=cbind(tc2[,1:2],tc2[,3:98]/1000)
tc_pr<-cbind(tc2[,1:2],(tc2[,3:18]+tc2[,35:50]+tc2[,67:82])/3)
tc_un<-cbind(tc2[,1:2],(tc2[,19:34]+tc2[,51:66]+tc2[,83:98])/3)
tc_un[is.na(tc_un)]<-1
tc_pr[is.na(tc_pr)]<-1
tc_un_pr<-cbind(tc_un[,1:2],100*(tc_un[,3:18]/tc_un[,3]))
tc_pr_pr<-cbind(tc_pr[,1:2],100*(tc_pr[,3:18]/tc_pr[,3]))
col1<-rgb(51/255, 51/255, 255/255, alpha=1)
col2<-rgb(255/255, 153/255, 0, alpha=1)
tc3<-cbind(tc_un[,1:2],tc_un[,3:18],tc_pr[,3:18])
tc3[,3:34][tc3[,3:34]<=1]<-1
tc3<-tc3[order(tc3$gam_3_cf_un_2001..SUM.+tc3$gam_3_cf_pr_2001..SUM.),]
tc3<-tc3[6:64,]
allspecs<-aggregate(cf_ass~cf_id,sum,data=biodiv)
endemics<-aggregate(cf_end~cf_id,sum,data=biodiv)
tc3 <-merge(tc3,allspecs,by.x="ID",by.y="cf_id")
tc3 <-merge(tc3,endemics,by.x="ID",by.y="cf_id")
tc3 <-tc3[order(tc3$gam_3_cf_un_2001..SUM.+tc3$gam_3_cf_pr_2001..SUM.),]
tt<-merge(tc3,full,by.x="ID",by.y="cf_id",all=F)
tt<-tt[!duplicated(tt$ID),]
tc3$continents<-tt$CFRegionContinent
tc3$specloss<-tc3$cf_ass*(1-((tc3$gam_3_cf_pr_2016..SUM.+tc3$gam_3_cf_un_2016..SUM.)/(tc3$gam_3_cf_pr_2001..SUM.+tc3$gam_3_cf_un_2001..SUM.))^0.75)
tc3 <-tc3[order(tc3$specloss, decreasing = FALSE),]
col_Af<-rgb(27/255,158/255,119/255)
col_Am<-rgb(217/255,95/255,2/255)
col_As<-rgb(0,0,0)
tc3$continents<-as.character(tc3$continents)
tc3$continents[tc3$continents=="Africa"]<-col_Af
tc3$continents[tc3$continents=="Americas"]<-col_Am
tc3$continents[tc3$continents=="Asia-Pacific"]<-col_As
}
df_ext_rates<-c()
for (CONT in unique(c("Americas","Africa","Asia-Pacific")))
{
{
nN<-nN+1 assosiated<-aggregate(full$scientificname[full$CFRegionContinent==CONT]~full$CFRegionContinent[full$CFRegionContinent==CONT],FUN=function(x) {length(unique(x))})
colnames(assosiated)<-c("cf_id","s0")
endemics<-aggregate(ends$scientificname[ends$CFRegionContinent==CONT]~ends$CFRegionContinent[ends$CFRegionContinent==CONT],FUN=function(x) {length(unique(x))})
colnames(endemics)<-c("cf_id","s0")
calcSlos<-function(s0,a,z,loss=FALSE)
{
a0<-a[1:length(a)-1]
a1<-a[2:length(a)]
sst<-c()
sst<-c(sst,s0[1])
for (n in 1:length(a0))
{
sst<-c(sst,sst[n]-sst[n]*(1-(a1[n]/a0[n])^z))
}
ss<-rep(s0,length(a))
if (loss==TRUE)
{
sx<-as.numeric(ss)-as.numeric(sst)
sx=sx*(-1)
}
if (loss==FALSE)
{
sx=sst
}
sx
}
#recalculate aa
contids<-tcX[,c(1,19)]
pa_aa<-pa_stat_all
pa_aa<-merge(pa_aa,contids,by.x="ID_org",by.y="ID")
pa_aa<-aggregate(pa_aa[,2:19],by=list(pa_aa$CFRegionContinent),FUN=function(x){sum(x)})
colnames(pa_aa)<-c("cont","a2001","a2002","a2003","a2004","a2005","a2006","a2007","a2008","a2009","a2010","a2011","a2012","a2013","a2014","a2015","a2016","a2017","a2018")
aa<-pa_aa
#colntca<-c(colnames(tc_all),c("cf2017","cf2018"))
tc_all<-aa
#colnames(tc_all)<-colntca
for (xx in unique(c("assosiated")))
{
# get the species numbers
x<-get(xx)
ax<-as.matrix(tc_all[tc_all$cont==CONT,2:19])[1,]
mx<-c()
percent=FALSE
for (nn in 1:length(x[,1]))
{
m2<-c()
for (z in seq(0.1,1,0.01))
{
a<-as.vector(ax[1:18])
s0<-x[nn,2]
m<-calcSlos(s0,a,z,loss=FALSE)
mm<-c(s0,m)
m2<-rbind(m2,mm)
}
mx<-rbind(mx,m2)
}
mx<-mx[,2:length(mx[1,])]
mx<-as.data.frame(mx,row.names = FALSE)
mxt<-c()
for (row in 1:length(mx[,1]))
{ y<-seq(2001,2018,1)
mxt<-rbind(mxt,cbind(y,as.vector(as.matrix(mx[row,]))))
}
colnames(mxt)<-c("y","s")
mxt<-as.data.frame(mxt)
assign(paste0(xx,"_mxt"),mxt)
}
assosiated_mxt$s<-assosiated$s0-assosiated_mxt$s
assosiated_mxt$s<-assosiated_mxt$s*(-1)
assosiated_mxt<-assosiated_mxt[complete.cases(assosiated_mxt),]
ass_mm<-aggregate(assosiated_mxt$s~assosiated_mxt$y,FUN=function(x){median(x)})
ass_sd<-aggregate(assosiated_mxt$s~assosiated_mxt$y,FUN=function(x){sd(x)})
upper<-ass_mm$`assosiated_mxt$s`+ass_sd$`assosiated_mxt$s`
lower<-ass_mm$`assosiated_mxt$s`-ass_sd$`assosiated_mxt$s`
# try a differnt way of plotting it
byy=2
if(CONT=="Americas"){byy<-10}
m=0
zts<-seq(0.1,1,0.1)
zts<-as.character(zts)
if(CONT=="Americas")
{
plot(0,ylim=c(-30,0),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="")
par(xpd=NA)
title(ylab="Expected species loss",line=3.0,cex.lab=1)
par(xpd=F)
for(n in seq(-30,0,5)){abline(n,0,col="grey85",lwd=0.5)}
}
if(CONT=="Africa"||CONT=="Asia-Pacific")
{
plot(0,ylim=c(-30,0),xlim=c(2001,2018),ylab="",xlab="",axes=FALSE,main="")
par(xpd=F)
for(n in seq(-30,0,5)){abline(n,0,col="grey85",lwd=0.5)}
par(xpd=NA)
}
if(CONT=="Americas"||CONT=="Asia-Pacific")
{
axis(side=1, at=c(2001:2018),cex=0.5,col="grey50")
}
if(CONT=="Africa")
{
axis(side=1, at=c(2001:2018),labels=T,col="grey50")
title(xlab="year",line=2.2,cex.lab=1.5)
}
if(CONT=="Americas")
{
axis(side=2, at=seq(-30,0, by=5),labels=c(" -30",""," -20",""," -10","","0"),col="grey50",las=1)
}
if(CONT=="Asia-Pacific"||CONT=="Africa")
{
axis(side=2, at=seq(-30,0, by=5),col="grey50",labels=F)
}
col=rgb(85/255, 25/255,139/255,alpha=0.01)
seq_rev<-seq(1603,1620,1)
for (n in 1:90)
{
seqq1<-seq(1,18,1)+m
y0<-assosiated_mxt$s[seqq1]#seq(1,16,1)
m<-m+18
seqq2<-seq(1,18,1)+m
y1<-assosiated_mxt$s[seqq1]
y2<-assosiated_mxt$s[seq_rev]
x <-seq(2001,2018,1) try(polygon(c(rev(x), x), c(rev(y2), y0), col = col, border = NA))
seq_rev<-seq_rev-18
}
par(xpd=TRUE)
aa<-aggregate(tc3$gam_3_cf_un_2001..SUM.~tc3$continents,FUN=function(x) {mean(x)})
aa$`tc3$gam_3_cf_un_2001..SUM.`<-aa$`tc3$gam_3_cf_un_2001..SUM.`+aggregate(tc3$gam_3_cf_pr_2001..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2002<-aggregate(tc3$gam_3_cf_un_2002..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2002..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2003<-aggregate(tc3$gam_3_cf_un_2003..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2003..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2004<-aggregate(tc3$gam_3_cf_un_2004..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2004..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2005<-aggregate(tc3$gam_3_cf_un_2005..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2005..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2006<-aggregate(tc3$gam_3_cf_un_2006..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2006..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2007<-aggregate(tc3$gam_3_cf_un_2007..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2007..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2008<-aggregate(tc3$gam_3_cf_un_2008..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2008..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2009<-aggregate(tc3$gam_3_cf_un_2009..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2009..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2010<-aggregate(tc3$gam_3_cf_un_2010..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2010..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2011<-aggregate(tc3$gam_3_cf_un_2011..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2011..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2012<-aggregate(tc3$gam_3_cf_un_2012..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2012..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2013<-aggregate(tc3$gam_3_cf_un_2013..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2013..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2014<-aggregate(tc3$gam_3_cf_un_2014..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2014..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2015<-aggregate(tc3$gam_3_cf_un_2015..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2015..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
aa$a2016<-aggregate(tc3$gam_3_cf_un_2016..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]+aggregate(tc3$gam_3_cf_pr_2016..SUM.~tc3$continents,FUN=function(x) {mean(x)})[,2]
colnames(aa)<-c("cont","a2001","a2002","a2003","a2004","a2005","a2006","a2007","a2008","a2009","a2010","a2011","a2012","a2013","a2014","a2015","a2016")
aa$cont<-c("Asia-Pacific","Africa","Americas")
#recalculate aa
contids<-tcX[,c(1,19)]
pa_aa<-pa_stat_all
pa_aa<-merge(pa_aa,contids,by.x="ID_org",by.y="ID")
pa_aa<-aggregate(pa_aa[,2:19],by=list(pa_aa$CFRegionContinent),FUN=function(x){sum(x)})
colnames(pa_aa)<-c("cont","a2001","a2002","a2003","a2004","a2005","a2006","a2007","a2008","a2009","a2010","a2011","a2012","a2013","a2014","a2015","a2016","a2017","a2018")
aa<-pa_aa
z=0.5
mam_ende<-calcSlos(specs_end[specs_end$cont==CONT&specs_end$taxogroup=="Mammals",3],aa[aa$cont==CONT,2:19],z,loss=T)
mam_asso<-calcSlos(specs_cont[specs_cont$cont==CONT&specs_cont$taxogroup=="Mammals",3],aa[aa$cont==CONT,2:19],z,loss=T)
mam_str <-calcSlos(specs_cfend[specs_cfend$cont==CONT&specs_cfend$taxogroup=="Mammals",3],aa[aa$cont==CONT,2:19],z,loss=T)
df_ext_rates<-rbind(df_ext_rates,mam_asso)
df_ext_rates<-rbind(df_ext_rates,mam_str)
df_ext_rates<-rbind(df_ext_rates,mam_ende)
amp_ende<-calcSlos(specs_end[specs_end$cont==CONT&specs_end$taxogroup=="Amphibians",3],a,z,loss=TRUE)
amp_asso<-calcSlos(specs_cont[specs_cont$cont==CONT&specs_cont$taxogroup=="Amphibians",3],a,z,loss=TRUE)
amp_str <-calcSlos(specs_cfend[specs_cfend$cont==CONT&specs_cfend$taxogroup=="Amphibians",3],a,z,loss=TRUE)
df_ext_rates<-rbind(df_ext_rates,amp_asso)
df_ext_rates<-rbind(df_ext_rates,amp_str)
df_ext_rates<-rbind(df_ext_rates,amp_ende)
tre_ende<-calcSlos(specs_end[specs_end$cont==CONT&specs_end$taxogroup=="Cyatheaceae",3],a,z,loss=TRUE)
tre_asso<-calcSlos(specs_cont[specs_cont$cont==CONT&specs_cont$taxogroup=="Cyatheaceae",3],a,z,loss=TRUE)
tre_str <-calcSlos(specs_cfend[specs_cfend$cont==CONT&specs_cfend$taxogroup=="Cyatheaceae",3],a,z,loss=TRUE)
df_ext_rates<-rbind(df_ext_rates,tre_asso)
df_ext_rates<-rbind(df_ext_rates,tre_str)
df_ext_rates<-rbind(df_ext_rates,tre_ende)
bir_ende<-calcSlos(specs_end[specs_end$cont==CONT&specs_end$taxogroup=="Birds",3],a,z,loss=TRUE)
bir_asso<-calcSlos(specs_cont[specs_cont$cont==CONT&specs_cont$taxogroup=="Birds",3],a,z,loss=TRUE)
bir_str <-calcSlos(specs_cfend[specs_cfend$cont==CONT&specs_cfend$taxogroup=="Birds",3],a,z,loss=TRUE)
df_ext_rates<-rbind(df_ext_rates,bir_asso)
df_ext_rates<-rbind(df_ext_rates,bir_str)
df_ext_rates<-rbind(df_ext_rates,bir_ende)
rx<-(-0.75) yyx<-(-20)
cols<-c(alpha("goldenrod1", 0.5), alpha("royalblue4", 0.5), alpha("darkgreen",0.5), alpha("firebrick",0.5))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
add.pie(z=zs,2006,yyx,radius=log(sum(zs)/1.5)*5,col=cols, labels="",label.dist = 0.5,border=FALSE)
cols<-c(alpha("darkorange1", 0.9), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_str[18])*rx,(bir_str[18])*rx,(tre_str[18])*rx,(amp_str[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[1],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0.9), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_ende[18])*rx,(bir_ende[18])*rx,(tre_ende[18])*rx,(amp_ende[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[1],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0.9), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.1
ze<-c((mam_str[18])*rx,(bir_str[18])*rx,(tre_str[18])*rx,(amp_str[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[2],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0.9), alpha("darkgreen",0), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_ende[18])*rx,(bir_ende[18])*rx,(tre_ende[18])*rx,(amp_ende[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[2],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0), alpha("darkgreen",0.9), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_str[18])*rx,(bir_str[18])*rx,(tre_str[18])*rx,(amp_str[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[3],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0), alpha("black",0.9), alpha("firebrick",0))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_ende[18])*rx,(bir_ende[18])*rx,(tre_ende[18])*rx,(amp_ende[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[3],col=cols, labels="",border=FALSE)
cols<-c(alpha("darkorange1", 0), alpha("royalblue4", 0), alpha("darkgreen",0), alpha("firebrick",0.9))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_str[18])*rx,(bir_str[18])*rx,(tre_str[18])*rx,(amp_str[18])*rx)
ze[is.na(ze)]<-0.001 p<-ze/zs
p[is.na(p)]<-0.1
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[4],col=cols, labels="",border=FALSE)
cols<-c(alpha("black", 0), alpha("black", 0), alpha("black",0), alpha("black",0.9))
zs<-c((mam_asso[18])*rx,(bir_asso[18])*rx,(tre_asso[18])*rx,(amp_asso[18])*rx)
zs[is.na(zs)]<-0.001
ze<-c((mam_ende[18])*rx,(bir_ende[18])*rx,(tre_ende[18])*rx,(amp_ende[18])*rx)
ze[is.na(ze)]<-0.001
p<-ze/zs
p[is.na(p)]<-0.001
add.pie(z=zs,2006,yyx,radius=log((sum(zs))/2)*5*p[4],col=cols, labels="",border=FALSE)
}
cexX<-1
w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
asp <- w/h
#yy=log(2/1.5)*5*rx
#text("-2 sp.",x=2006,y=-18.25,pos=3,col="black",cex=cexX,font=2)
#draw.ellipse(x=2006,y=-yyx,a = yy/asp,b=yy, nv=500,border="white",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(5/1.5)*5*rx
text("-5 sp.",x=2006,y=-16.5,pos=3,col="black",cex=cexX,font=1)
draw.ellipse(x=2006,y=yyx,a = yy/asp,b=yy, nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(10/1.5)*5*rx
text("-10 sp.",x=2006,y=-13.75,pos=3,col="black",cex=cexX,font=1)
draw.ellipse(x=2006,y=yyx,a = yy/asp,b=yy, nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
yy=log(20/1.5)*5*rx
text("-20 sp.",x=2006,y=-11,pos=3,col="black",cex=cexX,font=1)
draw.ellipse(x=2006,y=yyx,a = yy/asp,b=yy, nv=500,border="black",col=NA ,lty=1,density=NULL,lwd=0.5)
par(xpd=NA)
}
dev.off()
# write out extinction rates
cnames<-c(rep("Americas",3*4),rep("Africa",3*4),rep("Asia-Pacific",3*4))
df_ext_rates<-as.data.frame(df_ext_rates)
colnames(df_ext_rates)<-seq(2001,2018,1)
df_ext_rates$continent<-cnames
write.csv(df_ext_rates,file=paste0(OUTPUT,"df_extinction_rates.csv"))
#############################################################################################################################################################################################################################