above_atl08_airborne.Rmd

---
title: "above_atl08_airborne"
output:
  html_document:
    df_print: paged
---

```{r}
library(tidyverse)
library(plyr)
library(fs)
library(reshape2)
library(stringr)

give.n <- function(x){
  # https://stackoverflow.com/questions/28846348/add-number-of-observations-per-group-in-ggplot2-boxplot
  #return(c(y = min(x), label = length(x)))
  return(c(y = 0, label = length(x))) 
  # experiment with the multiplier to find the perfect position
}


```
```{r fig.width=8, fig.height=6}
#
# TERRAIN
#
# Finland results from Amy - 
# https://docs.google.com/spreadsheets/d/1kF1E0A9TJeAq9IEw94njXhtDz9ZbMtF6Dn3zmQMm8ZI/edit?usp=sharing
file='ATL08 validation Finland - terrain.csv'
results_finland_terrain = data.frame(read.csv(path("C:\\Users\\pmontesa\\Downloads", file)))
names(results_finland_terrain) = c("Beam","Day.Night" ,"Snow.Cover", "Bias","MAE","RMSE","N.samples")
results_finland_terrain.m = melt(results_finland_terrain , measure.vars = c("Bias", "RMSE", "MAE"))

ggplot(results_finland_terrain.m , aes(x=Day.Night, y=value, fill=Beam, size=N.samples)) + 
  facet_grid(variable~Snow.Cover, scales = "free_y") + 
  geom_point(shape=21, alpha=0.75) + 
  theme_bw() +
  labs(fill="ATL08 beam type", size="# obs.", y="Error (m)", x=NULL,  title = "Error estimates from validation of ATL08 median terrain elevation", subtitle = "ATL08 metric: h_te_median") +
  theme(text = element_text(size=14))+
  guides(fill = guide_legend(order = 1, override.aes = list(size=6)),
         size = guide_legend(order = 2))
  #scale_fill_manual(values = c("dodgerblue", "orange", "green"))
  #geom_text(aes(y = value, label=N.samples), parse = TRUE, hjust=-0.5, vjust=0, size=2, color = "black", fill="white", label.padding=unit(0.1,"lines"))

```
```{r fig.width=10, fig.height=5}
#
# CANOPY
#
# Finland results from Amy
# https://docs.google.com/spreadsheets/d/1kF1E0A9TJeAq9IEw94njXhtDz9ZbMtF6Dn3zmQMm8ZI/edit?usp=sharing
file='ATL08 validation Finland - canopy.csv'
results_finland = data.frame(read.csv(path("C:\\Users\\pmontesa\\Downloads", file)))
results_finland$Snow.Cover = mapvalues(factor(results_finland$Snow.Cover), from = levels(results_finland$Snow.Cover), to=c("No snow & Snow", "No snow", "Snow"))
results_finland$Day.Night = mapvalues(factor(results_finland$Day.Night), from = levels(results_finland$Day.Night), to=c("Day & Night", "Day", "Night"))
names(results_finland) = c("Beam","Day.Night" ,"Snow.Cover", "MODIS.Landcover", "Season","Canopy.N","Bias","RMSE","Bias_prc","RMSE_prc")
results_finland.m = reshape2::melt(results_finland %>% select(-c( "Bias_prc", "RMSE_prc")), measure.vars = c("Bias", "RMSE"))
results_finland_prc.m = reshape2::melt(results_finland %>% select(-c( "Bias", "RMSE")), measure.vars = c("Bias_prc", "RMSE_prc"), value.name="value_prc")

results_finland.m = cbind(results_finland.m, results_finland_prc.m %>% select(value_prc))

options(scipen=10000)

ggplot(results_finland.m %>% filter(Beam=="Strong"), aes(x=Day.Night, y=value, fill=Season, size=Canopy.N, group=Season)) + 
  facet_grid(variable~Snow.Cover, scales = "free_y") + 
  geom_point( alpha=0.75, position=position_dodge(width=0), shape=21) + 
  #geom_text(data=results_finland.m %>% filter(Snow.Cover=="No snow"), aes(label=MODIS.Landcover), color='black', size=3, hjust=-0.2, angle=45)+
    geom_text(data=results_finland.m %>% filter(Beam=="Strong"), aes(label=paste0(round(value_prc, dig=1),"%")), color='black', size=3,angle=0,  position = position_nudge(x = 0.4))+
  theme_bw() +
  labs(color="Season", size="# obs.",  y="Error (m)", x=NULL,  title = "Error estimates from validation of ATL08 canopy height", subtitle = "ATL08 metric: 'h_can'") +
  theme(text = element_text(size=14))+
  scale_fill_manual(values = c("dodgerblue", "orange"))+
  #scale_fill_brewer(palette = "Set1")+
  scale_size(range = c(4,10), breaks = c(50000,150000,350000))+
  #scale_shape_manual(values = c(21,24), name="Land cover")+
  guides(fill = guide_legend(order = 1, override.aes = list(size=6)),
         size = guide_legend(order = 2 )#,
         #shape = guide_legend(order = 3, override.aes = list(size=5))
         )



```

