In envima/GEDItools:

Load libraries, set paths and read data

library(sf)
# library(plyr)
# library(ranger)
# library(caret)
# library(GEDItools)
library(raster)
# library(mapview)
library(ggplot2)
library(reshape2)

inpath_sen2 <- "C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel"
inpath_sen1 <- "C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel1/snap"
GEDI_sf <- readRDS("C:/Users/Alice/Uni/Projekte/GEDI/data/rGEDI_2B_clipped/bv_df_sf.rds")
LRT_rst <- raster("C:/Users/Alice/Uni/Projekte/GEDI/data/LRT/LRT/conifer_deciduous_subclasses_beechF1.tif")
lut_LRT <- read.csv("C:/Users/Alice/Uni/Projekte/GEDI/data/LRT/LUT_LRT.csv")

extract Lebensraumtypen (LRT) (Jannis Gottwald)

LRT_df <- extract(LRT_rst, GEDI_sf, df = T)
names(LRT_df) <- c("ID", "LRT_class")
LRT_name <- unlist(lapply(LRT_df$LRT_class, function(x) lut_LRT$LRT_name[match(x, lut_LRT$class)]))
LRT_type <- unlist(lapply(LRT_df$LRT_class, function(x) lut_LRT$LRT_type[match(x, lut_LRT$class)]))
LRT_df$LRT_name <- LRT_name
LRT_df$LRT_type <- LRT_type

merge LRT and GEDI

GEDI_LRT_all <- sf:::cbind.sf(GEDI_sf, LRT_df)

mask GEDI to forests from Lebensraumtypen (LRT) Wald (Jannis Gottwald)

GEDI_LRT <- GEDI_LRT_all[complete.cases(GEDI_LRT_all$LRT_class), ]

extract Sentinel 2

sen2_dirs <- list.dirs(inpath_sen2, recursive = F)
sen2_ext_lst <- lapply(sen2_dirs, function(i){
  # i <- sen2_dirs[[1]]
  basenm <- basename(i)
  # granule_dir <- list.dirs(paste0(basenm, "/", basenm, ".SAFE/GRANULE/"), recursive = F, full.names = F)
  sen2_tif <- list.files(path = paste0(i, "/", "/clipped/"), pattern = "\\.tif$", full.names =T) #.tif at end of string...no tic.aux.xml
  # date_splt <- strsplit(sen2_tif, "_")
  # date_nm <- substr(date_splt[[1]][3], 1, 8)
  sen2_stck <- stack(sen2_tif)
  names(sen2_stck) <- paste0("S2_", substr(names(sen2_stck), 11, 18), "_", 
                             c("1", "2", "3", "4", "5", "6", "7", "8", "8a", "9", "11", "12", "NDVI"))
  sen2_df <- extract(sen2_stck, GEDI_LRT, df = T)
  return(sen2_df)
})

sen2_df <- do.call(cbind, sen2_ext_lst)

extract Sentinel 1

sen1_files <- list.files(inpath_sen1, pattern =  "\\.tif$", recursive = F, full.names = T)
sen1_ext_lst <- lapply(sen1_files, function(i){
  # i <- sen1_files[[1]]
    sen1_rst <- raster(i)
    names(sen1_rst) <- paste0("S1_", substr(names(sen1_rst), 9, 16))
    sen1_df <- extract(sen1_rst, GEDI_LRT, df = T)

return(sen1_df)
})

sen1_df <- do.call(cbind, sen1_ext_lst)

merge Sentinel 2, Sentinel 1, GEDI + Lebensraumtypen data frames

GEDI_LRT_sen12 <- sf:::cbind.sf(GEDI_LRT, sen1_df, sen2_df)
GEDI_LRT_sen12 <- GEDI_LRT_sen12[,-grep("ID.", colnames(GEDI_LRT_sen12))]

create sf object with matching sen1 and 2 values per GEDI by closest date

melt
subset
cast

df_melt <- melt(as.data.frame(GEDI_LRT_sen12), id = c("beam", "shot_number", "delta_time", "latitude_bin0", 
                                       "longitude_bin0", "elev_highestreturn", "elev_lowestmode", 
                                       "pai", "fhd_normal", "omega", "pgap_theta", "cover", 
                                       "date", "filename", "LRT_class", "LRT_name", "LRT_type", "geometry", "ID"))
df_melt$date_variable <- as.Date(substr(df_melt$variable, 4, 11), format = "%Y%m%d")
time_diff_all <- abs(df_melt$date_variable - df_melt$date) < 10
df_melt_flt <- df_melt[time_diff_all,]
df_melt_flt$sensor <- substr(df_melt_flt$variable, 1,2)
splt <- strsplit(as.character(df_melt_flt$variable), "_")
splt_df <- do.call(rbind, splt)
df_melt_flt$band <- paste0(splt_df[,1], "_", splt_df[,3])
#cast
cast_df <- dcast(df_melt_flt, ID~band, sum) 

