EVI download and processing code for PEPRMT model

CHANGELOG.md

@@ -34,6 +34,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Added `PEcAn.data.land::to_co2e()` for converting SOC change, CH4, and N2O to CO2-equivalent emissions using IPCC Global Warming Potential values.
 - Added `PEcAn.data.land::event_parquet_to_json` for generating PEcAn `event.json` files from well-formatted event parquet files, with support for ensembles of events.
 - Added statewide synthetic fertilization and compost amendment event workflows for CA ag parcels. Outputs share an ensemble naming so a downstream cleaner unions them into one fertilization event type for SIPNET.
+- Added `PEcAn.remote::get_EVI_HLS()` for loading harmonized landsat-sentinel data and calculating EVI
 
 ### Fixed
 - The median run's manifest row went in with the literal strings `"NA"` for pft and trait, which `read.csv` turns into real `NA`, so `read.sa.output` never matched the median quantile and `splinefun` silently dropped that knot. Sensitivity analysis output changes as a result: partial variances shift slightly, though rankings are unaffected in the cases checked.

models/peprmt/R/hls2model.PEPRMT.R

@@ -14,177 +14,93 @@
 ##' @param ... additional arguments, currently ignored
 ##' @author Abigail Lewis (add names)
 hls2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date,
+                             edl_username = Sys.getenv("edl_username"), 
+                             edl_password = Sys.getenv("edl_password"),
                              lat, lon,overwrite = FALSE, verbose = FALSE, ...) {
 
   PEcAn.logger::logger.info("START hls2model.PEPRMT")
 
   # Error if any of the soft dependencies we use here are not installed
   PEcAn.utils::need_packages(c("terra", "earthdatalogin", "dygraphs", "imager", "rstac", "xts"))
 
-  #https://github.com/nasa/HLS-Data-Resources/blob/main/r/HLS_Tutorial.Rmd
-  earthdatalogin::edl_netrc()
-  s = rstac::stac("https://cmr.earthdata.nasa.gov/stac/LPCLOUD/")
-  HLS_col <- list("HLSS30_2.0", "HLSL30_2.0")
-  #roi <- terra::vect(data.frame(latitude = lat, longitude = lon))
-  roi <- terra::vect("models/peprmt/demo_run/raw-data/Field_Boundary.geojson")
-  roi_extent <- terra::ext(roi)
-  bbox <- c(roi_extent$xmin, roi_extent$ymin, roi_extent$xmax, roi_extent$ymax)
-  roi_datetime <- '2021-08-01T00:00:00Z/2021-09-30T23:59:59Z'
-  items <- s %>%
-    rstac::stac_search(collections = HLS_col,
-                bbox = bbox,
-                datetime = roi_datetime,
-                limit = 100) %>%
-    rstac::post_request()
-  assets <- rstac::items_assets(items)
-  print(assets)
-  browse_image_url <- items$features[[1]]$assets$browse$href
-  browse <-imager::load.image(browse_image_url)
-  plot(browse)
-  sf_items <- rstac::items_as_sf(items)
-  # Retrieve Granule ID for each feature
-  granule_id <- sapply(items$features, function(feature) feature$id)
-  # Add as first column in sf_items
-  sf_items <- cbind(granule = granule_id, sf_items)
-  # Define a function to extract asset urls for selected bands
-  # # This also includes a check to ensure the correct bands are extracted
-  # # # depending on the collection (HLSL30 or HLSS30)
-  extract_asset_urls <- function(feature) {
-    collection_id <- feature$collection
-    if (collection_id == "HLSS30_2.0") {
-      bands = c('B02', 'B8A','B04','Fmask')
-    } else if (collection_id == "HLSL30_2.0") {
-      bands = c('B02', 'B05','B04','Fmask')}
-    sapply(bands, function(band) feature$assets[[band]]$href)
-  }
-  # Retrieve Asset URLs for each feature using our extract_asset_urls function and transpose them to columns
-  asset_urls <- t(sapply(items$features, extract_asset_urls))
-
-  colnames(asset_urls) <- c('blue', 'nir', 'red', 'fmask')
-  sf_items <- cbind(sf_items, asset_urls)
-
-  # Filter based on cloud cover
-  sf_items <- sf_items[sf_items$eo.cloud_cover < 30,]
-  # Reset Row Indices
-  row.names(sf_items) <- NULL
-
-  terra::setGDALconfig("GDAL_HTTP_UNSAFESSL", value = "YES")
-  terra::setGDALconfig("GDAL_HTTP_COOKIEFILE", value = ".rcookies")
-  terra::setGDALconfig("GDAL_HTTP_COOKIEJAR", value = ".rcookies")
-  terra::setGDALconfig("GDAL_DISABLE_READDIR_ON_OPEN", value = "EMPTY_DIR")
-  terra::setGDALconfig("CPL_VSIL_CURL_ALLOWED_EXTENSIONS", value = "TIF")
-
-  # This function reads an HLS scene from a URL, applies the scale factor if necessary, and optionally crops and
-  # masks the scene based on a polygon. It requries the above GDAL configurations and a .netrc file. A .netrc
-  # can be created by running `earthdatalogin::edl_netrc()`.
-  open_hls <- function(url, roi = NULL) {
-    # Add VSICURL prefix
-    url <- paste0('/vsicurl/', url)
-    # Retrieve metadata
-    meta <- terra::describe(url)
-    # Check if dataset is Quality Layer (Fmask) - no scaling this asset (int8 datatype)
-    is_fmask <- any(grep("Fmask", meta))
-    # Check if Scale is present in band metadata
-    will_autoscale <- any(grep("Scale:", meta))
-    # Read the raster
-    r <- terra::rast(url)
-    # Apply Scale Factor if necessary
-    if (!will_autoscale && !is_fmask){
-      print(paste("No scale factor found in band metadata. Applying scale factor of 0.0001 to", basename(url)))
-      r <- r * 0.0001
-    }
-    # Crop if roi specified
-    if (!is.null(roi)){
-      # Reproject roi to match crs of r
-      roi_reproj <- terra::project(roi, terra::crs(r))
-      r <- terra::mask(terra::crop(r, roi_reproj), roi_reproj)
-    }
-    return(r)
-  }
-
-  # Test opening and crop
-  red <- open_hls(sf_items$red[1], roi)
-  terra::plot(red)
-
-  red_stack <- lapply(sf_items$red, open_hls, roi = roi)
-  nir_stack <- lapply(sf_items$nir, open_hls, roi = roi)
-  blue_stack <- lapply(sf_items$blue, open_hls, roi = roi)
-  fmask_stack <- lapply(sf_items$fmask, open_hls, roi = roi)
-
-  calculate_EVI <- function(nir, red, blue){
-    evi <- 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1.0)
-    return(evi)
-  }
-
-  # Calculate EVI For All of our Scenes
-  evi_stack <- mapply(calculate_EVI, nir_stack, red_stack, blue_stack, SIMPLIFY = FALSE)
-
-  # Rename the Scenes in our List
-  names(evi_stack) <- sf_items$datetime
-  # Create a single Rast Object from our list
-  evi_stacks <- terra::rast(evi_stack)
-
-  # Test fmask
-  fmask <- fmask_stack[[23]]
-  plot(fmask)
-
-  selected_bit_nums <- c(1,2,3,4,5)
+  # Format dates for file names
+  start_date <- as.POSIXlt(start_date, tz = "UTC")
+  start_date_string <- as.character(strptime(start_date, "%Y-%m-%d"))
+  end_date <- as.POSIXlt(end_date, tz = "UTC")
+
+  out.file <- paste0(in.prefix,
+                     ".csv")
+  out.file.full <- file.path(outfolder, out.file)
 