```{r}
#
# Functions to use with the ATL08__*__Zonal CSV files from M. Wooten
#
lm_eqn = function(df){
      y=df[,1]
      x=df[,2]
      m = lm(y ~ x, df);
      
      eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
           list(a = format(coef(m)[1], digits = 2), 
                b = format(coef(m)[2], digits = 2), 
               r2 = format(summary(m)$r.squared, digits = 3)))
      int = format(coef(m)[1], digits = 2)
      slope = format(coef(m)[2], digits = 2)
      # https://stackoverflow.com/questions/43123462/how-to-obtain-rmse-out-of-lm-result
      
      rss   <- c(crossprod(m$residuals))
      mse   <- rss / length(m$residuals)
      
      rmse  <- format(sqrt(mse),dig=2)
      sig2  <- format(rss / m$df.residual, dig=2) # Pearson estimated residual variance (as returned by summary.lm):
      n     <- nobs(m)
      r2    <- format(summary(m)$r.squared, dig = 2)
      mae   <- format(mean(abs(m$residuals)), dig = 2)
      mape  <- format(mean(abs((m$residuals)/y))*100, dig=2)  #https://stats.stackexchange.com/questions/299712/what-are-the-shortcomings-of-the-mean-absolute-percentage-error-mape

      #return(as.character(as.expression(eq)))
      #return(data.frame(cbind(int, slope, r2, rmse, sig2, n, mae, mape)))
      df = data.frame(cbind(int, slope, r2, rmse, sig2, n, mae, mape)) #%>% dplyr::mutate_all(as.numeric)
      return(df)
}
lm_mod = function(df){
      y=df[,1]
      x=df[,2]
      m = lm(y ~ x, df);
      return(m)
}

do_height_plots = function(df, ref_name, filt_type, filt_str, y_var, x_var, outdir){
  
  title_info = paste0(" [ ", filt_type, " : ", filt_str," ] ")
  outplotdir = paste0("output_",format(Sys.time(), "%Y%m%d"))
  dir.create(path(outdir, outplotdir), showWarnings = FALSE)
  sz_txt = 3

  pal_rev_spectral9 = c("#3288BD","#66C2A5","#ABDDA4", "#E6F598","#FFFFBF", "#FEE08B", "#FDAE61", "#F46D43", "#D53E4F")
  ymax = 25
  leg_breaks = c(1,50,100,500)
  leg_lims = c(min(leg_breaks), max(leg_breaks))
  
  val_plot_geom_list = list(

          guides(colour = guide_legend(override.aes = list(alpha = 1))),
          geom_vline(xintercept = 0, linetype = "dashed"),
          geom_hline(yintercept = 0, linetype = "dashed"),
          geom_abline(intercept = 0, slope = 1, "dotted"),
          stat_smooth(method = "lm", col = "black", linetype = "dashed"),
          scale_x_continuous(breaks = seq(0,ymax,5), limits = c(0,ymax+1)),
          scale_y_continuous(breaks = seq(0,ymax,5), limits = c(0,ymax+1)),
          scale_fill_gradientn( labels = leg_breaks, breaks =leg_breaks, limits = leg_lims, colours=pal_rev_spectral9, trans = "log", name='# obs.'),
          theme_minimal()+theme(axis.text = element_text(size = 6)),
          coord_fixed()
          )


  # # Multiple Scattering warning flag. The multiple scattering warning flag (ATL09 parameter msw_flag) has values from -1 to 5 where zero means no multiple scattering and 5 the greatest. If no layers were detected, then msw_flag = 0. If blowing snow is detected and its estimated optical depth is greater than or equal to 0.5, then msw_flag = 5. If the blowing snow optical depth is less than 0.5, then msw_flag = 4. If no blowing snow is detected but there are cloud or aerosol layers detected, the msw_flag assumes values of 1 to 3 based on the height of the bottom of the lowest layer: < 1 km, msw_flag = 3; 1-3 km, msw_flag = 2; > 3km, msw_flag = 1. A value of -1 indicates that the signal to noise of the data was too low to reliably ascertain the presence of cloud or blowing snow. We expect values of -1 to occur only during daylight.
  # df$msw_flg   = mapvalues(factor(df$msw_flg), from = seq(-1,5,1), to=c("Low SNR", "No MS", "No snow; aer. > 3km", "No snow; aer. 1-3km", "No snow; aer. < 1km", "Snow; AOD<0.5", "Snow; AOD>0.5"))
  
  # Create the equation labels for the two groups
  #source(path("C:\\Users\\pmontesa\\Google Drive\\Work\\R","plot_functions.R"))
  
 
  # Facet by night_flg and beam_type
  df_mod_1  <- dlply(df %>% dplyr::select_(y_var, x_var, "msw_flg","cloud_flg", "night_flg", "tcc_flg", "gt", "beam_type" ), .( night_flg, beam_type), c("lm_mod"))
  df_eq_1   <- ddply(df %>% dplyr::select_(y_var, x_var, "msw_flg","cloud_flg", "night_flg", "tcc_flg", "gt", "beam_type" ), .( night_flg, beam_type), c("lm_eqn") )
  #df_eq_1$slope = round(as.numeric(as.character(df_eq_1$slope)), dig=2)
  
  plot_ht_val_beamtype = ggplot(df, aes_string(y=y_var, x=x_var)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("Reference height: ",x_var)) + 
    geom_bin2d(binwidth = c(0.75, 0.75))+
    #geom_point(alpha=0.25)+
    labs(y=paste0("Height (m); ATL08, ",y_var), x=paste0("Height (m); ", x_var))+
    val_plot_geom_list+
    geom_text(data=df_eq_1, aes(x = Inf, y = 0,    label=paste0("r^2 ==", r2)),      inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_1, aes(x = Inf, y = 0+3,  label=paste0("mae ==", mae)),     inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_1, aes(x = Inf, y = 0+3*2,label=paste0("rmse ==", rmse)),   inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_1, aes(x = Inf, y = 0+3*3,label=paste0("bias ==", int)),    inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_1, aes(x = Inf, y = 0+3*4,label=paste0("slope ==", slope)), inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_1, aes(x = Inf, y = Inf,label=paste0("n ==", n)),         inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    #facet_wrap(.~night_flg~)+
    facet_grid(night_flg~beam_type)+
    theme(axis.text = element_text(size = 10), strip.text = element_text(size = 12))
  ggsave(plot = plot_ht_val_beamtype,
         file = path(outdir,outplotdir, paste0(paste("plot_ht_val_beamtype",ref_name,filt_type,y_var,x_var,sep="_"),".png")),
         device = 'png', dpi = 300, width = 12, height = 6)
  
  
  # Facet by ref TCC bin and night_flg
  df_mod_2 <- dlply(df %>% dplyr::select_(y_var, x_var, "cloud_flg", "night_flg", "tcc_flg", "gt","ref_tcc_bin", "beam_type" ),.(ref_tcc_bin, night_flg, beam_type),c("lm_mod") )
  df_eq_2  <- ddply(df %>% dplyr::select_(y_var, x_var, "cloud_flg", "night_flg", "tcc_flg", "gt","ref_tcc_bin", "beam_type" ),.(ref_tcc_bin, night_flg, beam_type),c("lm_eqn") )
  #df_eq_2$slope = round(as.numeric(as.character(df_eq_2$slope)), dig=2)
  
  plot_ht_val_TCC = ggplot(df, aes_string(y=y_var, x=x_var)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("Reference height: ", x_var)) + 
    geom_bin2d(binwidth = c(0.75, 0.75))+
    #geom_point(alpha=0.25)+
    labs(y=paste0("Height (m); ATL08, ",y_var), x=paste0("Height (m); ", x_var))+
    val_plot_geom_list+
    geom_text(data=df_eq_2, aes(x = Inf, y = 0,     label=paste0("r^2 ==", r2)),     inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_2, aes(x = Inf, y = 0+3,   label=paste0("mae ==", mae)),    inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+    
    geom_text(data=df_eq_2, aes(x = Inf, y = 0+3*2, label=paste0("rmse ==", rmse)),  inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_2, aes(x = Inf, y = 0+3*3, label=paste0("bias ==", int)),   inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_2, aes(x = Inf, y = 0+3*4, label=paste0("slope ==", slope)),inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_2, aes(x = Inf, y = Inf, label=paste0("n ==", n)),          inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    facet_grid(night_flg+beam_type~ref_tcc_bin)+
    theme(axis.text = element_text(size = 10), strip.text = element_text(size = 12))
    ggsave(plot = plot_ht_val_TCC,
         file = path(outdir,outplotdir, paste0(paste("plot_ht_val_TCC",ref_name,filt_type,y_var,x_var,sep="_"),".png")),
         device = 'png', dpi = 300, width = 12, height = 6)
    
  # Facet by seg_snow and night_flg and group by ref_tcc_bin
  df_mod_3 <- dlply(df %>% dplyr::select_(y_var, x_var, "cloud_flg", "night_flg", "tcc_flg", "gt","ref_tcc_bin","seg_snow", "beam_type" ),.(seg_snow, night_flg, ref_tcc_bin, beam_type),c("lm_mod") )
  df_eq_3  <- ddply(df %>% dplyr::select_(y_var, x_var, "cloud_flg", "night_flg", "tcc_flg", "gt","ref_tcc_bin","seg_snow", "beam_type" ),.(seg_snow, night_flg, ref_tcc_bin, beam_type),c("lm_eqn") )
  #df_eq_3$slope = round(as.numeric(as.character(df_eq_3$slope)), dig=2)
  
  plot_ht_val_snow_daynight_TCC = ggplot(df, aes_string(y=y_var, x=x_var)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("Reference height: ", x_var)) + 
    geom_bin2d(binwidth = c(0.75, 0.75))+
    #geom_point(alpha=0.25)+
    labs(y=paste0("Height (m); ATL08, ",y_var), x=paste0("Height (m); ", x_var))+
    val_plot_geom_list+
    geom_text(data=df_eq_3, aes(x = Inf, y = 0,     label=paste0("r^2 ==", r2)),     inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_3, aes(x = Inf, y = 0+3,   label=paste0("mae ==", mae)),    inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+  
    geom_text(data=df_eq_3, aes(x = Inf, y = 0+3*2, label=paste0("rmse ==", rmse)),  inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_3, aes(x = Inf, y = 0+3*3, label=paste0("bias ==", int)),   inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_3, aes(x = Inf, y = 0+3*4, label=paste0("slope ==", slope)),inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    geom_text(data=df_eq_3, aes(x = Inf, y = Inf, label=paste0("n ==", n)),          inherit.aes=FALSE, parse = TRUE, hjust="inward", vjust="inward", size=sz_txt)+
    facet_grid(night_flg+beam_type~seg_snow+ref_tcc_bin)+
    theme(axis.text = element_text(size = 10), strip.text = element_text(size = 12))
    ggsave(plot = plot_ht_val_snow_daynight_TCC,
         file = path(outdir,outplotdir, paste0(paste("plot_ht_val_snow_daynight_TCC",ref_name,filt_type,y_var,x_var,sep="_"),".png")),
         device = 'png', dpi = 300, width = 12, height = 6)
      
  return(list(plot_ht_val_beamtype, plot_ht_val_TCC, plot_ht_val_snow_daynight_TCC, df_eq_1, df_eq_2, df_eq_3, df_mod_1, df_mod_2, df_mod_3))
}

# 
pal_ht    <- colorRampPalette(c('#636363','#fc8d59','#fee08b','#ffffbf','#d9ef8b','#91cf60','#1a9850'))(20)

```


