GEDI_PA_figs.Rmd

---
title: "GEDI PA Figures"
output: html_notebook
---

This notebook produces the figures for the Paper 'Global protected areas are effective for preserving forest structure and carbon'.

Set up environment
```{r}
#load any packages necessary
library(dplyr)
library(ggplot2)
local_fig_path <- '/gpfs/data1/duncansongp/GEDI_PA_paper/figures/'

#function for plotting up multiple plots with a single legend
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {

  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)

  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))

  grid.newpage()
  grid.draw(combined)

  # return gtable invisibly
  invisible(combined)

}
```

Fig 1: Methods flowchart
```{r}
#Created Externally by V Leitold
```

Table 1: Input data
```{r}
#Created Externally by M Liang
```

Fig 2: Global GEDI height, RH50, cover & biomass
Code written by M Liang
```{r}
source("/gpfs/data1/duncansongp/amberliang/trends.Earth/git/GEDI_PA/fig2_script.R")
```

Fig 3: Waveform Examples
Created externally by L Duncanson
```{r}
```

Fig 4: Global PAs for paper
Created externally by V Leitold
```{r}
#Created externally by V Leitold
```

Fig 5: Number of GEDI shots per 1 km pixel
Coded by M Liang.
```{r}
#fig 5 a and 5b
#Coded  by M Liang
#GEDI shots distribution per 1km cell  Figure 4a
fpath <- "/gpfs/data1/duncansongp/GEDI_global_PA/WDPA_GEDI_extract3/"
dir(fpath)

shotInCell <- data.frame()
shotf=list.files(paste(fpath,"cell_stats",sep=""), full.names = TRUE)
# vars=c("x","y","gshot","status","pid","iso3")

for (d in shotf){
  print(d)
  cout <- read.csv(d)
  cout <- cout[!is.na(cout$pid),] 
  iso3 <- cout$country %>% unique() %>% na.omit()
  print(iso3)
  tmp <-  data.frame(iso3=cout$country, status=as.factor(cout$status), gshot=cout$gshot_counts, stringsAsFactors=F)
  shotInCell <- shotInCell %>% rbind(tmp, stringsAsFactors=F)
}

print(dim(shotInCell))

# plot
p <- shotInCell %>%
  ggplot( aes(x=gshot)) +
  geom_histogram( binwidth=5, fill="#69b3a2", color="#e9ecef", alpha=0.9) +
  labs(title = "GEDI shots in each matched 1km cell")+
       # ,subtitle = paste("Binwidth = 5; n =", format(nrow(shotInCell), big.mark=","), "; ","iso3 =", length(unique(shotInCell$iso3)), sep="")) +
  labs(x = "Number of GEDI shots in 1km Pixel")+theme_bw()

pzoom <- shotInCell %>%  #zoom in plot
  ggplot( aes(x=gshot)) +
  geom_histogram( binwidth=5, fill="#69b3a2", color="#e9ecef", alpha=0.9) +coord_cartesian(xlim=c(0,100))+
  labs(x = "Number of GEDI shots in 1km Pixel")+theme_bw()


#plot for shot in cell with inset map
p1I <- p + annotation_custom(ggplotGrob(pzoom), xmin = 280, xmax = 820, 
                       ymin = 250000, ymax =820000)+
  theme(text=element_text(family="Times", face="bold", size=14),
        legend.title=element_text(size=12), 
        legend.text=element_text( size=12))+
  labs(tag = "A")

ggplot2::ggsave(filename=paste('/gpfs/data1/duncansongp/GEDI_global_PA/figures/fig4a.png', sep=''), plot=p1I, width=6, height=6, units = "in", bg = "transparent")


#GEDI shots distribution per pa  Figure 4b
shotInPA <- data.frame()
shotpaf=list.files(paste(fpath, "pa_stats",sep=""), full.names = TRUE)

for (d in shotpaf){
  print(d)
  cout <- read.csv(d,header = TRUE)
  cout <- cout[!is.na(cout$pa_id),] 
  iso3 <- cout$iso3 %>% unique() %>% na.omit()
  print(iso3)
  tmp <-  data.frame(iso3=cout$iso3,pid=cout$pa_id, gshot=cout$count_1, stringsAsFactors=F)
  shotInPA <- shotInPA %>% rbind(tmp, stringsAsFactors=F)
}  
print(dim(shotInPA))
# plot
p2 <- shotInPA %>%
  ggplot( aes(x=gshot)) +
  geom_histogram( binwidth=500, fill="#FA8072", color="#FA8072", alpha=0.6) +
  labs(title = "GEDI shots in each matched PA")+
       # ,subtitle = paste("Binwidth = 500; ","n =", format(nrow(shotInPA), big.mark=","), "; iso3 =", length(unique(shotInPA$iso3)), sep=" ")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=15)
  )+
  labs(x = "Number of GEDI shots in each PA")

p2zoom <- shotInPA %>%
  ggplot( aes(x=gshot)) +
  geom_histogram( binwidth=500, fill="#FA8072", color="#FA8072", alpha=0.6)+coord_cartesian(xlim=c(0,20000))+
  labs(x = "GEDI shots in each matched PA")+theme_bw()

#plot for shot in cell with inset map
p2I <- p2 + annotation_custom(ggplotGrob(p2zoom), ymin = 2800, ymax =10000, 
                      xmin = 450000, xmax =1400000)+
  theme(text=element_text(family="Times", face="bold", size=14),
        legend.title=element_text(size=12), 
        legend.text=element_text( size=12)) +
  labs(tag = "B")

ggplot2::ggsave(filename=paste('/gpfs/data1/duncansongp/GEDI_global_PA/figures/fig4b.png', sep=''), plot=p2I, width=6, height=6, units = "in", bg = "transparent")
```

