In mikejohnson51/geogrids: Create mosaiced geogrids with flexable structure

knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE)

hydro_path       = params$hydrofabric_file
data_dir         = params$data_dir
hydrofabric_file = params$hydrofabric_file

library(geogrids)
library(data.table)
library(sf)
library(dplyr)
library(zonal)
library(terra)
library(ggplot2)
library(foreach)
library(doParallel)

out_csv   <-  file.path(params$output_dir, paste0(params$version,"-hydrofabric_traits.csv"))

make_plot = function(traits, geom, var, varname, colors, outfile){
 df = merge(geom, traits) %>% 
   select(var  = !!var)

 gg = ggplot(df) + 
   geom_sf(aes(fill = var), color = NA) + 
     scale_fill_gradientn(colours=colors, guide = "colourbar") +
     labs(fill = varname) 

  ggsave(gg, filename = outfile)
}

cats  <- read_sf(hydro_path, "catchments")
conus <- AOI::aoi_get(state = "conus", union = TRUE) %>% 
  st_transform(st_crs(cats))

geom  = st_filter(cats, conus, .predicate = st_within)
years = params$years

files = geo_cache_list(data_dir)

Load PPT, TAVG, and Snow Thresholds

pr = filter(files, grepl('pr', fullname)) %>% 
     filter(grepl(".nc", fullname)) %>% 
     arrange(desc(fullname)) %>% 
     slice(1:years) %>% 
     arrange(fullname)

tm = filter(files, grepl('tavg', fullname))  %>%
     filter(grepl(".nc", fullname)) %>% 
     arrange(desc(fullname)) %>% 
     slice(1:years) %>% 
     arrange(fullname)

gridmet_w = zonal::weighting_grid(pr$fullname[1], geom, "ID")

pr_list = foreach::foreach(i = 1:nrow(pr), .combine = cbind) %dopar% { execute_zonal(file = pr$fullname[i], w = gridmet_w)[,-c("ID")] }
ppt = setnames(pr_list, paste0('day', 1:ncol(pr_list)))

tm_list = foreach::foreach(i = 1:nrow(tm), .combine = cbind) %dopar% { execute_zonal(file = tm$fullname[i], w = gridmet_w)[,-c("ID")] }
tavg = setnames(tm_list, paste0('day', 1:ncol(tm_list)))


traits = data.frame(ID = sort(geom$ID))

Mean PPT

meanPPT = rowMeans(ppt)

traits$meanPPT = meanPPT

if(params$plot){
  make_plot(traits, geom, 
            'meanPPT', "Mean PPT", 
            colors = pals::warmcool(100), outfile = "img/meanPPT.png")
}

jennings = filter(files, grepl("jennings", fullname))

snow_tif = geogrid_warp(jennings$fullname, make_grid(tm$fullname[1]))

jen = execute_zonal(file = snow_tif, w = gridmet_w)
jen$kelvins = jen$V1 + 273.15

snow_day = tavg[, .SD  < jen$kelvins]
snow = ppt * snow_day

traits$snowFrac  = rowSums(snow) / rowSums(ppt)

if(params$plot){
  make_plot(traits, geom, 
            'snowFrac', "Snow Fraction",
            colors = pals::warmcool(100), outfile = "img/snowfrac.png")
}

High PPT Freq

ppt5x   = 5 * meanPPT
highPPT = ppt[, .SD > ppt5x]
traits$high_ppt_freq = rowSums(highPPT) / years

lowPPT = ppt[, .SD < 1]
traits$low_ppt_freq = rowSums(lowPPT) / years

if(params$plot){
  make_plot(traits, geom, 
            'low_ppt_freq', "Low PPT Freq",
            colors = pals::coolwarm(100), outfile = "img/low_ppt_freq.png")

  make_plot(traits, geom, 
            'high_ppt_freq', "High PPT Freq",
            colors = pals::warmcool(100), outfile = "img/high_ppt_freq.png")

}

high <- data.table(t(highPPT))
mod_cols <- names(high) 