```{r}
#
# Read in ATL08 csv: !! huge file !!
#

atl08_dir = "D:/projects/above_icesat2/atl08"

# Read un main ATL08 filtered file for North America
if(!exists("df_atl08")){
  
  library(data.table)
  df_atl08 = fread(path(atl08_dir, "ATL08_merged_003_01_filt_45_90_-180_-13.csv")) 
  
}
dim(df_atl08)

# Get beam_type
df_atl08 = df_atl08 %>%
    filter(orb_orient != 2) %>%
    mutate(beam_type = case_when(
                  orb_orient == 1 & grepl("r$", df_atl08$gt) == TRUE  ~ "Strong",
                  orb_orient == 1 & grepl("r$", df_atl08$gt) == FALSE ~ "Weak",
                  orb_orient == 0 & grepl("r$", df_atl08$gt) == TRUE  ~ "Weak",
                  orb_orient == 0 & grepl("r$", df_atl08$gt) == FALSE ~ "Strong"
                  )
    )

```

```{r fig.height=4, fig.width=18}
#
# Note: These zonal tables should be updated with terrain elev from lidar stacks.
#
# Read in Zonal Stats files
dir_zonal = "D:\\projects\\3dsi\\ZonalStats\\"
library(data.table)
df_atl08_gliht_in = fread(path(dir_zonal, "ATL08__GLIHT__ZonalStats.csv"))
df_atl08_lvis_in  = fread(path(dir_zonal, "ATL08__LVIS__ZonalStats.csv"))

#/att/gpfsfs/briskfs01/ppl/mwooten3/3DSI/ZonalStats/README/GLiHT_layers.txt

# L1: <stackName>_CHM_warp.tif
# L2: <stackName>_tree_p10_warp.tif
# L3: <stackName>_tree_p20_warp.tif
# L4: <stackName>_tree_p30_warp.tif
# L5: <stackName>_tree_p40_warp.tif
# L6: <stackName>_tree_p50_warp.tif
# L7: <stackName>_tree_p60_warp.tif
# L8: <stackName>_tree_p70_warp.tif
# L9: <stackName>_tree_p80_warp.tif
# L10: <stackName>_tree_p90_warp.tif
# L11: <stackName>_tree_p100_warp.tif
# L12: <stackName>_tree_fcover_warp.tif
# L13: <stackName>_slope_warp.tif
# L14: <stackName>_aspect_warp.tif
# L15: NA_standage_v2_warp.tif
# L16: NA_focal_stats_std_v2_warp.tif
# L17: NA_focal_stats_mask_v2_warp.tif
# L18: <stackName>_LS7_CC_warp.tif
# L19: MCD12Q1_A2017001_LC_Type1_warp.tif
# L20: boreal_clust_25_100_2019_10_8_warp.tif
# L21: PCA_NaN_1_093019_warp.tif
# L22: PCA_NaN_2_093019_warp.tif
# L23: PCA_NaN_3_093019_warp.tif
# L24: wc2.0_bio_30s_01_warp.tif
# L25: wc2.0_bio_30s_10_warp.tif
# L26: wc2.0_bio_30s_18_warp.tif
# L27: UiO_PEX_PERPROB_3.0_20171201_2000_2016_warp.tif
# L28: BDTICM_M_250m_warp.tif

#/att/gpfsfs/briskfs01/ppl/mwooten3/3DSI/ZonalStats/README/LVIS_layers.txt

# L1: <stackName>_RH010_mean_15m_warp.tif
# L2: <stackName>_RH020_mean_15m_warp.tif
# L3: <stackName>_RH030_mean_15m_warp.tif
# L4: <stackName>_RH040_mean_15m_warp.tif
# L5: <stackName>_RH050_mean_15m_warp.tif
# L6: <stackName>_RH060_mean_15m_warp.tif
# L7: <stackName>_RH070_mean_15m_warp.tif
# L8: <stackName>_RH080_mean_15m_warp.tif
# L9: <stackName>_RH090_mean_15m_warp.tif
# L10: <stackName>_RH098_mean_15m_warp.tif
# L11: <stackName>_retge_1p37_15m_warp.tif
# L12: <stackName>_ZG_mean_15m_slope_warp.tif
# L13: <stackName>_ZG_mean_15m_aspect_warp.tif
# L14: NA_standage_v2_warp.tif
# L15: NA_focal_stats_std_v2_warp.tif
# L16: NA_focal_stats_mask_v2_warp.tif
# L17: NA_standage_v2_mmu_warp.tif
# L18: NA_disturbance_v2_mmu_warp.tif
# L19: <stackName>_LS7_CC_warp.tif
# L20: C2C_change_year_warp.tif
# L21: C2C_change_type_warp.tif
# L22: MCD12Q1_A2017001_LC_Type1_warp.tif
# L23: boreal_clust_25_100_2019_10_8_warp.tif
# L24: PCA_NaN_1_093019_warp.tif
# L25: PCA_NaN_2_093019_warp.tif
# L26: PCA_NaN_3_093019_warp.tif
# L27: wc2.0_bio_30s_01_warp.tif
# L28: wc2.0_bio_30s_10_warp.tif
# L29: wc2.0_bio_30s_18_warp.tif
# L30: UiO_PEX_PERPROB_3.0_20171201_2000_2016_warp.tif
# L31: BDTICM_M_250m_warp.tif


#
# Data cleaning: GLiHT
#
#Handle no data
df_atl08_gliht_in[df_atl08_gliht_in == -99] <- NA
# Set col names
names(df_atl08_gliht_in)[grep(".median", names(df_atl08_gliht_in))][1:12] = c("CHM.median", paste0("rh", seq(10,100,10), ".median"), "tcc1.37.median")
df_atl08_gliht_in$ref_ht = df_atl08_gliht_in$rh100.median
df_atl08_gliht_in$ref_tcc = df_atl08_gliht_in$tcc1.37.median * 100
df_atl08_gliht_in$lidar_name = "GLiHT"

# Data cleaning: LVIS
#
#Handle no data

# Set col names
names(df_atl08_lvis_in)[grep(".median", names(df_atl08_lvis_in))][1:11] = c(paste0("rh", c(seq(10,90,10), 98), ".median"), "tcc1.37.median")
df_atl08_lvis_in$ref_ht = df_atl08_lvis_in$rh98.median
df_atl08_lvis_in$ref_tcc = as.numeric(as.character(df_atl08_lvis_in$tcc1.37.median))
df_atl08_lvis_in$lidar_name = "LVIS"

#
# Combine GLiHT and LVIS
#
df_atl08_LIDAR = rbind(df_atl08_gliht_in %>% 
                         drop_na(ref_ht) %>%
                         #dplyr::select(ref_ht, ref_tcc, h_can ,lat, lon, d, m, yr, gt, n_toc_ph, ter_slp, h_te_unc, cloud_flg, seg_snow, night_flg, msw_flg, tcc_flg, stackName, lidar_name) %>%
                         dplyr::select(contains('ref'), lat, lon, stackName, lidar_name), #%>%
                         #dplyr::mutate_at(vars(contains('flg')), funs(as.factor)),
                       df_atl08_lvis_in  %>% 
                         drop_na(ref_ht) %>%
                         dplyr::select(contains('ref'), lat, lon, stackName, lidar_name) #%>%
                         #dplyr::mutate_at(vars(contains('flg')), funs(as.factor))
                       ) #%>%
            #dplyr::mutate_at(vars(lat, lon), funs(as.numeric)) %>%
            #dplyr::mutate_at(vars(lat, lon), funs(round(., 6)))

DO_TABULAR_JOIN=T
if(DO_TABULAR_JOIN){
  
  #
  # Tabular Join: *new* fields from filtered ATL08 not available in Maggie's Zonal db
  #
  df_atl08_LIDAR = left_join(df_atl08_LIDAR , 
                             df_atl08 %>% 
                              #dplyr::mutate_at(vars(lat, lon), funs(round(., 6))) %>%
                              #dplyr::select(contains('rh'), lat, lon, asr, h_dif_ref, ter_flg, ph_rem_flg, dem_rem_flg, seg_wmask, lyr_flg, beam_type, orb_num, orb_orient), by = c("lat", "lon"))
                              dplyr::select(contains('rh'), contains('_flg'), contains('mask'), contains('snow'), contains('can'), yr, d, n_toc_ph, h_te_unc, ter_slp, lat, lon, asr, h_dif_ref, beam_type, orb_num, orb_orient, gt), by = c("lat", "lon")) %>%
    dplyr::mutate_at(vars(contains('flg')), funs(as.factor)) %>%
    dplyr::mutate_at(vars(contains('gt')), funs(as.factor)) %>%
    dplyr::mutate_at(vars(contains('beam')), funs(as.factor)) %>%
    dplyr::mutate_at(vars(contains('snow')), funs(as.factor))
}

# Make treecover bins
cc_bin=20
df_atl08_LIDAR$ref_tcc_bin = cut(df_atl08_LIDAR$ref_tcc, seq(0,100,cc_bin), labels= as.character(seq(cc_bin,100,cc_bin)) )

# Bin the reference ht values
df_atl08_LIDAR$ref_ht_bin = cut(df_atl08_LIDAR$ref_ht, seq(0,30,1), labels= as.character(seq(1,30,1)) )

# Assign names to categorical variables that can be used to subset data
df_atl08_LIDAR$cloud_flg = mapvalues(factor(df_atl08_LIDAR$cloud_flg),from = c(0,1,2,3,4,5), to=c("High conf. clear","Med conf. clear", "Low conf. clear", "Low conf. cloudy ", "Med conf. cloudy", "High conf. cloudy"))
df_atl08_LIDAR$seg_snow  = mapvalues(factor(df_atl08_LIDAR$seg_snow), from = c(0,1,2,3),     to=c("ice free water", "No snow", "Snow", "ice"))
df_atl08_LIDAR$night_flg = mapvalues(factor(df_atl08_LIDAR$night_flg),from = c(0,1),         to=c("day","night"))
df_atl08_LIDAR$tcc_flg   = mapvalues(factor(df_atl08_LIDAR$tcc_flg),  from = c(0,1),         to=c("<= 5%","> 5%"))
df_atl08_LIDAR$msw_flg   = mapvalues(factor(df_atl08_LIDAR$msw_flg), from = seq(-1,5,1), to=c("Low SNR", "No MS", "No snow; aer. > 3km", "No snow; aer. 1-3km", "No snow; aer. < 1km", "Snow; AOD<0.5", "Snow; AOD>0.5"))

n_lidar = dim(df_atl08_LIDAR %>% drop_na(ref_ht, h_can))[1]

ggplot(df_atl08_LIDAR %>% drop_na(ref_ht, h_can), aes(factor(gt), fill=gt)) + 
    geom_histogram(stat="count") +
    labs(title='Summary of ICESAt-2 ATL08 observations across airborne lidar',
         subtitle = paste0('# crossover obs.: ', n_lidar ),
         fill='ATL08\nbeam',x=NULL)+
    facet_grid(beam_type~night_flg+msw_flg)

# Save the set of observations before filtering
write.csv(df_atl08_LIDAR, file=path(dir_zonal,"df_atl08_LIDAR.csv"))

```