Fig 6: GEDI coverage over a single PA
Created externally by V Leitold
```{r}
#Created externally by V Leitold
```

Fig 7: Distributions of AGBD per biome as an overlapping histogram
```{r}
library(dplyr)
library(tidyr)
library(gridExtra)
library(cowplot)

fpath <- "/gpfs/data1/duncansongp/GEDI_global_PA/WDPA_GEDI_extract3/"
dir(fpath)

sub_for_iso <- readRDS("/gpfs/data1/duncansongp/GEDI_global_PA/figures/all_fullds_2_9_fore_iso.rds")


```

Old Fig 7 - Comparison of four metrics by PFT
Coded by M Liang
```{r}
source("/gpfs/data1/duncansongp/amberliang/trends.Earth/git/GEDI_PA/fig5_script.R")
```

Fig 8: Comparison of four metrics for globe, divided by continent
Note: need to add whiskers for the SD 
```{r}
fpath <- "/Volumes/gsdata1/duncansongp/GEDI_global_PA/WDPA_GEDI_extract3/iso_full_nodup"
dir(fpath)

library(dplyr)
library(tidyr)
library(gridExtra)
library(cowplot)

#read in the full dataset
sub <-readRDS("/gpfs/data1/duncansongp/GEDI_global_PA/figures/all_iso_fullds.rds")

#change pa_status to text
status_factor <- rep('temp', nrow(metric_sub))
status_factor[which(metric_sub$pa_stat==0)] <- 'Control'
status_factor[which(metric_sub$pa_stat==1)] <- 'Protected Area'
status_factor <- as.factor(status_factor)
metric_sub <- cbind(metric_sub, status_factor)

#get the mean difference in the status groups by region.
#AGBD by Continent
sub$region <- fct_explicit_na(sub$region)
sub %>%
  group_by(status, region) %>%
  dplyr::summarise(mean = mean(AGBD, na.rm=TRUE)) -> regional_means_AGBD

mean_control <- regional_means_AGBD[which(regional_means_AGBD$status==0),]
mean_treatment <- regional_means_AGBD[which(regional_means_AGBD$status==1),]
mean_diffs <- 100*((mean_treatment$mean - mean_control$mean) / mean_treatment$mean)

mp <- data.frame(mean_diffs=mean_diffs[-length(mean_diffs)], continent=mean_control$region[-length(mean_diffs)])

mp$colour <- ifelse(mp$mean_diffs < 0, "firebrick1","steelblue")
mp$hjust <- ifelse(mp$mean_diffs > 0, 1.3, -0.3)

agbd_continent <- ggplot(mp, aes(y=mean_diffs, x=continent, hjust = hjust)) + 
  geom_bar(stat = "identity", aes(fill = colour)) +
  labs(y = 'Mean AGBD Difference (% Mg/ha)') +
  labs(x = 'Continent') +
  theme(legend.position = 'none', text = element_text(size=40))

#Height by Continent
sub$region <- fct_explicit_na(sub$region)
sub %>%
  group_by(status, region) %>%
  dplyr::summarise(mean = mean(rh_098, na.rm=TRUE)) -> regional_means_ht

mean_control <- regional_means_ht[which(regional_means_ht$status==0),]
mean_treatment <- regional_means_ht[which(regional_means_ht$status==1),]
mean_diffs <- mean_treatment$mean - mean_control$mean

mp <- data.frame(mean_diffs=mean_diffs[-length(mean_diffs)], continent=mean_control$region[-length(mean_diffs)])

mp$colour <- ifelse(mp$mean_diffs < 0, "firebrick1","steelblue")
mp$hjust <- ifelse(mp$mean_diffs > 0, 1.3, -0.3)

ht_continent <- ggplot(mp, aes(y=mean_diffs, x=continent, hjust = hjust)) + 
  geom_bar(stat = "identity") +
  labs(y = 'Mean Height Difference (m)') +
  labs(x = 'Continent') +
  theme(legend.position = 'none', text = element_text(size=40))

#PAI by Continent
sub$region <- fct_explicit_na(sub$region)
sub %>%
  group_by(status, region) %>%
  dplyr::summarise(mean = mean(pai, na.rm=TRUE)) -> regional_means_pai

mean_control <- regional_means_pai[which(regional_means_pai$status==0),]
mean_treatment <- regional_means_pai[which(regional_means_pai$status==1),]
mean_diffs <- 100*(mean_treatment$mean - mean_control$mean) / mean_treatment$mean

mp <- data.frame(mean_diffs=mean_diffs[-length(mean_diffs)], continent=mean_control$region[-length(mean_diffs)])

mp$colour <- ifelse(mp$mean_diffs < 0, "steelblue","firebrickred1")
mp$hjust <- ifelse(mp$mean_diffs > 0, 1.3, -0.3)

pai_continent <- ggplot(mp, aes(y=mean_diffs, x=continent, hjust = hjust)) + 
  geom_bar(stat = "identity") +
  labs(y = 'Mean PAI Difference (%)') +
  labs(x = 'Continent') +
  theme(legend.position = 'none', text = element_text(size=40))

#Cover by Continent
sub$region <- fct_explicit_na(sub$region)
sub %>%
  group_by(status, region) %>%
  dplyr::summarise(mean = mean(cover, na.rm=TRUE)) -> regional_means_cover

mean_control <- regional_means_cover[which(regional_means_cover$status==0),]
mean_treatment <- regional_means_cover[which(regional_means_cover$status==1),]
mean_diffs <- 100*(mean_treatment$mean - mean_control$mean) / mean_treatment$mean

mp <- data.frame(mean_diffs=mean_diffs[-length(mean_diffs)], continent=mean_control$region[-length(mean_diffs)])

mp$colour <- ifelse(mp$mean_diffs < 0, "firebrickred1","steelblue")
mp$hjust <- ifelse(mp$mean_diffs > 0, 1.3, -0.3)

cover_continent <- ggplot(mp, aes(y=mean_diffs, x=continent, hjust = hjust)) + 
  geom_bar(stat = "identity") +
  labs(y = 'Mean Cover Difference (%)') +
  labs(x = 'Continent') +
  theme(legend.position = 'none', text = element_text(size=40))

png_file<-paste(filename="/gpfs/data1/duncansongp/GEDI_global_PA/figures/fig6.png")
png(width=1800,height=2000,file=png_file)
grid.arrange(ht_continent, cover_continent, pai_continent, agbd_continent, ncol=2)
dev.off()


```

Fig 9: AGBD PA vs. Control overlapping histograms per biome
Coded by M Liang
```{r}

```

Fig XX: Relative importance of variables to determine biomass preservation efficacy 
```{r}

```

Tables and final statistical moments aggregation
```{r}
source("/gpfs/data1/duncansongp/amberliang/trends.Earth/git/GEDI_PA/final_stats_agg.R")

```