ff = function(x){
    cs = cumsum(x)
    cs = cs - cummax((x == 0) * cs)
    c(ifelse(diff(cs) < 0, cs, NA), cs[length(cs)])
}

high[ , (mod_cols) := lapply(.SD, ff), .SDcols = mod_cols]  

low<- data.table(t(lowPPT))
low[ , (mod_cols) := lapply(.SD, ff), .SDcols = mod_cols]  

traits$high_ppt_dur = colMeans(high, na.rm = TRUE)
traits$low_ppt_dur = colMeans(low, na.rm = TRUE)

if(params$plot){

  make_plot(traits, geom, 
            'high_ppt_dur', "High Duration PPT",
            colors = pals::warmcool(100), outfile = "img/high_dur_ppt.png")

  make_plot(traits, geom, 
            'low_ppt_dur', "Low Duration PPT",
            colors = pals::coolwarm(100), outfile = "img/low_dur_ppt.png")
}

Catchment

DEM = filter(files, grepl("elevation", fullname))$fullname
DEM = terra::rast(DEM)
t = c(terra::terrain(DEM), DEM)
terrain = execute_zonal(t, geom, 'ID') %>% 
  setNames(c("ID", "elevation", "slope"))

terrain$areasqkm = as.numeric(st_area(geom)/1e6)

traits = left_join(traits, terrain, by = "ID")

if(params$plot){
  make_plot(traits, geom, 
            'elevation', "Mean Elevation",
            colors = pals::parula(100), outfile = "img/mean_elevation.png")

  make_plot(traits, geom, 
            'slope', "Mean Slope",
            colors = pals::tol.rainbow(100), outfile = "img/mean_slope.png")
}

Soil

s = filter(files, grepl("sand-1m|silt-1m|clay-1m|rockdepm|GLIM_xx|permeability_permafrost", fullname)) %>% 
  filter(grepl(".tif", fullname)) %>% 
  dplyr::pull(fullname) %>% 
  terra::rast()

soils_w = weighting_grid(s$`clay-1m-percent`, geom, "ID")

soils = execute_zonal(s, w = soils_w) %>% 
  setNames(c("ID", names(s)))

soils = soils %>% 
  mutate(k = -0.60 + 
           (0.0126*`sand-1m-percent`) - 
           (0.0064*`clay-1m-percent`),
         porosity = 50.5 - 
           (0.142*`sand-1m-percent`)  - 
           (0.037*`clay-1m-percent`),
         mwc = porosity * rockdepm,
         carbonates = GLIM_xx,
         GLIM_xx = NULL)

traits = left_join(traits, soils, by = "ID")

if(params$plot){
  make_plot(traits, geom, 
            'clay-1m-percent', "Clay Percentage",
            colors = pals::parula(100), outfile = "img/clay.png")

  make_plot(traits, geom, 
            'sand-1m-percent', "Sand Percentage",
            colors = pals::parula(100), outfile = "img/sand.png")

  make_plot(traits, geom, 
            'silt-1m-percent', "Silt Percentage",
            colors = pals::parula(100), outfile = "img/silt.png")

  make_plot(traits, geom, 
            'rockdepm', "Rock Depth",
            colors = pals::parula(100), outfile = "img/rockdepth.png")

  make_plot(traits, geom, 
            'permeability_permafrost', "Permeability",
            colors = pals::parula(100), outfile = "img/permeability.png")

  make_plot(traits, geom, 
            'carbonates', "Carbonates",
            colors = pals::parula(100), outfile = "img/carbonates.png")

  make_plot(traits, geom, 
            'k', "K",
            colors = pals::parula(100), outfile = "img/k.png")

  make_plot(traits, geom, 
            'porosity', "Porosity",
            colors = pals::parula(100), outfile = "img/porosity.png")

  make_plot(traits, geom, 
            'mwc', "Max Water Content",
            colors = pals::parula(100), outfile = "img/mwc.png")
}

Energy

energy = filter(files, grepl("ai_normal|pet_normal", fullname)) %>% 
  dplyr::pull(fullname) %>% 
  terra::rast()