```{r}
#
# Base filtering
#
# Remove points with bad data
df_atl08_LIDAR_filt = df_atl08_LIDAR %>%
              #dplyr::filter(yr >= 2018 & yr < 2050) %>%
              drop_na(ref_ht) %>%
        dplyr::filter(n_toc_ph>0) %>%
        dplyr::filter(h_can<thresh_h_can) %>%
        dplyr::filter(h_te_unc!=max(h_te_unc))%>%
        dplyr::filter(ter_slp!=max(ter_slp)) %>%
        # These should remove a lot of data
        #dplyr::filter(ter_flg==0 | is.na(ter_flg)) %>% # 0 or 1; 0 is below the threshold for deviation from a ref DEM; dont remove NAs, b/c then hard to assess if joins are working as expected!
        dplyr::filter(msw_flg=='No MS') %>%
        dplyr::filter(lat > 45)

# CHeck NAs after join....
#ggplot(df_atl08_LIDAR, aes(x=factor(ter_flg), fill=lidar_name)) + geom_histogram(position = "stack", stat = "count")


n_lidar_filt = dim(df_atl08_LIDAR_filt %>% drop_na(ref_ht, h_can))[1]

ggplot(df_atl08_LIDAR_filt, aes(factor(gt), fill=gt)) + 
    geom_histogram(stat="count") +
    labs(title='Summary of filtered ICESAt-2 ATL08 observations across airborne lidar',
         subtitle = paste0("# crossover obs.: ", n_lidar_filt),
         fill='ATL08\nbeam',x=NULL)+
    facet_grid(beam_type~night_flg+msw_flg)


# Save the set of observations for validation
write.csv(df_atl08_LIDAR_filt, file=path(dir_zonal,"df_atl08_LIDAR_filt.csv"))

```