-  # Filter based on quality
-  build_mask <- function(fmask, selected_bit_nums){
-    # Create a mask of all zeros
-    mask <- terra::rast(fmask, vals=0)
-    for (b in selected_bit_nums){
-      # Apply Bitwise AND to fmask values and selected bit numbers
-      mask_temp <- terra::app(fmask, function(x) bitwAnd(x, bitwShiftL(1,b)) >0)
-      # Update Mask to maintain only 1 layer with bitwise OR
-      mask <- mask | mask_temp
-    }
-    return(mask)
-  }
+  roi <- terra::vect(in.path)
 
-  qmask <- build_mask(fmask[[1]], selected_bit_nums)
-  plot(qmask)
+  df <- get_HLS_EVI(edl_username,
+                    edl_password,
+                    start_date,
+                    end_date,
+                    roi)
 
-  # Create List of Masks
-  qmask_stack <- lapply(fmask_stack, build_mask, selected_bit_nums=selected_bit_nums)
+  fits_beck <- df |>
+    dplyr::mutate(year = lubridate::year(Date)) |>
+    dplyr::group_by(year) |>
+    calc_curve_beck() |>
+    dplyr::bind_rows()
 
-  # Apply Mask to EVI using NA Values
-  evi_masked <- mapply(function(x, y) {
-    terra::mask(x, y, maskvalue = TRUE, updatevalue = NA)
-  }, evi_stack, qmask_stack, SIMPLIFY = FALSE)
+  out <- fits_beck |>
+    dplyr::mutate(Year = lubridate::year(Date),
+                  DOY_disc = doy)|>
+    dplyr::rename(EVI = pred) |>
+    dplyr::select(Year, DOY_disc, EVI) 
 
