In RETURN-project/GEE.aux: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE, eval = F)

inputs

ifolder <- '/home/wanda/Documents/data/upscaleRecovery/test/sample_n10000'#'../inst/testdata/' # folder with all processed data
ofolder <- '/home/wanda/Documents/data/upscaleRecovery/test/sample_n10000'#'../inst/testdata/' # folder with all processed data
ifile <- 'samples_forest_fire_nat_n10000_' # base name of the input file

import data

load(file.path(ofolder,paste0('df_rec_tot')))
load(file.path(ofolder,paste0('df_out_tot')))
load(file.path(ofolder,paste0('fireDts_tot')))
load(file.path(ofolder,paste0('seg_tot')))

dts <- as.Date(names(df_out_tot)[-c(1:4)])

Spatial plots

# Spatial plot
library(GISTools)
library(stringr)
library(sf)
library(raster)

# df_rec_tot <- dplyr::filter(df_rec_tot,ecoReg != 599 & ecoReg != 587 & ecoReg != 588 & ecoReg != 589 & ecoReg != 567 & ecoReg != 596 & ecoReg != 590 & ecoReg != 591 & ecoReg != 569 & ecoReg != 611 & ecoReg != 575 & ecoReg != 602 & ecoReg != 582 & ecoReg != 570 & ecoReg != 592 & ecoReg != 571 & ecoReg != 576 & ecoReg != 603 & ecoReg != 572 & ecoReg != 577 & ecoReg != 597 & ecoReg != 613 & ecoReg != 598 & ecoReg != 593 & ecoReg != 583 & ecoReg != 584 & ecoReg != 606 & ecoReg != 585 & ecoReg != 578 & ecoReg != 594 & ecoReg != 608 & ecoReg != 615 & ecoReg != 595 & ecoReg != 616 & ecoReg != 586 & ecoReg != 565 & ecoReg != 574 & ecoReg != 566)

# Convert dataframe of recovery metrics to an sf object
# Extract lat and lon
temp <- str_split(df_rec_tot$coords, ", ")
lon <- as.numeric(unlist(lapply(temp, function(x) x[1]))) 
lat <- as.numeric(unlist(lapply(temp, function(x) x[2]))) 
rm(temp)
df_rec_tot$lat <- lat
df_rec_tot$long <- lon
df_rec_tot_sp <- df_rec_tot

# assign coordinates
coordinates(df_rec_tot_sp) <- ~long+lat
proj4string(df_rec_tot_sp) <- CRS("+init=epsg:4326")

df_rec_tot_sp <- st_as_sf(df_rec_tot_sp,coords = c('long','lat'))

# Vector of Brazilian borders
br <- st_as_sf(getData('GADM', country = 'BRA', level = 1))# retrieve data
br <- st_transform(br, crs =  crs(df_rec_tot_sp))# convert CRS 

# convert vector to raster
rstNA <- raster(xmn= -80, xmx = -25, ymn = -40, ymx = 10, resolution = 1, crs = CRS("+init=epsg:4326")) # empty raster 
rstR80p <- rasterize(df_rec_tot_sp, rstNA, 'R80P', mean) # fill raster with recovery metric

# convert the raster to a ggplot compatible dataframe
gplot_data <- function(x, maxpixels = 50000)  {
  x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)
  coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))
  ## Extract values
  dat <- utils::stack(as.data.frame(raster::getValues(x)))
  names(dat) <- c('value', 'variable')

  dat <- dplyr::as.tbl(data.frame(coords, dat))

  if (!is.null(levels(x))) {
    dat <- dplyr::left_join(dat, levels(x)[[1]],
                            by = c("value" = "ID"))
  }
  dat
}
rstR80p_rdf <- gplot_data(rstR80p)

# filter NA
rstR80p_rdfilt <- dplyr::filter(rstR80p_rdf, value>-8000)
# get quantile range for colorbar 
perc_r80p <- quantile(rstR80p_rdfilt$value, c(.025, .975)) 

#generate plot
png(filename = file.path(ifolder, 'Spat_R80p_rst.png'), width = 500, height = 500)
ggplot()  +
  geom_tile(data = rstR80p_rdfilt, aes(x, y, fill = value), alpha = 1) +
  scale_fill_gradientn("R80p",
                      colors = brewer.pal(8, "Spectral"),
                      na.value = NA,
                      limits = c(perc_r80p[1], perc_r80p[2]), 
                      oob = scales::squish)+
    geom_sf(data = br, color = 'black', fill=NA) + 
  theme_bw()+
  ggspatial::annotation_north_arrow(location = "br", which_north = "true"  )+
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text=element_text(size=18),
        legend.text=element_text(size=18),
        legend.title=element_text(size=18))
dev.off()

# plot individual points
# perc_r80p <- quantile(df_rec_tot[!is.na(df_rec_tot$R80P),]$R80P, c(.10, .90)) 
# 
# png(filename = file.path(ifolder, 'Spat_R80p.png'), width = 500, height = 500)
# ggplot() +
#     geom_sf(data = br, color = 'black') +
#     geom_sf(data = df_rec_tot_sp[!is.na(df_rec_tot_sp$R80P),],  aes(col = R80P))+
#     scale_colour_gradientn(colours = terrain.colors(10), 
#                         limits = c(perc_r80p[1], perc_r80p[2]), oob = scales::squish)+ 
#   theme_bw()+
#   ggspatial::annotation_north_arrow(location = "br", which_north = "true"  )
# 
# dev.off()