sf_cast_mrg <- merge(GEDI_LRT_sen12, cast_df, by = "ID")
model_sf <- sf_cast_mrg[!grepl("2019", colnames(sf_cast_mrg))]
saveRDS(object = model_sf, file = "C:/Users/Alice/Uni/Projekte/GEDI/data/modeling/model_sf.rds")

create data frame for modelling filter data frame, so each GEDI point refers to one Band value, depending on recording date

to do

subset for deciduos and coniferous

sf_con <- GEDI_LRT_sen12[GEDI_LRT_sen12$LRT_type == "coniferous", ]
sf_dec <- GEDI_LRT_sen12[GEDI_LRT_sen12$LRT_type == "deciduous", ]

boxplot pai[dec] ~ date

ggplot(sf_dec, aes(x = date, y = pai, group = date))+
  geom_boxplot()

boxplot sen2[dec][NDVI] ~ date

df_dec <- as.data.frame(sf_dec)
dec_NDVI <- df_dec[,c("LRT_type", "S2_20190420_NDVI", "S2_20190624_NDVI", "S2_20190629_NDVI", "S2_20190823_NDVI", "date")]
dec_NDVI_melt <- melt(dec_NDVI, id = c("date", "LRT_type"), value.name = "NDVI")

dec_NDVI_melt$date_NDVI <- as.Date(substr(dec_NDVI_melt$variable, 4, 11), format = "%Y%m%d")

###hard code
time_diff_NDVI <- abs(dec_NDVI_melt$date_NDVI - dec_NDVI_melt$date) < 10
dec_NDVI_flt <- dec_NDVI_melt[time_diff_NDVI,]

ggplot(dec_NDVI_flt, aes(x = date_NDVI, y = NDVI, group = date_NDVI))+
  geom_boxplot()

boxplot sen1[dec] ~ date

df_dec <- as.data.frame(sf_dec)
dec_s1 <- df_dec[,c("LRT_type", "S1_20190422", "S1_20190522", "S1_20190625", "S1_20190628", "S1_20190824", "date")]
dec_s1_melt <- melt(dec_s1, id = c("date", "LRT_type"), value.name = "gamma")

dec_s1_melt$date_s1 <- as.Date(substr(dec_s1_melt$variable, 4, 11), format = "%Y%m%d")

###hard code
time_diff_s1 <- abs(dec_s1_melt$date_s1 - dec_s1_melt$date) < 10
dec_s1_flt <- dec_s1_melt[time_diff_s1,]

ggplot(dec_s1_flt, aes(x = date_s1, y = gamma, group = date_s1))+
  geom_boxplot()

boxplot pai[con] ~ date

ggplot(sf_con, aes(x = date, y = pai, group = date))+
  geom_boxplot()

boxplot sen2[con][NDVI] ~ date

df_con <- as.data.frame(sf_con)
con_NDVI <- df_con[,c("LRT_type", "S2_20190420_NDVI", "S2_20190624_NDVI", "S2_20190629_NDVI", "S2_20190823_NDVI", "date")]
con_NDVI_melt <- melt(con_NDVI, id = c("date", "LRT_type"), value.name = "NDVI")

con_NDVI_melt$date_NDVI <- as.Date(substr(con_NDVI_melt$variable, 4, 11), format = "%Y%m%d")

###hard code
time_diff_NDVI <- abs(con_NDVI_melt$date_NDVI - con_NDVI_melt$date) < 10
con_NDVI_flt <- con_NDVI_melt[time_diff_NDVI,]

ggplot(con_NDVI_flt, aes(x = date_NDVI, y = NDVI, group = date_NDVI))+
  geom_boxplot()

boxplot sen1[con] ~ date

df_con <- as.data.frame(sf_con)
con_s1 <- df_con[,c("LRT_type", "S1_20190422", "S1_20190522", "S1_20190625", "S1_20190628", "S1_20190824", "date")]
con_s1_melt <- melt(con_s1, id = c("date", "LRT_type"), value.name = "NDVI")

con_s1_melt$date_s1 <- as.Date(substr(con_s1_melt$variable, 4, 11), format = "%Y%m%d")

###hard code
time_diff_s1 <- abs(con_s1_melt$date_s1 - con_s1_melt$date) < 10
con_s1_flt <- con_s1_melt[time_diff_s1,]

ggplot(con_s1_flt, aes(x = date_s1, y = NDVI, group = date_s1))+
  geom_boxplot()