PEPRMT: new parameters, met integration

.gitignore

@@ -81,6 +81,7 @@ tests/BC*
 compile_on_geo.sh
 documentation/tutorials/*/*.html
 pecan.Rproj
+.rcookies
 shiny/BenchmarkReport/*
 # Installation files from make
 .install/

models/peprmt/R/hls2model.PEPRMT.R

@@ -26,7 +26,7 @@ hls2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
   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/example_data/Field_Boundary.geojson")
+  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'

models/peprmt/R/met2model.PEPRMT.R

@@ -22,27 +22,13 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
 
   PEcAn.logger::logger.info("START met2model.PEPRMT")
 
-  ## PEPRMT requires the following inputs: 
-  ## Time_2 = d[,1] # day of year (1-infinite # of days)
-  #DOY_disc_2=d[,2] #discontinuous day of year that starts over every year (1-365 or 366)
-  #Year_2=d[,3] #year 
-  #TA_2 = d[,4] #Air temperature (C)
-  #WT_2 = d[,5] #water table depth (cm) equals 0 when water table at soil surface
-  #PAR_2 <- d[,6] #photosynthetically active radiation (umol m-2 d-1)
-  #LAI_2 <- d[,7] #Leaf area index (if not using can be 0s or NaN)
-  #GI_2 <- d[,8] #greeness index from Phenocam (GCC) or Landsat EVI etc (unitless)
-  #FPAR <- d[,9] #If using LAI data, set FPAR variable to 1's, if using a greenness index set FPAR to 0's
-  #LUE<- d[,10] #growing season LUE computed for each site using measured GPP in gC per umol
-  #wetland_age_2= d[,11] #Age of wetland in years
-  #Sal <- d[,12] #Salinity (ppt)
-  #NO3 <- d[,13] #Dissolved NO3 (mg/L)
-  #SOM_2 <-d[,14] #Decomposed Organic matter : all the decomposed soil organic matter in top meter of soil informed buy MEM inclusive of current year
-  #site_2 <-d[,15] #Site: if running more than 1 site, have 1s in this column for first site, 2s for 2nd site and so on
-
   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,".")
+  if(nchar(in.prefix)>0 & 
+     !substr(in.prefix, nchar(in.prefix), nchar(in.prefix)) %in% c(".", "_")) {
+    in.prefix = paste0(in.prefix,".")
+    }
 
   out.file <- paste0(in.prefix, start_date_string,".",
                     strptime(end_date, "%Y-%m-%d"),
@@ -77,11 +63,10 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
   end_year <- lubridate::year(end_date)
 
   ## Loop through and add air temp for each year
-  for (year in start_year:end_year) {
-    print(year)
-    Year_2=year #year 
+  for (Year in start_year:end_year) {
+    print(Year)
 
-    old.file <- file.path(in.path, paste(in.prefix, year, ".nc", sep = ""))
+    old.file <- file.path(in.path, paste(in.prefix, Year, ".nc", sep = ""))
     if(!file.exists(old.file)) PEcAn.logger::logger.error("file not found",old.file)
     ## open netcdf
     nc <- ncdf4::nc_open(old.file)
@@ -90,7 +75,7 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
     sec <- nc$dim$time$vals
     sec <- PEcAn.utils::ud_convert(sec, unlist(strsplit(nc$dim$time$units, " "))[1], "seconds")
     timestep.s <- 86400  # seconds in a day
-    dt <- PEcAn.utils::seconds_in_year(year) / length(sec)
+    dt <- PEcAn.utils::seconds_in_year(Year) / length(sec)
     tstep <- round(timestep.s / dt) #4 per day
     dt    <- timestep.s / tstep  #dt is now an integer
 
@@ -102,21 +87,21 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
     ncdf4::nc_close(nc)
 
     ## build day of year
-    diy <- PEcAn.utils::days_in_year(year)
+    diy <- PEcAn.utils::days_in_year(Year)
     doy <- rep(seq_len(diy), each = timestep.s / dt)[seq_along(sec)]
 
     ## Aggregate variables up to daily
-    TA_2         <- PEcAn.utils::ud_convert(tapply(Tair, doy, mean, na.rm = TRUE), "Kelvin", "Celsius")
+    TA_C            <- PEcAn.utils::ud_convert(tapply(Tair, doy, mean, na.rm = TRUE), "Kelvin", "Celsius")
     #Consider using a different met source for PAR (since this only has SW)
-    PAR_2          <- tapply(2.114 * SW * dt, doy, sum, na.rm = TRUE) / (24*60*60) 
+    PAR_umol_m2_day <- tapply(2.114 * SW * dt, doy, sum, na.rm = TRUE) / (24*60*60) 
     #https://rdrr.io/cran/LakeMetabolizer/man/sw.to.par.html
-    Time_2          <- tapply(doy, doy, mean)
+    DOY_disc        <- tapply(doy, doy, mean)
 
     ## build data matrix
-    tmp <- cbind(Time_2, PAR_2, TA_2)
+    tmp <- cbind(DOY_disc, PAR_umol_m2_day, TA_C, Year)
 
     ##filter out days not included in start or end date
-    if(year == start_year){
+    if(Year == start_year){
       start.row <- length(as.Date(paste0(start_year, "-01-01")):as.Date(start_date)) #extra days length includes the start date
       if (start.row > 1){
         PEcAn.logger::logger.info("Subsetting PEPRMT met to match start date ", as.Date(start_date))
@@ -125,8 +110,8 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
         tmp <- tmp[start.row:nrow(tmp),]
       }
     }
-    if (year == end_year){
-      if(year == start_year){
+    if (Year == end_year){
+      if(Year == start_year){
         end.row <- length(as.Date(start_date):as.Date(end_date))
         if (end.row < nrow(tmp)){
           PEcAn.logger::logger.info("Subsetting PEPRMT met to match end date")
@@ -151,45 +136,14 @@ met2model.PEPRMT <- function(in.path, in.prefix, outfolder, start_date, end_date
 
   ## Assuming default values for some variables
   Dates = seq.Date(as.Date(start_date), as.Date(end_date), by = "1 day")
-  Time_2 = as.integer(Dates - as.Date(start_date))+1
-  DOY_disc_2 = lubridate::yday(Dates) #discontinuous day of year that starts over every year (1-365 or 366)
-
-  # TO DO: currently missing all of the following drivers
-  WT_2 = NA
-  LAI_2 = NA
-  GI_2 = NA
-  FPAR = NA
-  LUE = NA
-  wetland_age_2 = NA
-  Sal = NA
-  NO3 = NA
-  SOM_2 = NA
-  site_2 = 1
+  DOY = as.integer(Dates - as.Date(start_date))+1
 
   #column order matters
-  final <- cbind(out, DOY_disc_2, Year_2, TA_2, WT_2, LAI_2, GI_2, FPAR, 
-                 LUE, wetland_age_2, Sal, NO3, SOM_2, site_2) %>%
+  final <- cbind(out, DOY, Year) %>%
     data.frame() %>%
-    select(all_of(c("Time_2", "DOY_disc_2", "Year_2", "TA_2", "WT_2", 
-                    "PAR_2", "LAI_2", "GI_2", "FPAR", 
-                    "LUE", "wetland_age_2", "Sal", "NO3", "SOM_2", "site_2")))
-
-  if(sum(is.na(final$WT_2))!=0) warning("Warning: missing water table depth. PEPRMT won't run.")
-  if(sum(is.na(final$GI_2))!=0 & 
-     sum(is.na(final$LAI_2))!=0) warning("Warning: missing greenness index and LAI both. PEPRMT won't run.")
-  if(sum(is.na(final$FPAR))!=0) warning("Warning: missing FPAR. PEPRMT won't run.")
-  if(sum(is.na(final$PAR))!=0) warning("Warning: missing PAR. PEPRMT won't run.")
-  if(sum(is.na(final$LUE))!=0) warning("Warning: missing LUE. PEPRMT won't run.")
-  if(sum(is.na(final$wetland_age_2))!=0) {
-    warning("Warning: missing wetland age. Age assumed to be > 2yr.")
-    final$wetland_age_2 <- 1000
-  }
-  if(sum(is.na(final$Sal))!=0) warning("Warning: missing salinity. PEPRMT won't run.")
-  if(sum(is.na(final$NO3))!=0) warning("Warning: missing NO3. PEPRMT won't run.")
-  if(sum(is.na(final$SOM_2))!=0) warning("Warning: missing SOM. PEPRMT won't run.")
-  if(sum(is.na(final$site_2))!=0) warning("Warning: missing site. Site set to 1.")
+    select(all_of(c("Year", "DOY_disc",  "DOY", "TA_C", "PAR_umol_m2_day")))
 
-  utils::write.table(final, out.file.full, quote = FALSE, sep = " ", row.names = FALSE, col.names = FALSE)
+  utils::write.table(final, out.file.full, quote = FALSE, sep = " ", row.names = FALSE, col.names = TRUE)
 
   return(invisible(results))
 

models/peprmt/R/model2netcdf.PEPRMT.R

@@ -34,11 +34,8 @@ model2netcdf.PEPRMT <- function(outdir, sitelat, sitelon, start_date, end_date,
     sub.PEPRMT.output.dims <- dim(sub.PEPRMT.output)
 
     # ******************** Declare netCDF variables ********************#
-    start.day <- 1
-    if (y == lubridate::year(start_date)){
-      start.day <- lubridate::yday(start_date)
-    } 
-    tvals <- (start.day:sub.PEPRMT.output.dims[1])-1
+    start.day <- min(sub.PEPRMT.output$DOY_disc)
+    tvals <- (start.day:((start.day + sub.PEPRMT.output.dims[1]) - 1))
     bounds <- array(data=NA, dim=c(length(tvals),2))
     bounds[,1] <- tvals
     bounds[,2] <- bounds[,1]+1
@@ -74,8 +71,8 @@ model2netcdf.PEPRMT <- function(outdir, sitelat, sitelon, start_date, end_date,
     output[[4]] <- (sub.PEPRMT.output[, "NEE_mod"] * 0.001)   # NEE in kgC/m2/s
 
     ## Pools
-    output[[5]]  <- (sub.PEPRMT.output[, fluxes[1]] * 0.001)  # Soil Carbon, kgC/m2
-    output[[6]]  <- (sub.PEPRMT.output[, fluxes[2]] * 0.001)  # Soil Carbon, kgC/m2
+    output[[5]]  <- (sub.PEPRMT.output[, pools[1]] * 0.001)  # Soil Carbon, kgC/m2
+    output[[6]]  <- (sub.PEPRMT.output[, pools[2]] * 0.001)  # Soil Carbon, kgC/m2
 
     ## time_bounds
     output[[length(fluxes) +
@@ -99,9 +96,11 @@ model2netcdf.PEPRMT <- function(outdir, sitelat, sitelon, start_date, end_date,
     nc_var[[5]]  <- PEcAn.utils::to_ncvar("slow_soil_pool_carbon_content", dims)
     nc_var[[6]]  <- PEcAn.utils::to_ncvar("fast_soil_pool_carbon_content", dims)
 
-    nc_var[[7]] <- ncdf4::ncvar_def(name="time_bounds", units='', 
-                                     longname = "history time interval endpoints", dim=list(time_interval,time = t), 
-                                     prec = "double")
+    nc_var[[7]] <- ncdf4::ncvar_def(name="time_bounds", 
+                                    units='', 
+                                    longname = "history time interval endpoints", 
+                                    dim=list(time_interval,time = t), 
+                                    prec = "double")
 
     ### Output netCDF data
     nc <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var)

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

@@ -66,14 +66,20 @@ write.config.PEPRMT <- function(defaults, trait.values, settings, run.id) {
       )
     }
   }
-  GPP_names <- c("GPP_a0", "GPP_a1", "GPP_Ha", "GPP_Hd")
-  Reco_names <- c("Reco_Ea_som", "Reco_kM_som",
-                  "Reco_Ea_labile", "Reco_kM_labile")
-  CH4_names <- paste("CH4", 1:8, sep = "_") # TODO get meaningful names from Patty
-  missing_traitnames <- setdiff(c(GPP_names, Reco_names, CH4_names), trait_names)
+
+  params <- c("wetland_type", "a0", "a1", "Ha", "Hd", "T_opt_GPP", "Ea_SOM", 
+              "kM_SOM", "Ea_labile", "kM_labile", "Ea_SOM_CH4", "kM_SOM_CH4", 
+              "Ea_labile_CH4", "kM_labile_CH4", "Ea_oxi_CH4", "kM_oxi_CH4", 
+              "kI_SO4", "kI_NO3", "k_plant_oxi")
+  provided_traitnames <- intersect(params, trait_names)
+  if (!"wetland_type" %in% provided_traitnames) {
+    trait_values["wetland_type"] <- 2
+    provided_traitnames <- c(provided_traitnames, "wetland_type")
+  }
+  missing_traitnames <- setdiff(params, trait_names)
   if (length(missing_traitnames) > 0) {
-    PEcAn.logger::logger.error(
-      "Parameters missing from trait.values",
+    PEcAn.logger::logger.warn(
+      "Parameters missing from trait.values. Will use default",
       sQuote(missing_traitnames)
     )
   }
@@ -94,15 +100,29 @@ write.config.PEPRMT <- function(defaults, trait.values, settings, run.id) {
 
   jobsh <- gsub("@BINARY@", settings$model$binary, jobsh)
   jobsh <- gsub("@DELETE_RAW@", as.logical(settings$model$delete.raw %||% FALSE), jobsh)
+
+  param_str <- paste0(
+    "list(",
+    paste(
+      paste0(provided_traitnames, " = ", trait_values[provided_traitnames]),
+      collapse = ", "
+    ),
+    ")"
+  )
+
+  jobsh <- gsub("@PARAMS@", param_str, jobsh, fixed = TRUE)
 
-  jobsh <- gsub("@GPP_THETA@", paste(trait_values[GPP_names], collapse = ", "), jobsh)
-  jobsh <- gsub("@RECO_THETA@", paste(trait_values[Reco_names], collapse = ", "), jobsh)
-  jobsh <- gsub("@CH4_THETA@", paste(trait_values[CH4_names], collapse = ", "), jobsh)
-
-
-  # yes, this will be replaced with real params once demo is working
-  run_data <- PEPRMT::example_data |>
-    dplyr::filter(.data$site == settings$run$site$id)
+  # MET
+  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")]
+
+  peprmt_specific_input_path <- settings$run$inputs$PEPRMT$path
+
+  run_data <- read.csv(peprmt_specific_input_path) |> 
+    dplyr::select(-any_of(met_vars)) |>
+    dplyr::right_join(met) |>
+    dplyr::mutate(site = settings$run$site$id)
 
   utils::write.csv(run_data, file.path(rundir, "run_data.csv"), row.names = FALSE)
   writeLines(jobsh, con = file.path(rundir, "job.sh"))

models/peprmt/demo_run/.gitignore

@@ -0,0 +1,4 @@
+input_demo_out
+out
+output
+.Rapp.history