Merge pull request #3274 from istfer/stics_packages

update to STICS wrappers to use stics R packages

Makefile.depends

@@ -25,7 +25,7 @@ $(call depends,models/maespa): | .install/base/logger .install/base/remote .inst
 $(call depends,models/preles): | .install/base/logger .install/base/utils .install/modules/data.atmosphere
 $(call depends,models/sibcasa): | .install/base/logger .install/base/utils
 $(call depends,models/sipnet): | .install/base/logger .install/base/remote .install/base/utils .install/modules/data.atmosphere
-$(call depends,models/stics): | .install/base/db .install/base/logger .install/base/remote .install/base/settings .install/base/utils
+$(call depends,models/stics): | .install/base/logger .install/base/remote .install/base/settings .install/base/utils
 $(call depends,models/template): | .install/base/db .install/base/logger .install/base/utils
 $(call depends,modules/allometry): | .install/base/db
 $(call depends,modules/assim.batch): | .install/base/db .install/base/logger .install/base/remote .install/base/settings .install/base/utils .install/base/workflow .install/modules/benchmark .install/modules/emulator .install/modules/meta.analysis .install/modules/uncertainty

docker/depends/pecan_deps_from_github.txt

@@ -5,3 +5,5 @@ ebimodeling/biocro@0.951
 MikkoPeltoniemi/Rpreles
 ropensci/geonames
 ropensci/nneo
+SticsRPacks/SticsOnR
+SticsRPacks/SticsRFiles

docker/depends/pecan_package_dependencies.csv

@@ -112,6 +112,7 @@
 "httr","*","modules/data.land","Imports",FALSE
 "IDPmisc","*","modules/assim.batch","Imports",FALSE
 "jsonlite","*","base/remote","Imports",FALSE
+"jsonlite","*","models/stics","Imports",FALSE
 "jsonlite","*","modules/data.atmosphere","Imports",FALSE
 "knitr",">= 1.42","base/db","Suggests",FALSE
 "knitr",">= 1.42","base/qaqc","Suggests",FALSE
@@ -267,7 +268,6 @@
 "PEcAn.DB","*","base/settings","Imports",TRUE
 "PEcAn.DB","*","base/workflow","Imports",TRUE
 "PEcAn.DB","*","models/biocro","Suggests",TRUE
-"PEcAn.DB","*","models/stics","Imports",TRUE
 "PEcAn.DB","*","models/template","Imports",TRUE
 "PEcAn.DB","*","modules/allometry","Imports",TRUE
 "PEcAn.DB","*","modules/assim.batch","Imports",TRUE
@@ -425,7 +425,6 @@
 "purrr","*","base/settings","Imports",FALSE
 "purrr","*","base/utils","Imports",FALSE
 "purrr","*","models/ed","Imports",FALSE
-"purrr","*","models/stics","Imports",FALSE
 "purrr","*","modules/assim.sequential","Imports",FALSE
 "purrr","*","modules/data.land","Imports",FALSE
 "purrr","*","modules/data.remote","Imports",FALSE
@@ -542,6 +541,7 @@
 "stats","*","modules/assim.batch","Imports",FALSE
 "stats","*","modules/assim.sequential","Suggests",FALSE
 "stats","*","modules/photosynthesis","Imports",FALSE
+"SticsRFiles","*","models/stics","Suggests",FALSE
 "stringi","*","base/logger","Imports",FALSE
 "stringi","*","base/utils","Imports",FALSE
 "stringr","*","models/fates","Imports",FALSE

models/stics/DESCRIPTION

@@ -10,17 +10,20 @@ Authors@R: c(
 Description: This module provides functions to link the STICS to PEcAn.
 Imports:
     PEcAn.settings,
-    PEcAn.DB,
     PEcAn.logger,
     PEcAn.utils (>= 1.4.8),
     PEcAn.remote,
+    jsonlite,
     lubridate,
     ncdf4,
-    purrr,
     XML,
     dplyr
 Suggests:
+    SticsRFiles,
     testthat (>= 1.0.2)
+Remotes:
+    github::SticsRPacks/SticsRFiles,
+    github::SticsRPacks/SticsOnR
 SystemRequirements: STICS
 OS_type: unix
 License: BSD_3_clause + file LICENSE

models/stics/R/met2model.STICS.R

@@ -43,8 +43,8 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
                         host = PEcAn.remote::fqdn(),
                         mimetype = "text/plain",
                         formatname = "climate",
-                        startdate = start_date,
-                        enddate = end_date,
+                        startdate = start_date, # these need fixing,not same for all climate files
+                        enddate = end_date,  # these need fixing
                         dbfile.name = out.files,
                         stringsAsFactors = FALSE)
   PEcAn.logger::logger.info("internal results")