# length(duplicated(as.numeric(df_rec_tot$id)))
# 
# getPlotData <- function(id){
#   fdts <- fireDts_tot[[id]][1]# only take the first fire
#   fdts[!is.na(fdts)] <- rollback(fdts[!is.na(fdts)], roll_to_first = T)
#   # tdist <- which(dts %in% fdts)# convert fire date to observation number
#   # distyr <- as.numeric(df_out$lossyr[id])
#   obsi <- df_out_tot[id,-c(1:4)]
#   trend <- seg_tot[[id]]$trend
#   
#   brkpts <- seg_tot[[id]]$breakpoints
#   dbr <- trend[brkpts+1]-trend[brkpts]
#   brkpts <- brkpts[which(dbr<0)]
#   if(length(brkpts)>0){tbp <- brkpts[1]}else{tbp <- seg_tot[[id]]$breakpoints[1]}
#   brkdts <- dts[tbp]
#   predt_start <- seq(brkdts, length = 2, by = "-24 months")[2]
#   distdt_end <- seq(brkdts, length = 2, by = "+12 months")[2]
#   postdt_start <- seq(brkdts, length = 2, by = "+48 months")[2]
#   postdt_end <- seq(brkdts, length = 2, by = "+72 months")[2]
# 
#   met_rri <- df_rec_tot$RRI[id]
#   met_r80p <- df_rec_tot$R80P[id]
#   met_yryr <- df_rec_tot$YrYr[id]
#   out <- list(fdts, fdts,fdts,predt_start,distdt_end,postdt_start,postdt_end,
#               obsi,trend,met_rri,met_r80p,met_yryr)
#   names(out)<- c('distyr','lossdt_start','lossdt_end','predt_start','distdt_end','postdt_start',
#                  'postdt_end','obsi','trend','met_rri','met_r80p','met_yryr')
#   return(out)
# }
# # plot time series with recovery metrics that are not equal to NA
# for (id in which(!is.na(df_rec$RRI))){
#   dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr)
# }
# 
# # fraction of time series with recovery metric that does not equal NA
# fracRRI <- sum(!is.na(df_rec_tot$RRI))/length(df_rec_tot$RRI)
# fracR80p <- sum(!is.na(df_rec_tot$R80P))/length(df_rec_tot$R80P)
# fracYrYr <- sum(!is.na(df_rec_tot$YrYr))/length(df_rec_tot$YrYr)
# 
# # plot time series with max, median and min RRI
# #id 19737
# library(lubridate)
# ylim <- c(0,1)
# sortval <- sort(df_rec_tot$RRI)
# id <- which(df_rec_tot$RRI == sortval[length(sortval)-1])
# # calcRecMetrics(as.numeric(df_out_tot[id,-c(1:4)]), fdts, 12, 2, 1, 2, 4, 2, 4, dts[seg_tot[[id]]$breakpoints[1]], selBreak = 'first', chkBrk = T, timeThres = 1.5)
# dat <- getPlotData(id[1])
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr, ylim=ylim)
# 
# id <- which(df_rec_tot$RRI == sortval[floor(length(sortval)/2)-1])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr, ylim=ylim)
# 
# id <- which(df_rec_tot$RRI == sortval[1])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr, ylim=ylim)
# 
# # plot time series with max, median and min R80p
#   png(filename = file.path(ifolder,'high_R80p.png'), width = 650, height = 300)
# ylim <- c(-0.3,1)
# sortval <- sort(df_rec_tot$R80P)
# id <- which(df_rec_tot$R80P == sortval[length(sortval)-4])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#   ylim = ylim)
#   dev.off()
# 
# png(filename = file.path(ifolder,'med_R80p.png'), width = 650, height = 300)
# R80P_vals <- df_rec_tot$R80P[!is.na(df_rec_tot$R80P)]
# R80P_med <- R80P_vals[rank(R80P_vals) == floor(length(R80P_vals)/2 +1)]
# id <- which(df_rec_tot$R80P == R80P_med)
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
#   dev.off()
# 
# png(filename = file.path(ifolder,'low_R80p.png'), width = 650, height = 300)
# id <- which(df_rec_tot$R80P == min(df_rec_tot$R80P, na.rm = T))
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
#   dev.off()
# 
# # plot time series with max, median and min YrYr
# sortval <- sort(df_rec_tot$YrYr)
# id <- which(df_rec_tot$YrYr == sortval[length(sortval)-1])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
# 
# id <- which(df_rec_tot$YrYr == sortval[floor(length(sortval)/2)])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
# 
# id <- which(df_rec_tot$YrYr == sortval[1])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
# 
# # wrong
#  png(filename = file.path(ifolder,'wrong_sec_dist.png'), width = 650, height = 300)
# 
#   id <- which(df_rec_tot$YrYr == sortval[1])
# dat <- getPlotData(id)
#   plot_aux(dts, dat$obsi,id, dat$lossdt_start,dat$lossdt_end,dat$predt_start,dat$distdt_end,
#            dat$postdt_start,dat$postdt_end,dat$trend,dat$met_rri,dat$met_r80p,dat$met_yryr,
#            ylim = ylim)
#   dev.off()