```{r fig.height=5, fig.width=16, message=FALSE, warning=FALSE}
COVER= FALSE
thresh_h_can = 30
ref_name = "LiDAR"
p_title = paste0("Crossover locations of ICESat-2 ATL08 and ", ref_name)
p_subtitle = "Reference lidar canopy height (m, GLiHT: rh100, LVIS: rh98)"

out_ext = '.png'
out_ext_type = 'png'

####
# Create a seequence of facet plots based on data with different filters

# Filter 1
filt_type = "No multiple scattering"
filt_str = ""
title_info = paste0(" [ ", filt_type, " : ", filt_str," ] ")

cent_lat = mean(df_atl08_LIDAR_filt$lat)
cent_lon = mean(df_atl08_LIDAR_filt$lon)

#ggplot(df_atl08_LIDAR_filt) +
map_lidar = ggplot(map_data("world") ) + 
    ggplot2::coord_map("ortho", orientation = c(cent_lat, cent_lon, 0), xlim = c(-165,-100), ylim = c(50,70) ) +
    geom_polygon(aes(x = long, y = lat, group = group), fill="grey40") +
    geom_point(data=df_atl08_LIDAR_filt, alpha=0.2, aes(x = lon, y = lat, col=ref_ht), size = 0.05, alpha = 0.05) +
    scale_colour_gradientn(colors=pal_ht, limits=c(0,20), name = "Height (m)") +
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = p_subtitle)+
    facet_wrap(~lidar_name)+
    theme( strip.text.x = element_text(size=20) ) + #, axis.ticks.y=element_blank(), axis.text = element_blank(), axis.title = element_blank())+
    #theme_minimal()+
    labs(x= NULL, y= NULL, caption = "above_atl08_airborne.Rmd")
print(map_lidar)



```


```{r fig.height=10, fig.width=18}
p_title = paste0("Zonal stats analysis of ATL08 and ",ref_name)
p_subtitle = "ICESat-2 ATL08 Canopy Height (m, h_can)" 

#
# FILTERING to examine snow/no-snow
#

list_out_gliht = do_height_plots(df_atl08_LIDAR_filt %>% 
                                   filter(str_detect(seg_snow, 'snow')) %>%
                        drop_na(ref_tcc_bin, ref_ht_bin) %>%
                        filter(lidar_name=='GLiHT')
                      , ref_name, paste0(filt_type,"_GLiHT"), filt_str, "h_can", "ref_ht", dir_zonal)
print(list_out_gliht[1:3])

list_out_lvis = do_height_plots(df_atl08_LIDAR_filt %>% 
                                  filter(str_detect(seg_snow, 'snow')) %>%
                        drop_na(ref_tcc_bin, ref_ht_bin) %>%
                        filter(lidar_name=='LVIS')
                      , ref_name, paste0(filt_type,"_LVIS"), filt_str, "h_can", "ref_ht", dir_zonal)
print(list_out_lvis[1:3])

list_out_ALL = do_height_plots(df_atl08_LIDAR_filt %>% 
                                  filter(str_detect(seg_snow, 'snow')) %>%
                        drop_na(ref_tcc_bin, ref_ht_bin) 
                      , ref_name, paste0(filt_type,"_All"), filt_str, "h_can", "ref_ht", dir_zonal)
print(list_out_ALL[1:3])

```

```{r fig.height=6, fig.width=6}
#
# Results of GLiHT validation
#
results_gliht = data.frame(list_out_gliht[6]) %>% 
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.character)) %>%   
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.numeric)) %>%
  droplevels()
results_gliht$ref_name = 'GLiHT'

results_lvis = data.frame(list_out_lvis[6]) %>% 
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.character)) %>%   
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.numeric)) %>%
  droplevels()
results_lvis$ref_name = 'LVIS'

results_all = data.frame(list_out_ALL[6]) %>% 
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.character)) %>%   
  mutate_at(vars("int","slope","r2","rmse","sig2","n","mae","mape","ref_tcc_bin"),funs(as.numeric)) %>%
  droplevels()
results_all$ref_name = 'All'

results = rbind(results_gliht, results_lvis, results_all)

results$seg_snow = mapvalues(factor(results$seg_snow), from = levels(results$seg_snow), to=c("No snow", "Snow"))
results$night_flg  = mapvalues(factor(results$night_flg ), from = levels(results$night_flg ), to=c("Day", "Night"))
names(results)[grep('int', names(results))] <- 'bias'
results.m = reshape2::melt(results, measure.vars = c("bias","slope","r2","rmse","sig2","mae","mape"))

options(scipen=10000)

results_plot = results.m %>% 
                     filter(ref_name=="All") %>%
                     filter(!is.na(beam_type)) %>%
                     filter(variable=="rmse" | variable=="mape" | variable=="bias")

p_results = ggplot(results_plot, 
                      aes(x=night_flg, y=value, color=ref_tcc_bin, fill=ref_tcc_bin ,size=n
                       , group=ref_name
                       #, shape=ref_name
                       )) + 
  facet_grid(variable~seg_snow+beam_type, scales = "free_y") + 

  labs(color="ref_tcc_bin", size="# obs.",  y="Error value", x=NULL, 
       title = "Error estimates of ICESat-2 canopy height", subtitle = "ATL08 ('h_can') validation in North American boreal with airborne lidar", caption = "Airborne lidar: GLiHT 2014,2018; LVIS 2017,2019") +
  geom_text(data = results_plot %>% 
            filter(n<30) , aes(y=value, label = n), position = position_dodge(width = 1), size=3, vjust=0, hjust=0.1, color='black') +
  #scale_fill_manual(values = c("dodgerblue", "orange"))+
  #scale_color_distiller(palette = "RdYlGn", direction = 0, name="Tree canopy\ncover bin (%)")+
  #scale_fill_distiller(palette = "RdYlGn", direction = 0, guide = F)+
  scale_color_viridis_c(direction=-1, name="Tree canopy\ncover bin (%)")+
  scale_fill_viridis_c(direction=-1, guide = F)+
  scale_size(breaks = c(30,300,3000), range = c(1,8))+
  #scale_shape_manual(values = c(21,24, 22), name="Reference\nlidar")+
  guides(color = guide_legend(order = 1, override.aes = list(size=6)),
         size = guide_legend(order = 2 )#,
         #shape = guide_legend(order = 3, override.aes = list(size=5))
         )+
  geom_point( alpha=0.75, position=position_dodge(width=0.5)) + 
  theme_bw()+
  theme(text = element_text(size=14),
        plot.subtitle = element_text(size=10),
        plot.caption = element_text(size=8))
print(p_results)
ggsave(plot = p_results,
         file = path(dir_zonal,paste0("output_",format(Sys.time(), "%Y%m%d")), "p_results_error_snow_daynight_tcc_LVIS_GLiHT.png"),
         device = 'png', dpi = 300, width = 6, height = 6, units='in')
```