@@ -63,6 +63,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
 
     if (file.exists(out.files.full[ctr]) && !overwrite) {
       PEcAn.logger::logger.debug("File '", out.files.full[ctr], "' already exists, skipping to next file.")
+      ctr <- ctr + 1
       next
     }
 
@@ -112,14 +113,24 @@ met2model.STICS <- 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")
 
-      dt <- PEcAn.utils::seconds_in_year(year) / length(sec)
+      dt <- diff(sec)[1]
       tstep <- round(86400 / dt)
       dt <- 86400 / tstep
 
       ind <- rep(simdays, each = tstep)
 
+      if(unlist(strsplit(nc$dim$time$units, " "))[1] %in% c("days", "day")){
+        #this should always be the case, but just in case
+        origin_dt <- (as.POSIXct(unlist(strsplit(nc$dim$time$units, " "))[3], "%Y-%m-%d", tz="UTC") + 60*60*24) - dt
+        ydays <- lubridate::yday(origin_dt + sec)
+
+      }else{
+        PEcAn.logger::logger.error("Check units of time in the weather data.")
+      }
+
       # column 6: minimum temperature (°C)
       Tair   <- ncdf4::ncvar_get(nc, "air_temperature")  ## in Kelvin
+      Tair   <- Tair[ydays %in% simdays]
       Tair_C <- PEcAn.utils::ud_convert(Tair, "K", "degC")
       t_dmin <- round(tapply(Tair_C, ind, min, na.rm = TRUE), digits = 2) # maybe round these numbers 
       weather_df[ ,6] <- t_dmin
@@ -131,23 +142,28 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
       # column 8: global radiation (MJ m-2. j-1)
       rad <- ncdf4::ncvar_get(nc, "surface_downwelling_shortwave_flux_in_air")
       gr  <- rad * 0.0864 # W m-2 to MJ m-2 d-1
+      gr  <- gr[ydays %in% simdays] 
       weather_df[ ,8] <- round(tapply(gr, ind, mean, na.rm = TRUE), digits = 2) # irradiation (MJ m-2 d-1)
 
       # column 9: Penman PET (mm.j-1) OPTIONAL, leave it as -999.9 for now
 
       # column 10: rainfall (mm.j-1)
       Rain  <- ncdf4::ncvar_get(nc, "precipitation_flux") # kg m-2 s-1
+      Rain  <- Rain[ydays %in% simdays]
       raini <- tapply(Rain * 86400, ind, mean, na.rm = TRUE) 
       weather_df[ ,10] <- round(raini, digits = 2) # precipitation (mm d-1)	
 
       # column 11: wind (m.s-1) 
       # OPTIONAL if you're not using the “Shuttleworth and Wallace” method or the “Penman calculate” method to calculate PET in the station file
       U <- try(ncdf4::ncvar_get(nc, "eastward_wind"))
       V <- try(ncdf4::ncvar_get(nc, "northward_wind"))