-  evi_masked <- terra::rast(evi_masked)
+  utils::write.table(out, out.file.full, quote = FALSE, sep = " ", row.names = FALSE, col.names = TRUE)
 
-  # Add Date Only Column
-  sf_items$date <- sapply(sf_items$datetime, function(x) strsplit(x, "T")[[1]][1])
+  } # hls2model.PEPRMT
 
-  evi_mean <- terra::global(evi_masked, 'mean', na.rm=TRUE)
-  evi_max <- terra::global(evi_masked, 'max', na.rm=TRUE)
-  evi_min <- terra::global(evi_masked, 'min', na.rm=TRUE)
-  evi_sd <- terra::global(evi_masked, 'sd', na.rm=TRUE)
-
-  stats <- data.frame(
-    EVI_Max = evi_max,
-    EVI_Min = evi_min,
-    EVI_mean = evi_mean,
-    EVI_SD = evi_sd
-  )
-  stats$Date <- lubridate::ymd_hms(sf_items$datetime) # convert string to date format (ISO 8601)
-  variables = xts::xts(x=stats[,-5], order.by=stats$Date) # Choose the cols with the variables
-  dygraphs::dygraph(variables) %>%
-    dygraphs::dyAxis("y",label = "EVI")
-
-  start_date <- as.POSIXlt(start_date, tz = "UTC")
-  start_date_string <- as.character(strptime(start_date, "%Y-%m-%d"))
-  end_date <- as.POSIXlt(end_date, tz = "UTC")
-  if(nchar(in.prefix)>0 & substr(in.prefix,nchar(in.prefix),nchar(in.prefix)) != ".") in.prefix = paste0(in.prefix,".")
-
-  out.file <- paste0(in.prefix, start_date_string,".",
-                     strptime(end_date, "%Y-%m-%d"),
-                     ".dat")
-  out.file.full <- file.path(outfolder, out.file)
-
-  utils::write.csv(stats, out.file.full)
 
-  } # hls2model.PEPRMT
+#' Wrapper function for double log calculation
+#'
+#' @param df data frame
+#'
+#' @returns modeled evi on each date
+#' @export
+calc_curve_beck <- function(df) {
+  year_to_run <- unique(lubridate::year(df$Date))
+  if(length(year_to_run) >1){
+    stop("EVI has to be calculated one year at a time")
+  }
+  # Add explicit NAs
+  x <- df |>
+    dplyr::mutate(Date = as.Date(Date),
+                  img_doy = lubridate::yday(Date)) |>
+    dplyr::group_by(img_doy) |>
+    dplyr::summarise(
+      Date = min(Date),
+      evi = mean(evi, na.rm = TRUE),
+      .groups = "drop"
+    ) |>
+    tidyr::complete(img_doy = 1:365)
+
+  # Run double log function
+  fit <- FitDoubleLogBeck(x$evi, t = x$img_doy, hessian = T, ninit = 100)
+
+  # Format output
+  out <- data.frame(
+    param_name = names(fit$params),
+    param_value = fit$params,
+    stdError = fit$stdError
+  )
+  rownames(out) <- NULL
+
+  pred_df_beck <- data.frame(doy = rep(1:365)) |>
+    dplyr::cross_join(out |>
+                        tidyr::pivot_wider(names_from = param_name,
+                                           values_from = c(param_value, stdError))) |>
+    dplyr::mutate(pred = param_value_mn + (param_value_mx - param_value_mn) *
+                    (1/(1 + exp(-param_value_rsp * (doy - param_value_sos))) +
+                       1/(1 + exp(param_value_rau * (doy - param_value_eos))))) |>
+    dplyr::left_join(df |> dplyr::rename(doy = img_doy)) |>
+    dplyr::mutate(method = "Beck",
+                  Date = as.Date(paste0(year_to_run,"-01-01"))+ lubridate::days(doy-1))
+
+  return(pred_df_beck)
+}

models/peprmt/R/write.config.PEPRMT.R

@@ -117,12 +117,19 @@ write.config.PEPRMT <- function(defaults, trait.values, settings, run.id) {
   met_path <- settings$run$inputs$met$path
   met <- utils::read.table(met_path, header = T)
   met_vars <- colnames(met)[!colnames(met) %in% c("Year", "DOY_disc")]
+
+  # EVI
+  evi_path <- settings$run$inputs$evi$path
+  evi <- utils::read.table(evi_path, header = T)
+  evi_vars <- colnames(evi)[!colnames(evi) %in% c("Year", "DOY_disc")]
 
   peprmt_specific_input_path <- settings$run$inputs$PEPRMT$path
 
   run_data <- utils::read.csv(peprmt_specific_input_path) |>
     dplyr::select(-dplyr::any_of(met_vars)) |>
     dplyr::right_join(met) |>
+    dplyr::select(-dplyr::any_of(evi_vars)) |>
+    dplyr::right_join(evi) |>
     dplyr::mutate(site = settings$run$site$id)
 
   # Event handling

models/peprmt/demo_run/02_input_demo.qmd

@@ -116,6 +116,12 @@ settings <- settings_raw |>
     path = file.path("data", "PEPRMT_specific_inputs"),
     path_template = "{path}/{id}_formatted.csv"
   ) |>