```{r}
#
# Model and Map: with ALL lidar: map residuals
#
frac_of_vals = c(0.01, 0.1)
# https://rstudio-pubs-static.s3.amazonaws.com/212793_f130ecc723da4ed98680d6d3d5c4aff9.html
# https://r-spatial.github.io/mapview/articles/articles/mapview_02-advanced.html
pal_tcc = rev(viridis(5)) #colorRampPalette(brewer.pal(5, "RdYlGn"))(5)

for(frac_of_val in frac_of_vals){
  

  z = data.frame(
        sample_frac(df_atl08_LIDAR_filt , frac_of_val)%>% 
          drop_na(ref_tcc_bin, ref_ht_bin) %>%
          dplyr::group_by(night_flg, seg_snow, beam_type, ref_tcc_bin) %>%
          dplyr::mutate(resid = round(residuals( lm(h_can ~ ref_ht )), dig=2)) %>%
          dplyr::mutate(resid_abs = round(abs(residuals( lm(h_can ~ ref_ht ))), dig=2)) %>%
          dplyr::mutate(resid_rel = round(100*abs(residuals( lm(h_can ~ ref_ht )))/ref_ht, dig=0))# %>%
          #filter(resid_rel < 200)
      )
  ggplot(z, aes(resid_rel, fill=ref_tcc_bin)) + scale_fill_manual(values = pal_tcc ) + geom_histogram(binwidth = 10, stat = "count", position = "stack")
  
  #
  # Plot proportion of ref_tcc_bin that is >100% error for each ref_tcc_bin
  #
  z_smry =  z %>%
    mutate(resid_rel_bin = cut( resid_rel, breaks = c(0,100,5000), include.lowest=T)) %>%
    dplyr::group_by(ref_tcc_bin, resid_rel_bin) %>%
    dplyr::summarise(n_resid = n_distinct(ref_ht))  
  
  z_smry = z_smry %>%  pivot_wider(names_from = resid_rel_bin, values_from = n_resid) 
  z_smry = data.frame(
      z_smry %>%
        mutate(prc_gt_100pct_resid_rel = unlist(100 *.[3]/rowSums(.[2:3])))
  )
    
  # Quick plot 
  ggplot(z_smry, aes(x=ref_tcc_bin, y= prc_gt_100pct_resid_rel )) + geom_point()
  
  
  #z$ref_tcc_bin = as.numeric(as.character(z$ref_tcc_bin))
  z.sf = sf::st_as_sf(z, coords = c("lon", "lat"), crs = 4326, agr = "constant")
  m = mapview(z.sf,
          zcol="ref_tcc",
          layer.name = "Tree canopy cover (%)",
          at = seq(0, 100, 20),
          
          #zcol="resid",
          #at = seq(-20,20,2),
          
          cex="resid_rel", # Adjust circle sized manually: https://stackoverflow.com/questions/56467194/plotting-point-data-in-mapview-non-linear-scale-for-point-symbol-sizes 
          #color="grey40",
          col.regions = pal_tcc,
           
          legend = TRUE,
          map.types = c("Esri.WorldShadedRelief", "Esri.WorldImagery")
          )
  ## create standalone .html
  mapshot(m, url = path(dir_zonal, paste0("df_atl08_LIDAR_filt_frac_",frac_of_val,".html")))
}
```


```{r}
# Residual boxplots
# https://stackoverflow.com/questions/41475577/saving-residuals-of-a-regression-to-the-original-dataframe?rq=1

plot_residbox_MSW = ggplot(
      df_atl08_LIDAR_filt %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, msw_flg) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))), 
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_wrap(.~msw_flg, nrow = 2) + 
    geom_hline(aes(yintercept = 0)) + 
    scale_y_continuous(limits = c(-20,20))+
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_MSW)
ggsave(plot = plot_residbox_MSW,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_MSW",ref_name,filt_type,"h_can","ref_ht",sep="_"), out_ext)),
         device = out_ext_type, dpi = 300, width = 14, height = 6)

plot_residbox_TCC = ggplot(
      df_atl08_LIDAR_filt %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, ref_tcc_bin) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))),
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_wrap(.~ref_tcc_bin, nrow = 2) + 
    geom_hline(aes(yintercept = 0)) + 
    scale_y_continuous(limits = c(-20,20))+
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_TCC)
ggsave(plot = plot_residbox_TCC,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_TCC",ref_name,filt_type,"h_can","ref_ht",sep="_"), out_ext)),
         device = out_ext_type, dpi = 300, width = 14, height = 6)

plot_residbox_TCC_gt = ggplot(
      df_atl08_LIDAR_filt %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, ref_tcc_bin, gt) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))),
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_grid(gt~ref_tcc_bin) + 
    geom_hline(aes(yintercept = 0)) + 
    scale_y_continuous(limits = c(-20,20))+
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_TCC_gt)
ggsave(plot = plot_residbox_TCC_gt,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_TCC_gt",ref_name,filt_type,"h_can","ref_ht",sep="_"), out_ext)),
         device = out_ext_type, dpi = 300, width = 16, height = 10)
```

```{r fig.height=5, fig.width=16, message=FALSE, warning=FALSE}
# Filter 2
filt_type = "DOY"
filt_str = "d > 150 & d < 255"
title_info = paste0(" [ ", filt_type, " : ", filt_str," ] ")
df_atl08_LIDAR_FILT = df_atl08_LIDAR %>% dplyr::filter_(filt_str) # June - mid-Sept

ggplot(df_atl08_LIDAR_FILT) +
    geom_point(alpha=0.2, aes(x = lon, y = lat, col=h_can), size = 1, alpha = 0.25) +
    scale_colour_gradientn(colors=pal_ht, limits=c(0,20), name = "Height (m)") +
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = p_subtitle)+
    facet_wrap(~lidar_name)+
    theme_minimal()

print(do_height_plots(df_atl08_LIDAR_FILT %>% 
                        drop_na(ref_tcc_bin, ref_ht_bin) %>%
                        filter(lidar_name=='GLiHT'), 
                      ref_name, filt_type, filt_str, "h_can", "ref_ht", dir_zonal))

print(do_height_plots(df_atl08_LIDAR_FILT %>% 
                        drop_na(ref_tcc_bin, ref_ht_bin) %>%
                        filter(lidar_name=='LVIS'), 
                      ref_name, filt_type, filt_str, "h_can", "ref_ht", dir_zonal))
```