-      if(is.numeric(U) & is.numeric(V)){
+      if(is.numeric(U) & is.numeric(V) & !all(is.nan(U)) & !all(is.nan(V))){
+        U  <- U[ydays %in% simdays]
+        V  <- V[ydays %in% simdays]
         ws <- sqrt(U ^ 2 + V ^ 2)      
       }else{
         ws <- try(ncdf4::ncvar_get(nc, "wind_speed"))
+        ws <- ws[ydays %in% simdays]
         if (is.numeric(ws)) {
           PEcAn.logger::logger.info("eastward_wind and northward_wind absent; using wind_speed")
         }else{
@@ -161,6 +177,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
 
       # column 13: CO2 content(ppm). 
       co2 <- try(ncdf4::ncvar_get(nc, "mole_fraction_of_carbon_dioxide_in_air"))
+      co2 <- co2[ydays %in% simdays]
       if(is.numeric(co2)){
         weather_df[ ,13] <- round(tapply(co2 * 1e6, ind, mean,  na.rm = TRUE), digits = 1)
       }else{

models/stics/R/model2netcdf.STICS.R

@@ -22,14 +22,18 @@
 ##' @export
 ##'
 ##' @author Istem Fer
+##' 
 model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, overwrite = FALSE) {
-
-
+
   ### Read in model output in STICS format
   out_files <- list.files(outdir)
 
   stics_out_file <- file.path(outdir, out_files[grepl("mod_s.*", out_files)])
-  stics_output   <- utils::read.table(stics_out_file, header = T, sep = ";")
+  stics_output   <- lapply(stics_out_file, utils::read.table, header = TRUE, sep = ";")
+  stics_output  <- do.call("rbind", stics_output)
+  # probably already ordered, but order by year and DoY
+  stics_output <-  stics_output[order(stics_output[,1], stics_output[,4]), ]
+
 
   simulation_years <- unique(stics_output$ian)
 
@@ -38,7 +42,8 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o
 
   # check that specified years and output years match
   if (!all(year_seq %in% simulation_years)) {
-    PEcAn.logger::logger.severe("Years selected for model run and STICS output years do not match ")
+    # if not fail altogether, so that it won't break ensemble analysis 
+    PEcAn.logger::logger.error("Years selected for model run and STICS output years do not match.")
   }
 
   # determine time step?
@@ -52,7 +57,38 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o
     thisyear <- stics_output[ , "ian"] == y
 
     outlist <- list()
-    outlist[[1]] <- stics_output[thisyear, "lai.n."]  # LAI in (m2 m-2)
+    outlist[[length(outlist)+1]] <- stics_output[thisyear, "lai.n."]  # LAI in (m2 m-2)
+
+    # daily amount of CO2-C emitted due to soil mineralisation (humus and organic residues) (kg ha-1 d-1)
+    HeteroResp <- PEcAn.utils::ud_convert(stics_output[thisyear, "CO2sol"],  "ha-1 day-1", "m-2 s-1")
+
+    outlist[[length(outlist)+1]] <- HeteroResp
+
+
+    # dltams(n): daily growth rate of the plant (t.ha-1.d-1) 
+    dltams      <- PEcAn.utils::ud_convert(stics_output[thisyear, "dltams.n."], "ton", "kg") * 0.48 # ton to kgC
+    # dltaremobil: daily amount of perennial reserves remobilised (t.ha-1.d-1)
+    dltaremobil <- PEcAn.utils::ud_convert(stics_output[thisyear, "dltaremobil"], "ton", "kg") * 0.48 # ton to kgC
+
+    NPP <- dltams - dltaremobil # kgC ha-1 d-1
+    NPP[NPP<0] <- 0
+
+    # double checking that this is all NPP (above and below)
+    ## this:
+    #stics_output[thisyear, "dltams.n."]  # t.ha-1.d-1
+    ## should be roughly equal to this:
+    #diff(stics_output[thisyear, "masec.n."])+ diff(stics_output[thisyear, "msrac.n."]) # t.ha-1
+
+    NPP <- PEcAn.utils::ud_convert(NPP, "ha-1 day-1", "m-2 s-1") # kg C m-2 s-1  
+    outlist[[length(outlist)+1]] <- NPP
+
+    NEE <- -1*(NPP-HeteroResp)
+    outlist[[length(outlist)+1]] <- NEE
+
+    # other vars
+    # Cr: amount of C in organic residues mixed with soil (kg.ha-1)
+    # Crac: amount of C in roots at harvest (kg.ha-1)
+    # Chumt: amount of C in humified organic matter (active + inert fractions) (kg.ha-1)
 
     # ******************** Declare netCDF dimensions and variables ********************#
     t <- ncdf4::ncdim_def(name = "time", 
@@ -68,7 +104,10 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o
     dims <- list(lon = lon, lat = lat, time = t)
 
     nc_var <- list()
-    nc_var[[1]] <- PEcAn.utils::to_ncvar("LAI", dims)
+    nc_var[[length(nc_var)+1]] <- PEcAn.utils::to_ncvar("LAI", dims)
+    nc_var[[length(nc_var)+1]] <- PEcAn.utils::to_ncvar("HeteroResp", dims)
+    nc_var[[length(nc_var)+1]] <- PEcAn.utils::to_ncvar("NPP", dims)
+    nc_var[[length(nc_var)+1]] <- PEcAn.utils::to_ncvar("NEE", dims)
 
     # ******************** Declare netCDF variables ********************#