energy = execute_zonal(energy, w = gridmet_w) %>% 
  setNames(c("ID", names(energy)))

traits = left_join(traits, energy, by = "ID")

if(params$plot){
  make_plot(traits, geom, 
            'ai_normal_2001-2020', "AI",
            colors = pals::jet(100), outfile = "img/ai.png")

  make_plot(traits, geom, 
            'pet_normal_2001-2020', "PET",
            colors = pals::jet(100), outfile = "img/pet.png")
}

Landuse

lc = filter(files, grepl('MCD12Q1.006/1km/2019-01-01.tif', fullname))

lc_tiff = rast(lc$fullname)

lu = zonal::execute_zonal_cat(lc_tiff, w = soils_w)

forest = lu %>% 
  filter(value %in% c(2:5)) %>% 
  group_by(ID) %>% 
  summarise(forest = sum(percentage)) %>% 
  ungroup()

traits = left_join(traits, forest, by = "ID")

if(params$plot){
  make_plot(traits, geom, 
            'forest', "Forest",
            colors = pals::brewer.greens(100), outfile = "img/forest.png")
}

GVF & LAI

modis_mapping = read.csv('..//modis_lc.csv')

modis_mapping$ndvi_inf = 
c(0.81,0.86,0.88,0.90,0.87,0.86,0.86,0.75,
  0.76,0.74,0.86,0.84,0.86,0.82,NA, 0.86,NA)

indices<-rep(1:12,times=years)

ndvi_denom = classify(lc_tiff, select(modis_mapping, is = Class, becomes = ndvi_inf)) - 0.05
########
ndvi = filter(files, grepl('MOD13A3.006/mosaics/landcover', fullname)) %>% 
  arrange(desc(fullname)) %>% 
  slice(1:(12*years)) %>% 
  arrange(fullname)

# apply scaling factor
ndvi_rast = (rast(ndvi$fullname) / 10000.0) 

gvf = (ndvi_rast - 0.05) / ndvi_denom
gvf.mean<-terra::tapp(gvf, indices, fun = mean)
gvf.max = max(gvf.mean)
gvf.min = min(gvf.mean)

#########
lai_files = filter(files, grepl('MOD15A2H.006/monthly_means', fullname)) %>% 
  arrange(desc(fullname)) %>% 
  slice(1:(12*years)) %>% 
  arrange(fullname)

lai_rast = rast(lai_files$fullname) * 0.1

lai.mean<-tapp(lai_rast, indices, fun = mean)
lai.max = max(lai.mean)
lai.min = min(lai.mean)

gvf_lai_rast = rast(list(gvf_max = gvf.max, 
                     gvf_diff = gvf.max-gvf.min,
                     lai_max = lai.max, 
                     lai_diff = lai.max-lai.min))

gvf_lai_stat = zonal::execute_zonal(gvf_lai_rast, w = soils_w) %>% 
  setNames(c("ID", names(gvf_lai_rast)))

traits = left_join(traits, gvf_lai_stat, by = "ID")

if(params$plot){
  make_plot(traits, geom, 
            gvf_lai_stat$gvf_max, "Max GVF",
            colors = pals::brewer.greens(100), outfile = "img/max_gvf.png")

  make_plot(traits, geom, 
            gvf_lai_stat$gvf_diff, "GVF Differnece",
            colors = pals::brewer.greens(100), outfile = "img/dif_gvf.png")

  make_plot(traits, geom, 
            gvf_lai_stat$lai_max, "Max LAI",
            colors = pals::brewer.greens(100), outfile = "img/max_lai.png")

  make_plot(traits, geom, 
            gvf_lai_stat$lai_diff, "LAI Differnece",
            colors = pals::brewer.greens(100), outfile = "img/dif_lai.png")
}

Output

data.table::fwrite(traits, out_csv)

if(!is.null(params$AWS_bucket)){

  aws.s3::put_object(
      file   = out_csv, 
      object = basename(out_csv), 
      bucket = params$AWS_bucket,
      multipart = TRUE)
}