```{r}
# Residual boxplots
# https://stackoverflow.com/questions/41475577/saving-residuals-of-a-regression-to-the-original-dataframe?rq=1

plot_residbox_MSW = ggplot(
      df_atl08_LIDAR_FILT %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, msw_flg) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))),
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_wrap(.~msw_flg, nrow = 2) + 
    geom_hline(aes(yintercept = 0)) + 
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    scale_y_continuous(limits = c(-20,20))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_MSW)
ggsave(plot = plot_residbox_MSW,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_MSW",ref_name,filt_type,"h_can","ref_ht",sep="_"), out_ext)),
         device = out_ext_type, dpi = 300, width = 14, height = 6)

plot_residbox_TCC = ggplot(
      df_atl08_LIDAR_FILT %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, ref_tcc_bin) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))),
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_wrap(.~ref_tcc_bin, nrow = 2) + 
    geom_hline(aes(yintercept = 0)) + 
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    scale_y_continuous(limits = c(-20,20))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_TCC)
ggsave(plot = plot_residbox_TCC,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_TCC",ref_name,filt_type,"h_can","ref_ht",sep="_"),out_ext)),
         device = out_ext_type, dpi = 300, width = 14, height = 6)

plot_residbox_TCC_gt = ggplot(
      df_atl08_LIDAR_FILT %>% 
        drop_na(ref_tcc_bin, ref_ht_bin) %>%
        dplyr::group_by(night_flg, ref_tcc_bin, gt) %>%
        dplyr::mutate(resid = residuals( lm(h_can ~ ref_ht ))),
    aes(x=ref_ht_bin, y=resid, fill=night_flg, outlier.color=night_flg)) + 
    labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("ATL08 height: ", "h_can"), y="Residuals", x="Reference height bin (m)") +
    geom_boxplot(outlier.size = 0.5) + 
    facet_grid(gt~ref_tcc_bin) + 
    geom_hline(aes(yintercept = 0)) + 
    scale_y_continuous(limits = c(-20,20))+
  scale_fill_manual(values = c("orange", "grey"))+
  scale_color_manual(values = c("orange", "grey"))+
    theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
print(plot_residbox_TCC_gt)
ggsave(plot = plot_residbox_TCC_gt,
         #file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("plot_residbox_TCC_gt",ref_name,filt_type,"h_can","ref_ht",sep="_"),out_ext)),
       file = paste0(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), "/", paste0(paste("plot_residbox_TCC_gt",ref_name,filt_type,"h_can","ref_ht",sep="_"),out_ext) ),
         device = out_ext_type, dpi = 300, width = 14, height = 10)
```







```{r fig.width=14}
#
# Get min dif NAME and VALUE of ATL08 rh metrics with LVIS rh98
#
# Link
#
# Note: ATL08 h_can --> rh98, h_max_can --> rh100
# https://stackoverflow.com/questions/37800877/calculate-the-minimum-difference-among-the-three-column-and-give-the-correspondi?rq=1
df_tmp = df_atl08_LIDAR_FILT %>% 
      select(contains('rh'),h_can, ref_ht)

names(df_tmp)[grep("_can", names(df_tmp))] <- c("rh98")


df_atl08_LIDAR_FILT$min_dif_value <- as.factor(apply(df_tmp, 1, function(x) {
  round(min(abs(x[1:10] - x[c(11)])), dig=1)
}))

df_atl08_LIDAR_FILT$min_dif_name <- as.factor(apply(df_tmp, 1, function(x) {
  names(x)[which.min(abs(x[1:10] - x[c(11)]))]
}))


df_atl08_LIDAR_FILT$min_dif_name = factor(df_atl08_LIDAR_FILT$min_dif_name, levels =c("rh25","rh50","rh60","rh70","rh75","rh80","rh85","rh90", "rh95","rh98")) #, "rh100"
ggplot(df_atl08_LIDAR_FILT, aes(x=as.numeric(min_dif_value))) + geom_histogram(binwidth = 0.5) + labs(title) + facet_wrap(~ref_tcc_bin)
ggplot(df_atl08_LIDAR_FILT, aes(x=min_dif_name, fill=min_dif_name)) + geom_histogram(stat = "count") + labs(title) + scale_fill_brewer(palette = "Spectral", direction=-1) + facet_wrap(~ref_tcc_bin)

labs_val_plots =  labs(
      title = "Validation of forest height from ATL08", 
      y = "Forest height (m, ATL08)",
      x = "Forest height (m, LVIS rh98)", 
      caption = "Facets: (horizontal) RH metric from LVIS of minimum dz with ATL08 estimate & (vertical) tree cover interval" ,
      subtitle = paste0("For minimum dz LVIS RH metrics across tree cover intervals ",
                        "\n# samples (", 1*100,"% of total): ", format(dim(df_atl08_LIDAR_FILT)[1], big.mark=",") )
       )
text_options = theme(
          legend.title = element_text(size = 14),
          legend.text  = element_text(size = 14), 
          axis.text    = element_text(size = 13),
          strip.text   = element_text(size = 14), 
          axis.title   = element_text(size = 14), 
          plot.title   = element_text(size = 16),
          plot.subtitle= element_text(size = 15, hjust=0, face="italic", color="black")
          )
p_scatmatrix_atl08_LIDAR = ggplot(df_atl08_LIDAR_FILT %>% drop_na(ref_tcc_bin), aes(y=h_can, x=ref_ht)) + 
  labs_val_plots +
  #geom_hex(bins=20) + 
  theme_minimal() + 
  scale_y_continuous(limits = ylims)+
  scale_x_continuous(limits = xlims)+
  #geom_point(alpha=0.05)+
  geom_bin2d(binwidth = c(1, 1))+
  geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x, size=1) +
  #geom_text_lm +
  #geom_text(data=df_eq, aes(x = quantile(xlims, probs=.5, names=F), y = quantile(ylims, probs=.8, names=F),label=paste0("n ==",n)), parse = TRUE, hjust=0.5, vjust=0, size=3)+
  scale_fill_gradientn( breaks=c(1,10,100,1000), colours=pal_rev_spectral9, trans = "log", name='# obs.')+
  geom_abline(alpha=0.5)+
  facet_grid(ref_tcc_bin~min_dif_name)+
  coord_equal()+
  text_options
print(p_scatmatrix_atl08_LIDAR)

ggsave(plot = p_scatmatrix_atl08_LIDAR,
         file = path(dir_zonal, paste0("output_",format(Sys.time(), "%Y%m%d")), paste0(paste("p_scatmatrix_atl08_LIDAR",ref_name,filt_type,"h_can","ref_ht",sep="_"),out_ext)),
         device = out_ext_type, dpi = 300, width = 8, height = 8)

# ggplot(df_atl08_LIDAR_FILT %>% drop_na(ref_tcc_bin), aes(x=ref_ht ) ) + geom_histogram() +
#     facet_grid(ref_tcc_bin~min_dif_name)+
#   text_options


```

