04_FiguresTables.Rmd

---
title: "Project#02 - Figures and Tables"
author: "Francesco Maria Sabatini"
date: "4/28/2020"
output:
  html_document:
    toc: true
    theme: united
---

<center>
  ![](https://www.idiv.de/fileadmin/content/Files_sDiv/sDiv_Workshops_Photos_Docs/sDiv_WS_Documents_sPlot/splot-long-rgb.png "sPlot Logo")
</center>
  
    
      
        
**Timestamp:** `r date()`  
**Drafted:** Francesco Maria Sabatini  
**Revised:**  
**Version:** 1.0
  
This report documents the creation of figures and tables for the sPlotOpen manuscript.
  
```{r results="hide", message=F, warning=F}
library(tidyverse)
#library(openxlsx)
library(bib2df)
library(knitr)
library(kableExtra)
library(viridis)
library(plotbiomes)
library(cowplot)

library(raster)
library(sp)
library(sf)
library(rgdal)
library(rnaturalearth)
library(dggridR)
# library(rgeos)

library(Taxonstand)

#save temporary files
write("TMPDIR = /data/sPlot/users/Francesco/_tmp", file=file.path(Sys.getenv('TMPDIR'), '.Renviron'))
write("R_USER = /data/sPlot/users/Francesco/_tmp", file=file.path(Sys.getenv('R_USER'), '.Renviron'))
#rasterOptions(tmpdir="/data/sPlot/users/Francesco/_tmp")
```

#Load sPlotOpen data and create spatial objects
```{r}
load(file.path("_sPlotOpenDB", "sPlotOpen.RData"))
#header.oa <- header.oa %>% 
#  filter(!is.na(SoilClim_PC1))
```

Data Preparation for spatial plotting
```{r, cache=T, results="hide", warning=F, message=F}
header.sf <- SpatialPointsDataFrame(coords= header.oa %>% 
                                      dplyr::select(Longitude, Latitude), 
                                    proj4string = CRS("+init=epsg:4326"), 
                                    data=data.frame(PlotObservationID=header.oa$PlotObservationID, 
                                                    Dataset=header.oa$Dataset)) %>% 
  st_as_sf() %>% 
  st_transform(crs = "+proj=eck4")
```

# Load ancillary geographic data and create figure templates

Country boundaries and world graticule
```{r, cache=T, results="hide", warning=F, message=F}
#data downloaded from rnaturalearth package
countries <- readOGR("/data/sPlot/users/Francesco/Ancillary_Data/naturalearth/ne_110m_admin_0_countries.shp") %>% 
  st_as_sf() %>% 
  st_transform(crs = "+proj=eck4") %>% 
  st_geometry()
graticules <- readOGR("/data/sPlot/users/Francesco/Ancillary_Data/naturalearth/ne_110m_graticules_15.shp") %>% 
  st_as_sf() %>% 
  st_transform(crs = "+proj=eck4") %>% 
  st_geometry()

bb <- readOGR("/data/sPlot/users/Francesco/Ancillary_Data/naturalearth/ne_110m_wgs84_bounding_box.shp") %>% 
  st_as_sf() %>% 
  st_transform(crs = "+proj=eck4") %>% 
  st_geometry()
```


Continent boundaries
```{r}
sPDF <- rworldmap::getMap(resolution="coarse")
continent <- sPDF[,"continent"]
crs(continent) <- CRS("+init=epsg:4326")
continent@data[243,"continent"] <- "South America" ## Manually correct missing data
# create clipped version of continent to avoid going beyond 180 lON
coords <- data.frame(x=c(-180,180,180,-180),
                     y=c(-90,-90,90,90))
bboxc = Polygon(coords)
bboxc = SpatialPolygons(list(Polygons(list(bboxc), ID = "a")), proj4string=crs(continent))
continent_clipped <- rgeos::gIntersection(continent[-137,], bboxc, byid=T) # polygon 137 gives problems... workaround
continent_clipped <- continent_clipped %>% 
  st_as_sf()
```

Template of Global map - with country borders
```{r, cache=T, results="hide", warning=F, message=F}
# create ggplot template of the world map
w3a <- ggplot() +
  geom_sf(data = bb, col = "grey20", fill = "white") +
  geom_sf(data = graticules, col = "grey20", lwd = 0.1) +
  geom_sf(data = countries, fill = "grey90", col = NA, lwd = 0.3) +
  coord_sf(crs = "+proj=eck4") +
  theme_minimal() +
  theme(axis.text = element_blank(), 
        legend.title=element_text(size=12), 
        legend.text=element_text(size=12),
        legend.background = element_rect(size=0.1, linetype="solid", colour = 1), 
        legend.key.height = unit(1.1, "cm"), 
        legend.key.width = unit(1.1, "cm")) +
  scale_fill_viridis()
```

Create template of Global Map - without country borders
```{r}
w4a <- ggplot() +
    geom_sf(data = bb, col = "grey20", fill = "white") +
    geom_sf(data = continent_clipped, fill = "grey90", col = NA, lwd = 0.3) +
    geom_sf(data = bb, col = "grey20", fill = NA) +
    #geom_sf(data = graticules, col = "grey20", lwd = 0.1) +
    coord_sf(crs = "+proj=eck4") +
    theme_minimal() +
    theme(axis.text = element_blank(), 
          legend.title=element_text(size=12), 
          legend.text=element_text(size=12),
          legend.background = element_rect(size=0.1, linetype="solid", colour = 1), 
          legend.key.height = unit(1.1, "cm"), 
          legend.key.width = unit(1.1, "cm"))

```



# Figures
Figures and tables for the manuscript in the [sPlotOpen_Manuscript](https://fmsabatini.github.io/sPlotOpen_Manuscript/) project.

## Figure 1 - Geographic distribution of plots

Map of plot distribution - Version 1 - Coloured points.  
Each colour represents a database. Please note there are not enough colours in the palette to represent all 105 datasets
```{r, fig.width=8, fig.height=6, fig.align="center", warning=F, message=F, cache=T}
Figure1a <- w3a + 
  geom_sf(data=header.sf, aes(color=Dataset), pch="+", size=1, alpha=0.8) + # aes(col=Dataset),
  geom_sf(data = countries, col = "grey20", fill=NA, lwd = 0.3) + 
  theme(legend.position = "none")
```

Version 2  - hexagons
```{r, fig.width=8, fig.height=6, fig.align="center", message=F, cache=T}
header2 <- header.oa %>% 
  dplyr::filter(Resample_1) %>% 
  dplyr::select(PlotObservationID, Latitude, Longitude) %>% 
  filter(!(abs(Longitude) >171 & abs(Latitude>70)))
dggs <- dgconstruct(spacing=300, metric=T, resround='down')

#Get the corresponding grid cells for each plot
header2$cell <- dgGEO_to_SEQNUM(dggs, header2$Longitude, header2$Latitude)$seqnum

#Calculate number of plots for each cell
header.dggs   <- header2 %>% 
  group_by(cell) %>% 
  summarise(value.out=log(n(), 10))

#Get the grid cell boundaries for cells 
grid   <- dgcellstogrid(dggs, header.dggs$cell, frame=F) %>%
  st_as_sf() %>% 
  mutate(cell = header.dggs$cell) %>% 
  mutate(value.out=header.dggs$value.out) %>% 
  st_transform("+proj=eck4") %>% 
  st_wrap_dateline(options = c("WRAPDATELINE=YES"))

## plotting
Figure1b <- w3a + 
    geom_sf(data=grid, aes(fill=value.out),lwd=0, alpha=0.9)    +
    geom_sf(data = countries, col = "grey20", fill=NA, lwd = 0.3) + 
    scale_fill_viridis(
      name="# plots", breaks=0:5, labels = c("1", "10", "100",
                                             "1,000", "10,000", "100,000"), option="viridis")

#ggsave("_output/figure1.png", plot=Figure1, width=8, height=4, units="in", dpi=300)
```
Panel with two versions of figure 1
```{r, fig.width=5.5, fig.height=5, fig.align="center", cache=T}
fig1.leg <- get_legend(Figure1b + 
                         guides(fill = guide_colourbar(barwidth = 1, barheight = 6)) +
                         theme(legend.title = element_text(size = 7), 
                         legend.text = element_text(size = 7)))
fig1_panel <- plot_grid(Figure1a, NULL,
                        Figure1b + 
                          theme(legend.position = "none"), fig1.leg,
                        nrow=2, ncol=2, byrow = T, #labels = c("a","", "b", ""),
                        rel_widths = c(0.84,0.16))
ggsave("_output/figure1.pdf", plot=fig1_panel, width=5.5, height=5, units="in", dpi=300)
fig1_panel
```


## Figure 2 - PCA graph + world map of selected plots
Only for Resample draw #1  
Import PCA data
```{r}
### load PCA ordination of the world
load("_data/pca3.RData")
path.sPlot <- "/data/sPlot2.0/"
load(paste(path.sPlot, "splot.world2.RData", sep="/"))
```


```{r}
plot_data <- header.oa %>% 
  dplyr::filter(Resample_1) %>% 
  dplyr::select(PlotObservationID, Longitude, Latitude)
## code adapted from @lenjon's 'resampling_2d_JL.R'
CRSlonlat <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
coords <- cbind(plot_data$Longitude, plot_data$Latitude)
coords <- SpatialPoints(coords, proj4string=CRSlonlat)
plot_data <- SpatialPointsDataFrame(coords, plot_data)#, proj4string=CRSlonlat)


# Create world rasters of PCA values and extract plot values by geographic intersection
# raster at half a degree resolution (cf. 30 arc minute resolution)
rgeo <- raster(nrows=360, ncols=720, xmn=-180, xmx=180, ymn=-90, ymx=90) 
rgeo <- disaggregate(rgeo, fact=12) # raster at 2.5 arc minute resolution
splot.world2$cellID <- cellFromXY(rgeo, cbind(splot.world2$RAST_X, splot.world2$RAST_Y))

### create rasters from PCA
posit <- splot.world2$cellID
temp <- getValues(rgeo)

temp[posit] <- pca3$x[, 1]
PC1_r <- setValues(rgeo, temp)
temp[posit] <- pca3$x[, 2]
PC2_r <- setValues(rgeo, temp)

# Extract pca valus from PCA rasters
plot_data@data$pc1_val <- extract(PC1_r, coordinates(plot_data))
plot_data@data$pc2_val <- extract(PC2_r, coordinates(plot_data)) 


# Compute the density of environmental conditions available at the global scale across the entire bivariate (PC1-PC2) environmental space
res <- 100 # Setting the number of bins per PCA axis to 100
reco <- raster(nrows=res, ncols=res, xmn=min(pca3$x[, 1]), xmx=max(pca3$x[, 1]), 
               ymn=min(pca3$x[, 2]), ymx=max(pca3$x[, 2]))
PC1_PC2_r <- rasterize(pca3$x[, 1:2], reco, fun="count")
plot_data <- plot_data@data
plot_data$pc_cellID <- cellFromXY(reco, cbind(plot_data$pc1_val, plot_data$pc2_val))


# Compute the sampling effort (number of vegetation plots) per environmental unit (cell) across the entire bivariate (PC1-PC2) environmental space
sPlot_reco <- rasterize(plot_data[, c("pc1_val", "pc2_val")], reco, fun="count")
# Put zero values for the empty cells (cf. there is no vegeteation plots available for those environmental conditions: gaps) 
temp1 <- getValues(PC1_PC2_r)
temp1[!is.na(temp1)] <- 0
temp2 <- getValues(sPlot_reco)
temp2[which(temp1==0&is.na(temp2))] <- 0
sPlot_reco <- setValues(reco, temp2)

plot_data <- plot_data %>% 
  rename(PC1=pc1_val, PC2=pc2_val)
```

Transform to tibbles
```{r}
PC1_PC2.tbl <- data.frame(xyFromCell(PC1_PC2_r, cell = 1:10000)) %>% 
  rename(PC1=x, PC2=y) %>% 
  mutate(values=getValues(PC1_PC2_r)) %>% 
  as_tibble()

sPlot_reco.tbl <- data.frame(xyFromCell(sPlot_reco, cell = 1:10000)) %>% 
  rename(PC1=x, PC2=y) %>% 
  mutate(values=getValues(sPlot_reco)) %>% 
  as_tibble()
```



Figure2 - Bottom left.  
Make gridded heatmap plot
```{r}
pca.heatmap <- ggplot() + 
  geom_tile(data=PC1_PC2.tbl %>% 
              dplyr::filter(!is.na(values)), aes(x=PC1, y=PC2), 
            col = gray(0.8),
            fill = gray(0.8)) + 
  geom_tile(data=sPlot_reco.tbl %>% 
              dplyr::filter(!is.na(values)) %>% 
              dplyr::filter(values>0), ### !!!
            aes(x=PC1, y=PC2, col=values, fill=values)) +
  scale_fill_viridis("Number\nof plots", option = "magma") +
  scale_color_viridis("Number\nof plots", option = "magma") +
  theme_bw() + 
  coord_equal() + 
  theme(legend.position = "bottom", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

```

Figure 2 - Bottom right  
Randomly select four pixels with >45 plots, extract center coordinates and add to heatmap.
```{r}
set.seed(558)
ABC <- sPlot_reco.tbl %>% 
  filter(!is.na(values)) %>% 
  filter(values>45 & values<=50) %>% 
  mutate(upper=PC2 > 1) %>% 
  mutate(right=PC1 > 0) %>% 
  group_by(upper, right) %>% 
  sample_n(1) %>% 
  ungroup() %>% 
  mutate(label=c("D", "C", "A", "B"))
  

### add labels to heatmap
pca.heatmap2 <- pca.heatmap + 
  geom_point(data=ABC, aes(x=PC1, y=PC2), pch=21,  size=2, col="black", fill="white") + 
  ggrepel::geom_label_repel(data=ABC, aes(x=PC1, y=PC2, label=label), size = 2.5, fill = alpha(c("white"),0.7)) + 
  theme(legend.position = "left")
```

For each of the selected grid cells, get all plots belonging to that cell, and show their geographical distribution.
```{r}
inset.list <- list()
rangex <- max(PC1_PC2.tbl$PC1) - min(PC1_PC2.tbl$PC1)
rangey <- max(PC1_PC2.tbl$PC2) - min(PC1_PC2.tbl$PC2)
for(i in c("A", "B", "C", "D")){
  tmp.ABC <- ABC %>% 
    filter(label==i)
  xmin <- tmp.ABC$PC1 - 0.5*rangex/100
  xmax <- tmp.ABC$PC1 + 0.5*rangex/100
  ymin <- tmp.ABC$PC2 - 0.5*rangey/100
  ymax <- tmp.ABC$PC2 + 0.5*rangey/100
  coords.ABC <- plot_data %>% 
    mutate(sel=(PC1 > xmin & PC1 < xmax &
                PC2 > ymin & PC2 < ymax)) %>% 
    filter(sel) %>% 
    dplyr::select(Longitude, Latitude) %>% 
    rowwise() %>% 
    mutate_at(.vars=vars(Longitude, Latitude),
              .funs=~jitter(.)) %>%  
    ungroup() %>% 
    SpatialPoints(proj4string = CRS("+init=epsg:4326")) %>% 
    st_as_sf() %>% 
    st_transform(crs = "+proj=eck4")
  inset.list[[i]] <- w4a + 
    geom_sf(data=coords.ABC, col=2, pch="+", size=2)
}
```

Figure 2 - Top  

```{r}
splot.world2.eckert <- SpatialPointsDataFrame(coords=splot.world2 %>% 
                         dplyr::select(RAST_X, RAST_Y), 
                       data = splot.world2 %>% 
                         dplyr::select(RAST_ID), 
                       proj4string = CRS("+init=epsg:4326")) %>% 
  spTransform(CRS("+proj=eck4"))

splot.world2.eckert <- data.frame(splot.world2.eckert@coords, 
                                   splot.world2.eckert@data)

PCA_tbl <- as_tibble(splot.world2.eckert) %>%
  dplyr::select(RAST_ID, RAST_X, RAST_Y) %>%
  left_join(as_tibble(as.data.frame(pca3$x[,1:2]) %>%
              rownames_to_column(var="RAST_ID")) %>%
              mutate(RAST_ID=as.integer(RAST_ID)),
      by="RAST_ID") %>%
    mutate(PC0=1:n())

ggpc1 <- w4a + 
  geom_tile(data=PCA_tbl %>% 
              mutate(PC1=ifelse(PC1> 6,  6, PC1)) %>% 
              mutate(PC1=ifelse(PC1< -6, -6, PC1)), 
            aes(x=RAST_X,y=RAST_Y, fill=PC1, color=PC1)) + 
  geom_sf(data = bb, col = "grey20", fill = NA) +
  scale_fill_distiller("PC1", type = "seq", palette = "Spectral", 
                       direction=-1, limits = c(-6.1,6.1), breaks=seq(-6,6, by=3), labels=c("<-6", -3, 0, 3,">6" )) +
  scale_color_distiller("PC1", type = "seq", palette = "Spectral", 
                        direction=-1, limits = c(-6.1,6.1), breaks=seq(-6,6, by=3), labels=c("<-6", -3, 0, 3,">6" )) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        axis.text = element_blank(),
        panel.border = element_blank(), 
        axis.title = element_blank())



ggpc2 <- w4a + 
  geom_tile(data=PCA_tbl %>% 
              mutate(PC2=ifelse(PC2> 9,  9, PC2)) %>% 
              mutate(PC2=ifelse(PC2< -6, -6, PC2)), 
            aes(x=RAST_X,y=RAST_Y, fill=PC2, color=PC2)) + 
  geom_sf(data = bb, col = "grey20", fill = NA) +
  scale_fill_distiller("PC2", type = "seq", palette = "Spectral", 
                       direction=+1, limits = c(-6.1,9.1), breaks=seq(-6,9, length.out = 6), 
                       labels=c("<-6", -3, 0, 3, 6,">9" )) +
  scale_color_distiller("PC2", type = "seq", palette = "Spectral", 
                       direction=+1, limits = c(-6.1,9.1), breaks=seq(-6,9, length.out = 6), 
                       labels=c("<-6", -3, 0, 3, 6,">9" )) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        axis.text = element_blank(),
        panel.border = element_blank(), 
        axis.title = element_blank())
```


Figure 2 - Panel
```{r, fig.width=8, fig.height=5.6, fig.align="center", cache=T, cache.lazy = FALSE}
library(cowplot)
varexpl <- round(pca3$sdev^2/sum(pca3$sdev^2)*100,1)
panel.out1 <- 
  plot_grid(plot_grid(NULL, 
                      ggpc1 + theme(plot.margin=margin(r=.5, unit="cm")), 
                      ggpc2 + theme(plot.margin=margin(r=.5, unit="cm")), 
                      nrow=1, rel_widths = c(0.13,1,1)),
            plot_grid(
                pca.heatmap2 + 
                theme(legend.position=c(0.13, 0.68), 
                      legend.background = element_blank()) +
                xlab(paste("PC1 (", varexpl[1], "%) - cold/seasonal to hot/stable", sep="")) + 
                ylab(paste("PC2 (", varexpl[2], "%) - dry to wet", sep="")) + 
                xlim(c(-15,15)) + ylim(c(-5,22)),  
              plot_grid(plot_grid(inset.list[[1]], inset.list[[4]], NULL, 
                                  nrow=3, labels = c("A", "D", ""), rel_heights = c(1,1,0.1)),
                        plot_grid(inset.list[[2]], inset.list[[3]], NULL, 
                                  nrow=3, labels = c("B", "C", ""), rel_heights = c(1,1,0.1)))), 
            nrow=2, rel_heights = c(1,1.5), align="v")
          

ggsave(filename = "_output/figure2.tiff", width=8, height=5.6, units = "in", dpi=300, plot=panel.out1)

panel.out1
```


## Figure 3 - Whittaker Biome Graph
Get climatic data
```{r}
load("/data/sPlot/releases/sPlot2.1/sPlot_header_chelsa_20161124.RData")
climate.oa <- climate %>% 
  filter(PlotID %in% (header.oa %>% 
                        filter(Resample_1) %>% 
                        pull(PlotObservationID))) %>% 
  dplyr::select(-POINT_X, -POINT_Y) %>% 
  rename(PlotObservationID=PlotID)
```

Figure 3 - Left
Create plot of Schultz' biomes
```{r}
biome.order <- c('Polar and subpolar zone' ,'Alpine' ,'Boreal zone' ,'Temperate midlatitudes' ,
                 'Dry midlatitudes' ,'Dry tropics and subtropics' ,'Subtropics with year-round rain' ,
                 'Subtropics with winter rain' ,'Tropics with summer rain' ,'Tropics with year-round rain')
biome.labs <- c('Polar & subpolar' ,'Alpine' ,'Boreal zone' ,'Temperate midlatitudes' ,
                 'Dry midlatitudes' ,'Dry tropics & subtropics' ,'Subtropics - year-round\n rain' ,
                 'Subtropics - winter\n rain' ,'Tropics - summer rain' ,'Tropics - year-round\n rain')


mypalette <- palette(c('#CAB2D6','#6A3D9A', #violets
                        '#A6CEE3','#1F78B4', #blues
                        '#FDBF6F','#FF7F00', #orange
                        '#B2DF8A','#33A02C', #greens
                        '#FB9A99','#E31A1C' #reds
                        ))

# Plot of Temp vs Prec + sBiomes
biome.schu <- ggplot() + 
  theme_classic() +
  #geom_path(aes(x, y), data=contour_95) +
  geom_point(data=climate.oa %>% 
               left_join(header.oa %>% 
                           dplyr::filter(Resample_1) %>% 
                           dplyr::select(PlotObservationID, Biome), 
                         by="PlotObservationID") %>%  
              filter(bio12<4500 & bio01>-11) %>% 
               mutate(Biome=factor(Biome, levels=biome.order, labels=biome.labs))
              # mutate(alpha=ifelse(Biome=="Tropics with summer rain", 1, 1/3))
               , #filter out for plotting reasons 
              aes(x=bio01, y=bio12, col=Biome), 
             alpha=1/3, 
             cex=1/15) + 
  xlab("Temperature (°C)") + 
  ylab("Precipitation (mm)") + 
  scale_color_manual(values = c('#CAB2D6','#6A3D9A', #violets
                        '#A6CEE3','#1F78B4', #blues
                        '#FDBF6F','#FF7F00', #orange
                        '#B2DF8A','#33A02C', #greens
                        '#FB9A99','#E31A1C' #reds
                        ), name="sBiomes") + 
  guides(color = guide_legend(override.aes = list(size=5, shape=15, alpha=1))) + 
  theme(#plot.margin=margin(r=-0.1, unit="cm"),
        axis.title = element_text(size=9),
        axis.text = element_text(size=9),
        legend.text = element_text(size=8), 
        legend.title = element_text(size=9))
```

Figure 3 - Right  
Create Whittaker plot
```{r}
whitt.biome <- whittaker_base_plot() +
  theme_classic() + 
  geom_point(data=climate.oa %>% 
                filter(bio12<4500 & bio01>-11), #filter out for plotting reasons 
              aes(x=bio01, y=bio12/10), 
             alpha=1/4, 
             cex=1/20) + 
    theme(axis.text.y = element_blank(),
          axis.text = element_text(size=9),
          axis.title.y = element_blank(),
          axis.title = element_text(size=9),
          plot.margin=margin(l=-0.2, r=-0.1, unit="cm"), 
          legend.text = element_text(size=8),
          legend.title = element_text(size=9)) 
```

Figure 3 - Panel
```{r, fig.width=8, fig.height=4, fig.align="center", message=F, cache=T}
#Make a panel with the two plots
panel.biomes <- cowplot::plot_grid( biome.schu, whitt.biome, nrow=1, 
                                    align = "h", rel_widths = c(0.5, .5))
ggsave(filename="_output/figure3.tiff", plot = panel.biomes, width = 10, height=4, units="in", dpi=300)
panel.biomes
```

# Supplementary figures



## Figure S1 - Biplot
```{r, fig.width=5.7, fig.height=7, fig.align="center", warning=F}
library(ggrepel)
varexpl <- round(pca3$sdev^2/sum(pca3$sdev^2)*100,1)

mydata <- pca3$x %>% 
         as_tibble() %>% 
         dplyr::select(1:2) %>% 
         mutate(count=1)
myarrows <- pca3$rotation[,1:2] %>% 
  as.data.frame() %>% 
  rownames_to_column("mylab") %>%
  mutate_at(.vars=vars(-mylab), .funs = list(~.*25)) %>% 
  filter(!mylab %in% c("T_ANN", "P_ANN"))
  

ggpca3 <- ggplot(data=mydata) +
    stat_summary_2d(
    aes(x = PC1, y = PC2, z = count),
    bins = 100,
    fun = function(x) {log10(sum(x)+1)},
    alpha=1/2
    ) +
  geom_segment(data=myarrows, 
               aes(x=0, xend=PC1, y=0, yend=PC2), 
               arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
    geom_label_repel(data=myarrows, 
             aes(x=PC1, y=PC2, label=mylab), size=2, 
             position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) + 
  scale_fill_viridis("Number of 2.5 arcmin\n terrestrial grid cells", option = "magma", alpha=1/2,
                     limits=c(0,5), breaks=log10(c(9,99,999,9999)), labels=10^(1:4)) +
  theme_bw() + 
  coord_equal() + 
  theme(legend.position = c(0.2, 0.85), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  #guides(fill = guide_legend(override.aes = list(alpha=1/3))) + 
  xlab(paste("PC1 (", varexpl[1], "%)\n (cold and seasonal to hot and stable)", sep="")) + 
  ylab(paste("PC2 (", varexpl[2], "%)\n (dry to wet)", sep="")) + 
  ylim(c(-6,22))

ggsave("_output/figureS1.pdf", plot = ggpca3, width=5.7, height=7, unit="in", dpi=300)
ggpca3
```


# Tables 
## Table 1 - Database level information
Import databases and create reference tags
```{r, message=F, warning=F}
#load(file.path("_sPlotOpenDB", "sPlotOpen.RData"))
#Import BibTex
bib.db <- bib2df("/data/sPlot/users/Francesco/_sPlot_Management/Consortium/sPlot_References.bib")
#Import database-level information
databases <- read_csv("/data/sPlot/users/Francesco/_sPlot_Management/Consortium/Databases.out.csv")

# create citation tags that can be picked up by Manubot
databases <- databases %>% 
  left_join(bib.db %>% 
              dplyr::select(BIBTEXKEY, DOI, URL), 
            by="BIBTEXKEY") %>% 
  mutate(tag=NA) %>% 
  rowwise() %>% 
  mutate(tag=ifelse(!is.na(DOI), 
                    paste0("@doi:", DOI), 
                    tag)) %>% 
  mutate(tag=ifelse( (is.na(tag) & `GIVD ID` %in% unique(header.oa$GIVD_ID) & !is.na(Citation)), 
                     paste0("@", word(Citation, 1)), 
                     tag)) %>% 
  dplyr::select(-DOI, -URL, -BIBTEXKEY)
  

```

Create Table 1

```{r}
table1 <- databases %>% 
  filter(`Still in sPlot`==T, 
         Via!="Aggregator") %>% 
  dplyr::select(-Via, -`Still in sPlot`, -label) %>% 
  distinct() %>% 
  left_join(header.oa %>% 
              group_by(GIVD_ID) %>% 
              summarize(contributed_plots=n(), .groups = 'drop'), 
            by=c("GIVD ID"="GIVD_ID")) %>% 
  filter(!is.na(contributed_plots)) %>% 
  replace_na(list(tag="",
                  `Deputy custodian`="")) %>% 
  dplyr::select(`GIVD ID`, `Dataset name`=`DB_name GIVD`,  Custodian, `Deputy custodian`, `Nr. open-access plots` = contributed_plots, Ref=tag) %>% 
  arrange(`GIVD ID`) %>% 
  # replace citation rendered wrongly using doi
  mutate(Ref=replace(Ref,
                     list=`GIVD ID`=="AU-AU-002", 
                     values="@isbn:9781315368252"))


write_csv(table1, "_output/Table1_Databases.csv")
```

```{r, echo=F}
knitr::kable(table1%>% 
               slice(1:20),
             caption="Table 1 - Database level information [only first 20 rows shown]") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
                  full_width = F, position = "center")
```


## Table 2 - Metadata and environmental data included in header  
Define unit of measurements for columns
```{r}
um <- c('Latitude'='° (WGS84)',
'Longitude'='° (WGS84)',
'Location_uncertainty'='m',
'Releve_area'='m^2^',
'Elevation'='m a.s.l.',
'Aspect'='°',
'Slope'='°',
'Date_of_recording' = 'dd-mm-yyyy',
'Cover_total'='%',
'Cover_tree_layer'='%',
'Cover_shrub_layer'='%',
'Cover_herb_layer'='%',
'Cover_moss_layer'='%',
'Cover_lichen_layer'='%',
'Cover_algae_layer'='%',
'Cover_litter_layer'='%',
'Cover_bare_rocks'='%',
'Cover_cryptogams'='%',
'Cover_bare_soil'='%',
'Height_trees_highest'='m',
'Height_trees_lowest'='m',
'Height_shrubs_highest'='m',
'Height_shrubs_lowest'='m',
'Height_herbs_average'='cm',
'Height_herbs_lowest'='cm',
'Height_herbs_highest' = 'cm')
um <- data.frame(Variable=names(um), `Unit of Measurement`=um)

```

Create table 2 
```{r}
table2 <- header.oa %>% 
  dplyr::summarize_at(.vars=vars(!starts_with("PlotObservationID")),
               .funs = list(xxxNo.records=~sum(!is.na(.)),
                            xxxType.of.variable=~ifelse("logical" %in% class(.), "b",
                                           ifelse("ordered" %in% class(.), 
                                                 "o", 
                                                  ifelse(any(class(.) %in% c("character", "factor")), 
                                                         "n",
                                                         ifelse(class(.)=="Date",
                                                                "d",
                                                                "q")))),
                            xxxLevels=~(ifelse(is.numeric(.)|lubridate::is.Date(.), 
                                               paste(range(., na.rm=T), collapse=" - "),
                                               ifelse(is.ordered(.), 
                                                      paste(paste(1:nlevels(.), 
                                                                  levels(.), sep=" = "), collapse=", "),
                                                      ifelse(is.factor(.), 
                                                              paste(levels(.), collapse=", "),
                                                              ifelse(is.logical(.),
                                                                     paste(names(table(.)), "=",table(.),
                                                                           collapse="; "),
                                                                     ""))))))) %>%
  gather(key="Variable") %>% 
  separate(Variable, into = c("Variable", "feature"), sep="_xxx") %>% 
  spread(key=feature, value = value) %>% 
  rename(`Range/Levels`=Levels) %>% 
  mutate(Variable=factor(Variable, levels=colnames(header.oa))) %>% 
  arrange(Variable) %>% 
  left_join(um, by="Variable") %>% 
  mutate(Unit.of.Measurement=as.character(Unit.of.Measurement)) %>% 
  replace_na(list(Unit.of.Measurement="")) %>% 
  dplyr::select(Variable, `Range/Levels`, 
                `Unit of Measurement`=Unit.of.Measurement, 
                `Nr. of plots with information`=No.records, `Type`=Type.of.variable)

write_csv(table2, "_output/Table2_header.csv")
```

```{r, echo=F}
knitr::kable(table2,
             caption="Table 2 - Variables in header") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
                  full_width = F, position = "center")
```


## Sink tables for Manubot
```{r}
out.file <- "_output/91.ManubotTables.md"
readr::write_lines("## Supplementary Material {.page_break_before}\n", file = out.file)
readr::write_lines("Table: List of databases contributing to sPlotOpen, the environmentally-balanced, open-access, global dataset of vegetation plots. Databases are ordered based on their ID in the Global Index of Vegetation Databases (GVID ID). {#tbl:Table1 tag='1'}\n", file = out.file, append=T)
readr::write_lines("\n \n", file = out.file, append=T)
kable1 <- kable(table1, format = "markdown")
## fix header table 1
kable1[2] <- "|:------------|:--------------------------------------------|:--------------------|:--------------------|--------:|:--------|"
readr::write_lines(kable1, file = out.file, append=T)
readr::write_lines("\n \n \n", file = out.file, append=T)
readr::write_lines("Table: Description of the variables contained in the ‘header’ matrix, together with their range (if numeric) or possible levels (if nominal or binary) and the number of non-empty (i.e., non NA) records. Variable types can be n - nominal (i.e., qualitative variable), o - ordinal, q - quantitative, or b - binary (i.e., boolean), or d - date {#tbl:Table2 tag='2'}. Additional details on the variables are in Bruelheide et al. (2019) [@doi:10.1111/jvs.12710]. GIVD codes derive from Dengler et al. (2011) [@doi:10.1111/j.1654-1103.2011.01265.x]. Biomes refer to Schultz 2005 [@doi:10.1007/3-540-28527-x], modified to include also the world mountain regions by Körner et al. (2017)[@doi:10.1007/s00035-016-0182-6]. The column ESY refers to the EUNIS Habitat Classification Expert system described in Chytrý et al. (2020) [@doi:10.1111/avsc.12519].\n", file = out.file, append=T)
kable2 <- kable(table2, format = "markdown")
## fix header table 2
kable2[2] <- "|:---------------------------|:-------------------------------------------------------------|:-------------------|:-----------|:-----|"
readr::write_lines(kable2, file = out.file, append=T)
```


# SessionInfo
```{r}
sessionInfo()
```