# install.packages(c('rpostgis','exactextractr','randomForest'))
library(sf)
library(RPostgreSQL)
library(rpostgis)
library(raster)
library(magrittr)
library(exactextractr)
library(data.table)
library(ggplot2)
library(randomForest)
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Loading required package: DBI
Loading required package: sp
Attaching package: ‘magrittr’
The following object is masked from ‘package:raster’:
extract
Attaching package: ‘data.table’
The following object is masked from ‘package:raster’:
shift
randomForest 4.6-14
Type rfNews() to see new features/changes/bug fixes.
Attaching package: ‘randomForest’
The following object is masked from ‘package:ggplot2’:
margin
#Load feature dataset
#Extended project boundary of a region in Panama that might be interesting as a re-forestation project
# eco.sf <- st_read('/work/input_data/EcoVanao_ExtendedAOI.shp')
pnw.stack <- stack('/work/input_data/PNW_Disturbance_Analysis.tif')
# terrain.stack <- list.files('/work/input_data/B17/','*.tif', full.names = TRUE) %>% stack
biomass.stack <- list.files('/work/input_data/Global_Maps_C_Density_2010_1763/data/','*.tif', full.names = TRUE) %>% stack
#crop and resample terrain and reference biomass
# terrain.stack <- crop(terrain.stack,pnw.stack)
# terrain.stack <- resample(terrain.stack,eco.stack,method = 'ngb')
# biomass.stack <- crop(biomass.stack,pnw.stack)
# biomass.stack <- resample(biomass.stack,pnw.stack,method = 'ngb')
#Stack raster data
# reference.stack <- stack(pnw.stack,biomass.stack)
#For security reasons, this section of the script will be redacted.
# conn <- RPostgreSQL::dbConnect("PostgreSQL",
# host = '***',
# dbname = "gedi-research",
# user = "postgres",
# password = "***")
# and get data
# panama.sf <- st_read(conn, "gedi_2a_data")
# panama_2b.sf <- st_read(conn, "gedi_2b_data")
#For security purposes, we'll work off static versions of the above data
# panama.sf <- st_read('/work/gedi_2a_data.geojson')
pnw_2b.sf <- st_read('/work/input_data/pnw_gedi_2b_data.geojson')
Reading layer `pnw_gedi_2b_data' from data source `/work/input_data/pnw_gedi_2b_data.geojson' using driver `GeoJSON'
Simple feature collection with 33035 features and 24 fields
geometry type: POINT
dimension: XY
bbox: xmin: -122.4007 ymin: 44.02857 xmax: -122.2453 ymax: 44.13755
geographic CRS: WGS 84
# #drop cover and pai into the 2a data
# panama.sf$cover <- panama_2b.sf$cover[match(panama_2b.sf$shot_number,panama.sf$shot_number)]
# panama.sf$pai <- panama_2b.sf$pai[match(panama_2b.sf$shot_number,panama.sf$shot_number)]
#Pull RH_101
# panama.sf$RH_101 <- panama.sf$rh %>% lapply( ., function(x){
# (x %>% gsub('\\{','',.) %>% gsub('\\}','',.) %>% strsplit(.,','))[[1]][[101]] %>% as.numeric
# }) %>% unlist
# #Pull RH_50
# panama.sf$RH_50 <- panama.sf$rh %>% lapply( ., function(x){
# (x %>% gsub('\\{','',.) %>% gsub('\\}','',.) %>% strsplit(.,','))[[1]][[50]] %>% as.numeric
# }) %>% unlist
#Extract all raster bands, add to data.table, convert back to sf
pnw_2b.sf <- data.table(pnw_2b.sf, raster::extract(pnw.stack,pnw_2b.sf) %>% data.table) %>% st_as_sf
#make yod consistent since some come in NaN
names(pnw_2b.sf)
#to dt
pnw_2b.dt <- pnw_2b.sf[,] %>% data.table
#We'll need to calculate an upper bound on what we think is possible regarding forest cover.
#This can be thought of as 'If left to its own devices, what is the local maxima we can expect to acheive in regards to cover?'
#This will represent the upper bounds that we expect regions to be able to recover to.
pnw_2b.dt[cover != -9999,][['cover']] %>% qplot()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Looking at the distribution of cover, we can see a pretty clean break at around ~65% cover.
#Here we'll interpret this as the lower bound of undistiturbed regions of forest
#Working backwards we can see that this corresponds appoximately to the 95th percentile of cover data.
#This implies that regionally, around 5% of the spaces sampled were in an 'undisturbed state.
print(quantile(pnw_2b.dt[cover != -9999,][['cover']],.95))
asymtote.param <- pnw_2b.dt[cover != -9999 & cover > quantile(cover,.95),list(cover = mean(cover),sd = sd(cover)),]
print(asymtote.param)
95%
0.9274914
cover sd
1: 0.9253503 0.02555526
pnw_2b.dt[cover != -9999 & yod != 0,][['yod']] %>% table
pnw_2b.dt[cover != -9999 & yod != 0,][['mag']] %>% qplot
pnw_2b.dt[cover != -9999 & yod != 0,][['yod']] %>% qplot
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Here, we use change in cover over time as a proxy for growth after catastrophic removal.
#Remove any bad cover, no observed destruction, or bad quality flag
#We might be assuming an amount of catastrophic deforestation that isnt warranted.
#We know a disturbance happened; we dont know how compeltely.
# pnw_2b.dt[cover != -9999 & yod != 0,][['yod']] %>% qplot
pnw_2b.dt$yad <- 2021 - pnw_2b.dt$yod
cover.lm <- lm(cover~yad,pnw_2b.dt[cover != -9999& yod !=0 & mag >250,])
years_to_recovery <- sqrt((predict(cover.lm,data.table(yad =1:100)) - 0.9253503)^2) %>% which.min
print('Estimated time to recover: 81')
ggplot(pnw_2b.dt[cover != -9999& yod !=0 & mag >250,
list(cover = mean(cover),
sd = sd(cover),
per_10 = quantile(cover,.1),
per_90 = quantile(cover,.90)),by = yod],
#Years since disturbance
aes(x=2021-yod,y=cover))+
ylab('biomass as % of maximum')+
xlab('years after disturbance')+
geom_point()+
# geom_errorbar(aes(ymin = cover-sd,ymax = cover+sd))+
geom_segment(aes(x = 70, y = asymtote.param$cover, xend = 100, yend = asymtote.param$cover))+
geom_segment(aes(x = 70, y = asymtote.param$cover - asymtote.param$sd , xend = 100, yend = asymtote.param$cover - asymtote.param$sd ),linetype = 'dashed')+
geom_segment(aes(x = 70, y = asymtote.param$cover + asymtote.param$sd , xend = 100, yend = asymtote.param$cover + asymtote.param$sd ),linetype = 'dashed')+
geom_smooth(method = 'lm')+
geom_abline(slope = 0.01017 , intercept = 0.10242)
`geom_smooth()` using formula 'y ~ x'
KernelInterrupted: Execution interrupted by the Jupyter kernel.