```{r}
# Exploratory nonsense
# h_dif_ref vs asr   
ggplot(sample_frac(df_atl08_LIDAR %>% filter(abs(h_dif_ref) < 10), 0.1),           aes(x=h_can-rh95, y=h_dif_ref, color=h_dif_ref)) + geom_point(alpha=1) + scale_colour_gradient2() + facet_wrap(~ref_tcc_bin)
ggplot(sample_frac(df_atl08_LIDAR %>% filter(abs(h_dif_ref) < 10 & asr < 1), 0.1), aes(x=h_can-rh95, y=asr,       color=h_dif_ref)) + geom_point(alpha=0.2) + scale_colour_gradient2() + facet_wrap(~ref_tcc_bin)
ggplot(sample_frac(df_atl08_LIDAR %>% filter(abs(h_dif_ref) < 10 & asr < 1), 0.1), aes(x=h_can-rh95, y=h_te_unc,  color=h_dif_ref)) + geom_point(alpha=0.2) + scale_colour_gradient2() + facet_wrap(~ref_tcc_bin)
ggplot(sample_frac(df_atl08_LIDAR, 0.5),aes(x=d, fill=as.factor(msw_flg))) + geom_histogram(position="stack") + scale_colour_gradient2() + facet_wrap(~ref_tcc_bin) + coord_polar() + scale_x_continuous(breaks = seq(0,360,30))
ggplot(sample_frac(df_atl08_LIDAR, 1),aes(x=m+0.5, fill=as.factor(msw_flg))) + geom_histogram(position="stack") + scale_colour_gradient2() + facet_wrap(~yr) + coord_polar(theta = "x") + scale_x_continuous(breaks = seq(1,12,1))
```



```{r fig.height=5, fig.width=14, message=FALSE, warning=FALSE}

# COVER= TRUE
# thresh_h_can = 30
# ref_name = "GLiHT"
# p_title = paste0("Zonal stats analysis of ATL08 and ",ref_name)
# p_subtitle = "ICESat-2 ATL08 Canopy Height (m, h_can)"
# 
# 
# 
# # Sample the data
# df_atl08_gliht = sample_frac(df_atl08_gliht, 1)
# print(dim(df_atl08_gliht))
# 
# 
# 
# # Setup tcc
# df_atl08_gliht$tcc_prc_air = 100 * df_atl08_gliht$tcc1.37.median
# 
# # Make treecover classes
# cc_bin=20
# df_atl08_gliht$tcc_bin_air = cut(df_atl08_gliht$tcc_prc_air, seq(0,100,cc_bin), labels= as.character(seq(cc_bin,100,cc_bin)) )
# 
# # Base filtering
# df_atl08_gliht = df_atl08_gliht %>%
#     # Filter ATL08 data
#         dplyr::filter(n_toc_ph>0)%>%
#         dplyr::filter(h_can<thresh_h_can)%>%
#         dplyr::filter(h_te_unc!=max(h_te_unc))%>%
#         dplyr::filter(ter_slp!=max(ter_slp)) %>%
#         dplyr::filter(lat > 45)
# # Bin the reference rh100 values
# df_atl08_gliht$rh100.median_bin = cut(df_atl08_gliht$rh100.median, seq(0,30,1), labels= as.character(seq(1,30,1)) )
# 
# ####
# # Create a seequence of facet plots based on data with different filters
# 
# # Filter 1
# filt_type = "None"
# filt_str = ""
# title_info = paste0(" [ ", filt_type, " : ", filt_str," ] ")
# ggplot(df_atl08_gliht) +
#     geom_point(alpha=0.2, aes(x = lon, y = lat, col=h_can), size = 1, alpha = 0.25) +
#     scale_colour_gradientn(colors=pal_ht, limits=c(0,20), name = "Height (m)") +
#     labs(title = subset_title, subtitle = p_subtitle)
# print(dim(df_atl08_gliht))
# print(do_height_plots(df_atl08_gliht, ref_name, filt_type, filt_str, "h_can", "rh100.median", dir_zonal))
# 
# # Residual boxplots
# # https://stackoverflow.com/questions/41475577/saving-residuals-of-a-regression-to-the-original-dataframe?rq=1
# df_atl08_gliht = df_atl08_gliht %>% 
#     #dplyr::select("h_can", "rh100.median", "msw_flg","cloud_flg", "night_flg", "tcc_flg", "gt","tcc_bin_air", "rh100.median_bin" ) %>%
#     drop_na_(x_var) %>%
#     dplyr::group_by(night_flg, msw_flg) %>%
#     dplyr::mutate(resid = residuals( lm(h_can ~ rh100.median )))
# plot_residbox_MSW = ggplot(df_atl08_gliht, aes(x=rh100.median_bin, y=resid)) + 
#     labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("Reference height: ", x_var), y="Residuals", x="Reference height bin (m)") +
#     geom_boxplot() + 
#     facet_grid(night_flg~msw_flg) + 
#     geom_hline(aes(yintercept = 0)) + 
#     theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
# print(plot_residbox_MSW)
# # ggsave(plot = plot_residbox_MSW,
# #          file = path(outdir, paste0("output_",format(Sys.time(), "_%Y%m%d")), paste0(paste("plot_residbox_MSW",ref_name,filt_type,"h_can","rh100.median",sep="_"),".png")),
# #          device = 'png', dpi = 300, width = 14, height = 6)
# # Filter 2
# filt_type = "DOY"
# filt_str = "d > 150 & d < 255"
# title_info = paste0(" [ ", filt_type, " : ", filt_str," ] ")
# df_atl08_gliht_tmp = df_atl08_gliht %>% dplyr::filter_(filt_str) # June - mid-Sept
# 
# ggplot(df_atl08_gliht_tmp) +
#     geom_point(alpha=0.2, aes(x = lon, y = lat, col=h_can), size = 1, alpha = 0.25) +
#     scale_colour_gradientn(colors=pal_ht, limits=c(0,20), name = "Height (m)") +
#     labs(title = paste0(p_title, " [ ", filt_type, " : ", filt_str," ] "), subtitle = p_subtitle)
# print(dim(df_atl08_gliht_tmp))
# print(do_height_plots(df_atl08_gliht_tmp, ref_name, filt_type, filt_str, "h_can", "rh100.median", dir_zonal))
# # Residual boxplots
# # https://stackoverflow.com/questions/41475577/saving-residuals-of-a-regression-to-the-original-dataframe?rq=1
# df_atl08_gliht_tmp = df_atl08_gliht_tmp %>% 
#     #dplyr::select("h_can", "rh100.median", "msw_flg","cloud_flg", "night_flg", "tcc_flg", "gt","tcc_bin_air", "rh100.median_bin" ) %>%
#     drop_na_(x_var) %>%
#     dplyr::group_by(night_flg, msw_flg) %>%
#     dplyr::mutate(resid = residuals( lm(h_can ~ rh100.median )))
# plot_residbox_TCC = ggplot(df_atl08_gliht, aes(x=rh100.median_bin, y=resid)) +
#     labs(title = paste0("Zonal stats analysis of ATL08 and ", ref_name, " : ", title_info), subtitle = paste0("Reference height: ", x_var), y="Residuals", x="Reference height bin (m)") + 
#     geom_boxplot() + 
#     facet_grid(night_flg~tcc_bin_air) + 
#     geom_hline(aes(yintercept = 0)) + 
#     theme_minimal()+theme(axis.text = element_text(size = 6, angle = 90))#+coord_fixed()
# print(plot_residbox_TCC)
# # ggsave(plot = plot_residbox_TCC,
# #          file = path(outdir, paste0("output_",format(Sys.time(), "_%Y%m%d")), paste0(paste("plot_residbox_TCC",ref_name,filt_type,"h_can","rh100.median",sep="_"),".png")),
# #          device = 'png', dpi = 300, width = 14, height = 6)
```