+  setEnsemblePaths(
+    n_reps = 1,
+    input_type = "evi",
+    path = file.path("data", "evi"),
+    path_template = "{path}/{id}_evi.csv"
+  ) |>
   # helper functions from peprmt/R/helpers.R to address having multiple grid 
   # cells and date ranges in the same site_info file
   papply(id2grid) |>
@@ -178,9 +184,58 @@ for(i in 1:nrow(site_info)){
 }
 ```
 
+## EVI
+
+To provide EVI input data for each run, we can either manually load remote sensing data and interpolate to a daily timestep, or just use provided inputs. Here is code to do this manually, which requires wetland boundaries in `raw-data/boundaries` and login information for NASA earth data. This is not run by default because it takes a long time to run
+
+#```{r load-EVI}
+library(PEcAn.remote)
+
+# Set up folder for formatted EVI data
+dir.create(here::here("models","peprmt","demo_run","data","evi"))
+
+# Fill in processed files
+for(i in 1:nrow(site_info)){
+  for(ens in 1:ensemble_size){
+    in.path <- here::here("models","peprmt","demo_run","raw-data","boundaries",
+                          paste0(site_info$id[i], ".geojson"))
+    in.prefix <- paste0(site_info$id[i], "_evi")
+    start_date <- site_info$met.start[i]
+    end_date <- site_info$met.end[i]
+    outfolder <- paste0("data/evi")
+
+    # Format and save files
+    result <- hls2model.PEPRMT(in.path, in.prefix, outfolder, 
+                               start_date, end_date,
+                               edl_username = Sys.getenv("edl_username"), 
+                               edl_password = Sys.getenv("edl_password"))
+  }
+}
+#```
+
+Alternatively, we can load pre-processed files from raw-data/site-data-formatted.
+
+```{r}
+dir.create(here::here("models","peprmt","demo_run","data","evi"))
+
+for(i in 1:nrow(site_info)){
+  in.path <- here::here("models","peprmt","demo_run","raw-data","site-data-formatted",
+                        paste0(site_info$id[i], "_formatted.csv"))
+
+  data <- read_csv(in.path) %>%
+    dplyr::select(Year, DOY_disc, EVI)
+
+  outfolder <- here::here("models","peprmt","demo_run","data", "evi",
+                        paste0(site_info$id[i], "_evi.csv"))
+
+  readr::write_csv(data, outfolder)
+}
+```
+
+
 ## Other inputs
 
-PEPRMT requires a variety of inputs that we are currently supplying manually. These input data are stored in `raw-data/site-data-formatted`, so all we need to do is copy then over to the `data` folder, which is used in our demo run.
+PEPRMT requires a variety of inputs that we are currently supplying manually (specifically, depth, salinity, nitrate, wetland age, soil organic matter, and light use efficiency). These input data are stored in `raw-data/site-data-formatted`, so all we need to do is copy then over to the `data` folder, which is used in our demo run.
 
 ```{r load-other-inputs}
 # Set up folder

models/peprmt/demo_run/raw-data/template.xml

@@ -56,6 +56,10 @@
     <id></id>
     <path><!-- inserted at config time --></path>
    </PEPRMT>
+   <evi>
+    <id></id>
+    <path><!-- inserted at config time --></path>
+   </evi>
   </inputs>
   <start.date><!-- inserted at config time --></start.date>
   <end.date><!-- inserted at config time --></end.date>

models/peprmt/inst/xml_build.R

@@ -34,6 +34,12 @@ options <- list(
       "Should contain subdirs named by site id"
     )
   ),
+  optparse::make_option("--evi_dir",
+                        default = "data/evi",
+                        help = paste(
+                          "Directory containing csvs of PEPRMT-formatted evi data."
+                        )
+  ),
   optparse::make_option("--peprmt_input_dir",
     default = "data/PEPRMT_specific_inputs",
     help = paste(
@@ -146,6 +152,12 @@ settings <- settings_init |>
     d1 = "DATES-HERE",
     path_template = "{path}/ERA5_{id}/ERA5.{n}.{d1}.dat"
   ) |>
+  setEnsemblePaths(
+    n_reps = 1,
+    input_type = "evi",
+    path = file.path("data", "evi"),
+    path_template = "{path}/{id}_evi.csv"
+  ) |>
   setEnsemblePaths(
     n_reps = 1,
     input_type = "PEPRMT",