From 8265b29a3623d95105dfb1212e001f43bac7f83b Mon Sep 17 00:00:00 2001 From: istfer Date: Tue, 26 Oct 2021 13:20:16 +0300 Subject: [PATCH 01/68] registry file for site level CF format FO inputs --- .../inst/registration/register.FieldObservatory.xml | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 modules/data.atmosphere/inst/registration/register.FieldObservatory.xml diff --git a/modules/data.atmosphere/inst/registration/register.FieldObservatory.xml b/modules/data.atmosphere/inst/registration/register.FieldObservatory.xml new file mode 100644 index 00000000000..6c22d38e184 --- /dev/null +++ b/modules/data.atmosphere/inst/registration/register.FieldObservatory.xml @@ -0,0 +1,10 @@ + + +site + + 33 + CF Meteorology + application/x-netcdf + nc + + From 8cab98a1a3c27574f94f4573201b4acf0e6f21a6 Mon Sep 17 00:00:00 2001 From: istfer Date: Tue, 26 Oct 2021 16:11:04 +0300 Subject: [PATCH 02/68] update met2model --- models/stics/R/met2model.STICS.R | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/models/stics/R/met2model.STICS.R b/models/stics/R/met2model.STICS.R index eba287815bb..8f125ff5e34 100644 --- a/models/stics/R/met2model.STICS.R +++ b/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") @@ -112,14 +112,24 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date, sec <- nc$dim$time$vals sec <- udunits2::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 <- udunits2::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,12 +141,14 @@ 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) @@ -145,9 +157,12 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date, U <- try(ncdf4::ncvar_get(nc, "eastward_wind")) V <- try(ncdf4::ncvar_get(nc, "northward_wind")) if(is.numeric(U) & is.numeric(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 +176,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{ From 11568e2caacbd5462d8697a090057ef266206396 Mon Sep 17 00:00:00 2001 From: istfer Date: Thu, 6 Jan 2022 14:14:11 +0200 Subject: [PATCH 03/68] first case of set_param_xml --- models/stics/R/met2model.STICS.R | 1 + models/stics/R/write.config.STICS.R | 21 ++++++++++++++++----- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/models/stics/R/met2model.STICS.R b/models/stics/R/met2model.STICS.R index 8f125ff5e34..c7b3452eea8 100644 --- a/models/stics/R/met2model.STICS.R +++ b/models/stics/R/met2model.STICS.R @@ -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 } diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 83cc7b747cc..81f7cdcfcce 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -59,6 +59,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) usm_list <- XML::xmlToList(usm_xml) + # NOTE: it's the text files, not the xml files that are read by the STICS executable. ################################# Prepare Plant File ####################################### @@ -66,34 +67,44 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template plt file, has all the formalisms plt_xml <- XML::xmlParse(system.file("crop_plt.xml", package = "PEcAn.STICS")) - plt_list <- XML::xmlToList(plt_xml) + #plt_list <- XML::xmlToList(plt_xml) for (pft in seq_along(trait.values)) { pft.traits <- unlist(trait.values[[pft]]) pft.names <- names(pft.traits) + plant_file <- file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")) + + # save the template + XML::saveXML(plt_xml, file = plant_file) + # go over each formalism and replace params following the order in crop_plt # for now I vary only one parameter under roots. + # TODO: vary more params # plant name and group # effect of atmospheric CO2 concentration # phasic development # emergence and starting # leaves + # radiation interception + # to see parameters per formalism + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "radiation interception") + # unlist(values) + # shoot biomass growth # partitioning of biomass in organs # yield formation # roots + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "roots") - # specific root length (cm g-1) - # plt_list[[10]][[6]][[2]][[4]] position + # longsperac - specific root length (cm g-1) if ("SRL" %in% pft.names) { srl_val <- udunits2::ud.convert(pft.traits[which(pft.names == "SRL")], "m", "cm") - plt_list <- plt_list %>% purrr::modify_depth(-1, ~if(all(.x == "@longsperac@")) srl_val else .x) - + SticsRFiles::set_param_xml(plant_file, "longsperac", srl_val, overwrite = TRUE) } # frost From 24d5c203893e3593f88b0180ed5b483f9ad20e41 Mon Sep 17 00:00:00 2001 From: istfer Date: Thu, 6 Jan 2022 14:33:18 +0200 Subject: [PATCH 04/68] xml to txt --- models/stics/R/write.config.STICS.R | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 81f7cdcfcce..78efd82b165 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -59,6 +59,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) usm_list <- XML::xmlToList(usm_xml) + # stics and javastics path + stics_path <- settings$model$binary + javastics_path <- gsub("bin","", dirname(stics_path)) + # NOTE: it's the text files, not the xml files that are read by the STICS executable. ################################# Prepare Plant File ####################################### @@ -113,13 +117,20 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # correspondance code BBCH # cultivar parameters + # delete - SticsRFiles::set_param_xml overwrites already # write back + # if(names(trait.values)[pft] != "env"){ + # + # XML::saveXML(PEcAn.settings::listToXml(plt_list, "fichierplt"), + # file = file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")), + # prefix = '\n') + # + # } + + # convert xml2txt if(names(trait.values)[pft] != "env"){ - - XML::saveXML(PEcAn.settings::listToXml(plt_list, "fichierplt"), - file = file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")), - prefix = '\n') - + SticsRFiles::convert_xml2txt(xml_file = plant_file, java_dir = javastics_path) + # do I also need to move the file out of the plant folder to main rundir? } } # pft-loop ends @@ -347,7 +358,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { } # stics path - stics_path <- settings$model$binary + # stics_path <- settings$model$binary # symlink to binary file.symlink(stics_path, bindir) @@ -357,6 +368,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { usm_name <- defaults$pft$name + # if this script can already create the txts, bypass this step cmd_generate <- paste("java -jar", jexe,"--generate-txt", rundir, usm_name) # copy *.mod files From 409100df7453a6ca0f6446ca50025db5e5f80284 Mon Sep 17 00:00:00 2001 From: istfer Date: Tue, 11 Jan 2022 16:06:39 +0200 Subject: [PATCH 05/68] vary more parameters --- models/stics/R/write.config.STICS.R | 110 ++++++++++++++++++++++++---- 1 file changed, 97 insertions(+), 13 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 78efd82b165..8aba27e4861 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -80,26 +80,96 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { plant_file <- file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")) - # save the template + # save the template, will be overwritten below XML::saveXML(plt_xml, file = plant_file) + # to learn the parameters in a plant file + # SticsRFiles::get_param_info(file_path = plant_file) + # go over each formalism and replace params following the order in crop_plt - # for now I vary only one parameter under roots. # TODO: vary more params # plant name and group # effect of atmospheric CO2 concentration + # phasic development + # to see parameters per formalism + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "phasic development") + # unlist(values) + + # minimum temperature below which development stops (degree C) + if ("tdmin" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "tdmin", pft.traits[which(pft.names == "tdmin")], overwrite = TRUE) + } + + # maximum temperature above which development stops (degree C) + if ("tdmax" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "tdmax", pft.traits[which(pft.names == "tdmax")], overwrite = TRUE) + } + + # emergence and starting # leaves + # phyllotherme, thermal duration between the apparition of two successive leaves on the main stem (degree day) + # assuming this is the same as phyllochron + if ("phyllochron" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "phyllotherme", pft.traits[which(pft.names == "phyllochron")], overwrite = TRUE) + } + + # minimal density above which interplant competition starts (m-2) + if ("dens_comp" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "bdens", pft.traits[which(pft.names == "dens_comp")], overwrite = TRUE) + } + + # LAI above which competition between plants starts (m2 m-2) + if ("lai_comp" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "laicomp", pft.traits[which(pft.names == "lai_comp")], overwrite = TRUE) + } + + # basal height of crop (m) + if ("height" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "hautbase", pft.traits[which(pft.names == "height")], overwrite = TRUE) + } + + # maximum height of crop + if ("HTMAX" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "hautmax", pft.traits[which(pft.names == "HTMAX")], overwrite = TRUE) + } + # radiation interception - # to see parameters per formalism # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "radiation interception") - # unlist(values) # shoot biomass growth + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "shoot biomass growth") + + # optimal temperature (1/2) for plant growth + if ("teopt" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "teopt", pft.traits[which(pft.names == "teopt")], overwrite = TRUE) + } + + # optimal temperature (2/2) for plant growth + if ("teoptbis" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "teoptbis", pft.traits[which(pft.names == "teoptbis")], overwrite = TRUE) + } + # partitioning of biomass in organs + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "partitioning of biomass in organs") + + # maximum SLA (specific leaf area) of green leaves (cm2 g-1) + if ("SLAMAX" %in% pft.names) { + slamax <- pft.traits[which(pft.names == "SLAMAX")] + slamax <- udunits2::ud.convert(udunits2::ud.convert(slamax, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 + SticsRFiles::set_param_xml(plant_file, "slamax", slamax, overwrite = TRUE) + } + + # minimum SLA (specific leaf area) of green leaves (cm2 g-1) + if ("SLAMIN" %in% pft.names) { + slamin <- pft.traits[which(pft.names == "SLAMIN")] + slamin <- udunits2::ud.convert(udunits2::ud.convert(slamin, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 + SticsRFiles::set_param_xml(plant_file, "slamin", slamin, overwrite = TRUE) + } + # yield formation # roots @@ -115,17 +185,31 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # water # nitrogen # correspondance code BBCH + # cultivar parameters + # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "cultivar parameters") - # delete - SticsRFiles::set_param_xml overwrites already - # write back - # if(names(trait.values)[pft] != "env"){ - # - # XML::saveXML(PEcAn.settings::listToXml(plt_list, "fichierplt"), - # file = file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")), - # prefix = '\n') - # - # } + # there are multiple cultivars (varietes) in plt file + # for now I assume we will always use #1 in simulations and modify its parameter values + # hence, _tec file will always say variete==1, if you change the logic don't forget to update handling of the _tec file accordingly + + # maximal lifespan of an adult leaf expressed in summation of Q10=2 (2**(T-Tbase)) + if ("leaf_lifespan_max" %in% pft.names) { + # this will modify the first variete by default + SticsRFiles::set_param_xml(plant_file, "durvieF", pft.traits[which(pft.names == "leaf_lifespan_max")], overwrite = TRUE) + # see example for setting the Grindstad cultivar param + # SticsRFiles::set_param_xml(plant_file, "durvieF", pft.traits[which(pft.names == "leaf_lifespan_max")], select = "Grindstad", overwrite = TRUE) + } + + # cumulative thermal time between the stages LEV (emergence) and AMF (maximum acceleration of leaf growth, end of juvenile phase) + if ("cum_thermal_juvenile" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "stlevamf", pft.traits[which(pft.names == "cum_thermal_juvenile")], overwrite = TRUE) + } + + # cumulative thermal time between the stages AMF (maximum acceleration of leaf growth, end of juvenile phase) and LAX (maximum leaf area index, end of leaf growth) + if ("cum_thermal_growth" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "stamflax", pft.traits[which(pft.names == "cum_thermal_growth")], overwrite = TRUE) + } # convert xml2txt if(names(trait.values)[pft] != "env"){ From 7b99903b7f21f7a6ce35403941eccd58e1d433fa Mon Sep 17 00:00:00 2001 From: istfer Date: Tue, 11 Jan 2022 16:52:53 +0200 Subject: [PATCH 06/68] handle until ini file --- models/stics/R/write.config.STICS.R | 48 ++++++++++++++++------------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 8aba27e4861..4a25fe41dbc 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -62,8 +62,13 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # stics and javastics path stics_path <- settings$model$binary javastics_path <- gsub("bin","", dirname(stics_path)) + + + # Per STICS development team, there are two types of STICS inputs + # Global input: _plt.xml, param_gen.xml, param_newform.xml + # Local input: _ini.xml (initialization), sols.xml (soils), _tec.xml (crop management), (climate files) _sta.xml, *.year - # NOTE: it's the text files, not the xml files that are read by the STICS executable. + # NOTE: however, it's the text files, not the xml files that are read by the STICS executable. ################################# Prepare Plant File ####################################### @@ -220,6 +225,20 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { } # pft-loop ends + ############################## Param gen / newform #################################### + + ## DO NOTHING FOR NOW + gen_xml <- XML::xmlParse(system.file("param_gen.xml", package = "PEcAn.STICS")) + gen_file <- file.path(rundir, "param_gen.xml") + XML::saveXML(gen_xml, file = gen_file) + SticsRFiles::convert_xml2txt(xml_file = gen_file, java_dir = javastics_path) + # may delete the xml after this + + newf_xml <- XML::xmlParse(system.file("param_newform.xml", package = "PEcAn.STICS")) + newf_file <- file.path(rundir, "param_newform.xml") + XML::saveXML(newf_xml, file = newf_file) + SticsRFiles::convert_xml2txt(xml_file = newf_file, java_dir = javastics_path) + ############################ Prepare Initialization File ################################## @@ -228,15 +247,16 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template ini file ini_xml <- XML::xmlParse(system.file("pecan_ini.xml", package = "PEcAn.STICS")) - ini_list <- XML::xmlToList(ini_xml) + ini_file <- file.path(rundir, paste0(defaults$pft$name, "_ini.xml")) - # DO NOTHING FOR NOW + # write the ini file + XML::saveXML(ini_xml, file = ini_file) - # write the ini file - XML::saveXML(PEcAn.settings::listToXml(ini_list, "initialisations"), - file = file.path(rundir, paste0(defaults$pft$name, "_ini.xml")), - prefix = '\n') + # DO NOTHING FOR NOW + # but when you do note that this also has multiple options, e.g. + # SticsRFiles::set_param_xml(xml_file = ini_file, param_name = "lai0", param_value = 1, select = "plante", value = "1", overwrite = TRUE) + SticsRFiles::convert_xml2txt(xml_file = ini_file, java_dir = javastics_path) ############################ Prepare Soils ################################## @@ -285,20 +305,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## skipping for now - ############################## Param gen / newform #################################### - - ## DO NOTHING - gen_xml <- XML::xmlParse(system.file("param_gen.xml", package = "PEcAn.STICS")) - gen_list <- XML::xmlToList(gen_xml) - XML::saveXML(PEcAn.settings::listToXml(gen_list, "fichierpar"), - file = file.path(rundir, "param_gen.xml"), - prefix = '\n') - - newf_xml <- XML::xmlParse(system.file("param_newform.xml", package = "PEcAn.STICS")) - newf_list <- XML::xmlToList(newf_xml) - XML::saveXML(PEcAn.settings::listToXml(newf_list, "fichierparamgen"), - file = file.path(rundir, "param_newform.xml"), - prefix = '\n') ############################ Prepare Technical File ################################## From 05a9e24223b1e5737b3feb91ef3247ef668a7301 Mon Sep 17 00:00:00 2001 From: istfer Date: Wed, 12 Jan 2022 16:56:22 +0200 Subject: [PATCH 07/68] until weather --- models/stics/R/write.config.STICS.R | 44 +++++++++++++++++------------ models/stics/inst/param.sol | 8 ++++++ 2 files changed, 34 insertions(+), 18 deletions(-) create mode 100644 models/stics/inst/param.sol diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 4a25fe41dbc..91a8e254091 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -54,10 +54,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { file = file.path(cfgdir, "preferences.xml"), prefix = '\n\n') - # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. - # TODO: more than one usm - usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) - usm_list <- XML::xmlToList(usm_xml) # stics and javastics path stics_path <- settings$model$binary @@ -85,8 +81,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { plant_file <- file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")) - # save the template, will be overwritten below - XML::saveXML(plt_xml, file = plant_file) + if(names(trait.values)[pft] != "env"){ + # save the template, will be overwritten below + XML::saveXML(plt_xml, file = plant_file) + } # to learn the parameters in a plant file # SticsRFiles::get_param_info(file_path = plant_file) @@ -195,14 +193,14 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # values = SticsRFiles::get_param_xml(plant_file, select = "formalisme", value = "cultivar parameters") # there are multiple cultivars (varietes) in plt file - # for now I assume we will always use #1 in simulations and modify its parameter values + # for now I assume we will always use only #1 in simulations # hence, _tec file will always say variete==1, if you change the logic don't forget to update handling of the _tec file accordingly # maximal lifespan of an adult leaf expressed in summation of Q10=2 (2**(T-Tbase)) if ("leaf_lifespan_max" %in% pft.names) { - # this will modify the first variete by default + # this will modifies all varietes' durvieFs by default SticsRFiles::set_param_xml(plant_file, "durvieF", pft.traits[which(pft.names == "leaf_lifespan_max")], overwrite = TRUE) - # see example for setting the Grindstad cultivar param + # see example for setting a particular (the Grindstad) cultivar param # SticsRFiles::set_param_xml(plant_file, "durvieF", pft.traits[which(pft.names == "leaf_lifespan_max")], select = "Grindstad", overwrite = TRUE) } @@ -263,19 +261,25 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## this is where we modify soil characteristics - # read in template sols file - sols_xml <- XML::xmlParse(system.file("sols.xml", package = "PEcAn.STICS")) - sols_list <- XML::xmlToList(sols_xml) + #### THERE IS SOME BUG IN SticsRFiles::convert_xml2txt FOR SOLS.XML + #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol - sols_list$sol$.attrs[["nom"]] <- paste0("sol", defaults$pft$name) + # read in template sols file (xml) + # sols_xml <- XML::xmlParse(system.file("sols.xml", package = "PEcAn.STICS")) + sols_file <- file.path(rundir, "param.sol") - # DO NOTHING FOR NOW + # cp template sols file (txt) + file.copy(system.file("param.sol", package = "PEcAn.STICS"), sols_file) - # write the tec file - XML::saveXML(PEcAn.settings::listToXml(sols_list, "sols"), - file = file.path(rundir, "sols.xml"), - prefix = '\n') + # check param names + # sols_vals <- SticsRFiles::get_soil_txt(sols_file) + SticsRFiles::set_soil_txt(filepath = sols_file, param="typsol", value=paste0("sol", defaults$pft$name)) + + # DO NOTHING ELSE FOR NOW + + # this has some bug for sols.xml + # SticsRFiles::convert_xml2txt(xml_file = sols_file, java_dir = javastics_path) ######################### Prepare Weather Station File ############################### @@ -375,6 +379,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ################################ Prepare USM file ###################################### + # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. + usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) + usm_list <- XML::xmlToList(usm_xml) + # TODO: more than 1 USM and PFTs (STICS can run 2 PFTs max: main crop + intercrop) # pft name diff --git a/models/stics/inst/param.sol b/models/stics/inst/param.sol new file mode 100644 index 00000000000..de82c8ed781 --- /dev/null +++ b/models/stics/inst/param.sol @@ -0,0 +1,8 @@ + 1 SOLPFT 17.0 0.1500 15.0000 1.0000 6.5000 0.0100 0.2000 4.0000 0.0000 200.0000 50.0000 0.5000 60.0000 5.0000 0.0100 0.0000 0.3300 + 1 2 2 2 2 2 1 1 + 1 10.0000 0.0000 0.0000 0.0000 0.0000 1.0000 10 8.0000 + 1 20.00 16.00 4.00 1.20 0.00 1 50.00 5 + 1 20.00 15.30 9.30 1.30 0.00 1 50.00 2 + 1 60.00 15.30 9.30 1.30 0.00 1 50.00 1 + 1 0.00 0.00 0.00 0.00 0.00 1 50.00 1 + 1 0.00 0.00 0.00 0.00 0.00 1 50.00 1 \ No newline at end of file From 5786115660f468b0648668ebc6685ee50e1c061d Mon Sep 17 00:00:00 2001 From: istfer Date: Wed, 12 Jan 2022 17:55:48 +0200 Subject: [PATCH 08/68] until tec --- models/stics/R/write.config.STICS.R | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 91a8e254091..251c29dcffb 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -263,6 +263,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### THERE IS SOME BUG IN SticsRFiles::convert_xml2txt FOR SOLS.XML #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol + #### TODO: revise other to have txt templates directly in hte package # read in template sols file (xml) # sols_xml <- XML::xmlParse(system.file("sols.xml", package = "PEcAn.STICS")) @@ -287,24 +288,20 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_list <- XML::xmlToList(sta_xml) + sta_file <- file.path(rundir, "pecan_sta.xml") + + XML::saveXML(sta_xml, file = sta_file) + SticsRFiles::convert_xml2txt(xml_file = sta_file, java_dir = javastics_path) + + sta_txt <- file.path(rundir, "station.txt") # change latitute - sta_list[[1]][[3]]$text <- settings$run$site$lat + SticsRFiles::set_station_txt(sta_txt, param = "latitude", value = settings$run$site$lat) # DO NOTHING ELSE FOR NOW - # Should these be prepared by met2model.STICS? - - # write the sta file - XML::saveXML(PEcAn.settings::listToXml(sta_list, "fichiersta"), - file = file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")), - prefix = '\n') - - - ############################## Prepare LAI forcing #################################### ## skipping for now From 8fb2ef3b64cf7e7ed1294a96015c14dc53cd3c01 Mon Sep 17 00:00:00 2001 From: istfer Date: Tue, 25 Jan 2022 14:15:25 +0200 Subject: [PATCH 09/68] first pass at tec file --- models/stics/R/write.config.STICS.R | 99 ++++++++++++++++------------- 1 file changed, 54 insertions(+), 45 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 251c29dcffb..b0405ce4043 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -263,7 +263,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### THERE IS SOME BUG IN SticsRFiles::convert_xml2txt FOR SOLS.XML #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol - #### TODO: revise other to have txt templates directly in hte package + #### TODO: revise others to have txt templates directly in the package # read in template sols file (xml) # sols_xml <- XML::xmlParse(system.file("sols.xml", package = "PEcAn.STICS")) @@ -312,66 +312,75 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## this is where we modify management practices - # read in template tec file - tec_xml <- XML::xmlParse(system.file("pecan_tec.xml", package = "PEcAn.STICS")) - tec_list <- XML::xmlToList(tec_xml) + ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml + tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases - # If harvest file is given, use given dates - # this will need more complicated checks and file formats + SticsRFiles::gen_tec_xml(out_path = rundir, param_table = tec_df) + + # the following formalisms exist in the tec file: + ## supply of organic residus + ## soil tillage + ## sowing + ## phenological stages + ## irrigation + ## fertilisation + ## harvest + ## special techniques + ## soil modification by techniques (compaction-fragmentation) + + # if a harvest file is given, most (all?) of our harvest cases are actually fall under special techniques - cut crop if(!is.null(settings$run$inputs$harvest)){ h_days <- as.matrix(utils::read.table(settings$run$inputs$harvest$path, header = TRUE, sep = ",")) - # probably should use nicer list manipulation techniques, but these template files are static for now - # save last list to add at the end - attr_sublist <- tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[4]] - tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[4]] <- NULL - - list_no <- 2 + # param names + h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d + "hautcoupe" , # cut height for forage crops, m + "lairesiduel", # residual LAI after each cut of forage crop, m2 m-2 + "msresiduel" , # residual aerial biomass after a cut of a forage crop, t.ha-1 + "anitcoupe") # amount of mineral N added by fertiliser application at each cut of a forage crop, kg.ha-1 + + harvest_list <- list() for(hrow in seq_len(nrow(h_days))){ - - harvest_list <- tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[2]] # refreshing the "template" - intervention_names <- names(tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[2]]) - - # If given harvest date is within simulation days - if(as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01")) %in% dseq){ - - # STICS needs cutting days in cumulative julian days - # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) - # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 - # the following code should give that - harvest_list$colonne$text <- which(dseq == as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01"))) + lubridate::yday(dseq[1]) - 2 - - tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[list_no]] <- harvest_list - - list_no <- list_no + 1 - } - - } - - # this means we have prescribed more than 2 cutting days - if(list_no > 4){ - names(tec_list[[8]][[1]][[1]][[2]][[2]][[1]])[3:(list_no-1)] <- "intervention" + # empty + harvest_df <- data.frame(julfauche = NA, hautcoupe = NA, lairesiduel = NA, msresiduel = NA, anitcoupe = NA) - #add the last sublist back - attr_sublist["nb_interventions"] <- list_no - 2 - tec_list[[8]][[1]][[1]][[2]][[2]][[1]][[".attrs"]] <- attr_sublist - + + # If given harvest date is within simulation days + # probably need to break down >2 years into multiple usms + if(as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01")) %in% dseq){ + + # STICS needs cutting days in cumulative julian days + # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) + # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 + # the following code should give that + harvest_df$julfauche <- which(dseq == as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01"))) + lubridate::yday(dseq[1]) - 2 + harvest_df$lairesiduel <- h_days[hrow, 3] + # also pass cutting height + } + + colnames(harvest_df) <- paste0(h_param_names, "_", hrow) + harvest_list[[hrow]] <- harvest_df } + harvest_tec <- do.call("cbind", harvest_list) - + harvest_tec$codefauche <- 2 # use calendar days } + # TODO: if no harvest file is given modify harvest parameters, e.g. harvest decision - # OTHERWISE DO NOTHING FOR NOW + # DO NOTHING ELSE FOR NOW + # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION - # write the tec file - XML::saveXML(PEcAn.settings::listToXml(tec_list, "fichiertec"), - file = file.path(rundir, paste0(defaults$pft$name, "_tec.xml")), - prefix = '\n') + # same usm -> continue columns + tec_df <- cbind(tec_df, harvest_tec) + SticsRFiles::gen_tec_xml(out_path = rundir, param_table = tec_df) + # TODO: more than 1 USM, rbind + SticsRFiles::convert_xml2txt(xml_file = file.path(rundir, paste0(defaults$pft$name, "_tec.xml")), + java_dir = javastics_path) ################################ Prepare USM file ###################################### From 02e7b1fdd3e1649b02f4476184f66a1e12bb833a Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 28 Jan 2022 15:17:08 +0200 Subject: [PATCH 10/68] save status --- models/stics/R/write.config.STICS.R | 32 ++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index b0405ce4043..ef26c5544c0 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -288,15 +288,18 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_file <- file.path(rundir, "pecan_sta.xml") + sta_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) XML::saveXML(sta_xml, file = sta_file) - SticsRFiles::convert_xml2txt(xml_file = sta_file, java_dir = javastics_path) - sta_txt <- file.path(rundir, "station.txt") + # change latitude + SticsRFiles::set_param_xml(sta_file, "latitude", settings$run$site$lat, overwrite = TRUE) + + SticsRFiles::convert_xml2txt(xml_file = sta_file, java_dir = javastics_path) - # change latitute - SticsRFiles::set_station_txt(sta_txt, param = "latitude", value = settings$run$site$lat) + # another way to change latitute + # sta_txt <- file.path(rundir, "station.txt") + # SticsRFiles::set_station_txt(sta_txt, param = "latitude", value = settings$run$site$lat) # DO NOTHING ELSE FOR NOW # Should these be prepared by met2model.STICS? @@ -387,9 +390,24 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) - usm_list <- XML::xmlToList(usm_xml) + usm_file <- file.path(rundir, "usms.xml") + XML::saveXML(usm_xml, file = usm_file) + + # This may also be easier to generate from a data frame instead of overwriting a template + usms_param_df <- data.frame(usm_name = defaults$pft$name, # pft name + datedebut = lubridate::yday(settings$run$start.date), # beginning day of the simulation (julian.d) + datefin = usm_list$usm$datedebut + length(dseq) - 1, # end day of the simulation (julian.d) (at the end of consecutive years, i.e. can be greater than 366) + finit = paste0(defaults$pft$name, "_ini.xml"), # name of the initialization file + nomsol = paste0("sol", defaults$pft$name), # name of the soil in the sols.xml file + fstation = paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml"), # name of the weather station file + ) - # TODO: more than 1 USM and PFTs (STICS can run 2 PFTs max: main crop + intercrop) + + + # TODO: more than 1 USM and PFTs + # STICS can run 2 PFTs max: main crop + intercrop + + # pft name usm_list$usm$.attrs[["nom"]] <- defaults$pft$name From 3979d195277714c5208c56bc8d996c2791178d6c Mon Sep 17 00:00:00 2001 From: istfer Date: Thu, 24 Feb 2022 12:42:00 +0200 Subject: [PATCH 11/68] fixes for tec file --- models/stics/R/write.config.STICS.R | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index ef26c5544c0..4fbb44b373f 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -288,7 +288,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) + sta_file <- file.path(rundir, settings$pfts$pft$name, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) XML::saveXML(sta_xml, file = sta_file) @@ -318,8 +318,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases - SticsRFiles::gen_tec_xml(out_path = rundir, param_table = tec_df) - + # the following formalisms exist in the tec file: ## supply of organic residus ## soil tillage @@ -336,6 +335,9 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { h_days <- as.matrix(utils::read.table(settings$run$inputs$harvest$path, header = TRUE, sep = ",")) + #temporary dev hack, give 2-years of intervention until multiple usms + h_days <- h_days[1:4,] + # param names h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d "hautcoupe" , # cut height for forage crops, m @@ -378,11 +380,24 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # same usm -> continue columns tec_df <- cbind(tec_df, harvest_tec) - SticsRFiles::gen_tec_xml(out_path = rundir, param_table = tec_df) + + # these shouldn't be empty even if we don't use them (values from timothy example in STICS) + tec_df$iplt0 <- 999 # date of sowing + tec_df$profsem <- 2 # depth of sowing + tec_df$densitesem <- 100 # plant sowing density + tec_df$variete <- 1 # cultivar number corresponding to the cultivar name in the plant file (could be passed via a field activity file) + tec_df$irecbutoir <- 999 #latest date of harvest (imposed if the crop cycle is not finished at this date) + tec_df$profmes <- 120 # depth of measurement of the soil water reserve (cm) + tec_df$engrais <- 1 # fertilizer type + tec_df$concirr <- 0.11 # concentration of mineral N in irrigation water (kg ha-1 mm-1) + tec_df$ressuite <- 'straw+roots' # type of crop residue + tec_df$h2ograinmax <- 0.32 # maximal water content of fruits at harvest + + SticsRFiles::gen_tec_xml(out_path = file.path(rundir, defaults$pft$name), param_table = tec_df) # TODO: more than 1 USM, rbind - SticsRFiles::convert_xml2txt(xml_file = file.path(rundir, paste0(defaults$pft$name, "_tec.xml")), + SticsRFiles::convert_xml2txt(xml_file = file.path(rundir, defaults$pft$name, paste0(defaults$pft$name, "_tec.xml")), java_dir = javastics_path) @@ -452,7 +467,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { usm_list[[1]][[11]]$ftec <- paste0(defaults$pft$name, "_tec.xml") # name of the LAI forcing file for main plant (null if none) - usm_list[[1]][[11]]$flai <- "null" # hardcode for now + usm_list[[1]][[11]]$flai <- "default.lai" # hardcode for now, doesn't matter when codesimul==0 # name of the plant file for associated plant (intercropping) usm_list[[1]][[12]]$fplt <- "null" # hardcode for now From 6e698e6f67ec797325f405e11a3562256f98ccf2 Mon Sep 17 00:00:00 2001 From: istfer Date: Thu, 24 Feb 2022 16:09:52 +0200 Subject: [PATCH 12/68] use new templates --- models/stics/R/write.config.STICS.R | 123 ++++++++++++++-------------- models/stics/inst/template.job | 9 +- models/stics/inst/template.usm | 36 ++++++++ 3 files changed, 103 insertions(+), 65 deletions(-) create mode 100644 models/stics/inst/template.usm diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 4fbb44b373f..e46934c22f0 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -265,8 +265,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol #### TODO: revise others to have txt templates directly in the package - # read in template sols file (xml) - # sols_xml <- XML::xmlParse(system.file("sols.xml", package = "PEcAn.STICS")) sols_file <- file.path(rundir, "param.sol") # cp template sols file (txt) @@ -288,7 +286,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_file <- file.path(rundir, settings$pfts$pft$name, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) + sta_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) XML::saveXML(sta_xml, file = sta_file) @@ -315,6 +313,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## this is where we modify management practices + # probably need to do this before in the code + usmdir <- file.path(rundir, defaults$pft$name) + dir.create(usmdir, showWarnings = FALSE, recursive = TRUE) + ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases @@ -373,7 +375,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { harvest_tec$codefauche <- 2 # use calendar days } - # TODO: if no harvest file is given modify harvest parameters, e.g. harvest decision + # TODO: if no harvest file is given modify other harvest parameters, e.g. harvest decision # DO NOTHING ELSE FOR NOW # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION @@ -393,106 +395,100 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { tec_df$ressuite <- 'straw+roots' # type of crop residue tec_df$h2ograinmax <- 0.32 # maximal water content of fruits at harvest - SticsRFiles::gen_tec_xml(out_path = file.path(rundir, defaults$pft$name), param_table = tec_df) + SticsRFiles::gen_tec_xml(out_path = usmdir, param_table = tec_df) # TODO: more than 1 USM, rbind - SticsRFiles::convert_xml2txt(xml_file = file.path(rundir, defaults$pft$name, paste0(defaults$pft$name, "_tec.xml")), + SticsRFiles::convert_xml2txt(xml_file = file.path(usmdir, paste0(defaults$pft$name, "_tec.xml")), java_dir = javastics_path) ################################ Prepare USM file ###################################### # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. - usm_xml <- XML::xmlParse(system.file("usms.xml", package = "PEcAn.STICS")) - usm_file <- file.path(rundir, "usms.xml") - XML::saveXML(usm_xml, file = usm_file) - - # This may also be easier to generate from a data frame instead of overwriting a template - usms_param_df <- data.frame(usm_name = defaults$pft$name, # pft name - datedebut = lubridate::yday(settings$run$start.date), # beginning day of the simulation (julian.d) - datefin = usm_list$usm$datedebut + length(dseq) - 1, # end day of the simulation (julian.d) (at the end of consecutive years, i.e. can be greater than 366) - finit = paste0(defaults$pft$name, "_ini.xml"), # name of the initialization file - nomsol = paste0("sol", defaults$pft$name), # name of the soil in the sols.xml file - fstation = paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml"), # name of the weather station file - ) + usm_file <- file.path(usmdir, "new_travail.usm") - + # cp template usm file + file.copy(system.file("template.usm", package = "PEcAn.STICS"), usm_file) - # TODO: more than 1 USM and PFTs - # STICS can run 2 PFTs max: main crop + intercrop + # Type of LAI simulation + # 0 = culture (LAI calculated by the model), 1 = feuille (LAI forced) + SticsRFiles::set_usm_txt(usm_file, "codesimul", 0, add = FALSE) # hardcode for now - + # use optimization + # 0 = no; 1 = yes main plant; 2 = yes associated plant + SticsRFiles::set_usm_txt(usm_file, "codeoptim", 0, add = FALSE) - # pft name - usm_list$usm$.attrs[["nom"]] <- defaults$pft$name + # option to simulate several + # successive USM (0 = no, 1 = yes) + SticsRFiles::set_usm_txt(usm_file, "codesuite", 0, add = FALSE) # for now + # number of simulated plants (sole crop=1; intercropping=2) + SticsRFiles::set_usm_txt(usm_file, "nbplantes", 1, add = FALSE) # hardcode for now + + # pft name + SticsRFiles::set_usm_txt(usm_file, "nom", defaults$pft$name, add = FALSE) + # beginning day of the simulation (julian.d) - usm_list$usm$datedebut <- lubridate::yday(settings$run$start.date) + SticsRFiles::set_usm_txt(usm_file, "datedebut", lubridate::yday(settings$run$start.date), add = FALSE) # end day of the simulation (julian.d) (at the end of consecutive years, i.e. can be greater than 366) - usm_list$usm$datefin <- usm_list$usm$datedebut + length(dseq) - 1 + SticsRFiles::set_usm_txt(usm_file, "datefin", (lubridate::yday(settings$run$start.date) + length(dseq) - 1), add = FALSE) # name of the initialization file - usm_list$usm$finit <- paste0(defaults$pft$name, "_ini.xml") + SticsRFiles::set_usm_txt(usm_file, "finit", paste0(defaults$pft$name, "_ini.xml"), add = FALSE) + + # soil number + SticsRFiles::set_usm_txt(usm_file, "numsol", 1, add = FALSE) # name of the soil in the sols.xml file - usm_list$usm$nomsol <- paste0("sol", defaults$pft$name) + SticsRFiles::set_usm_txt(usm_file, "nomsol", paste0("sol", defaults$pft$name), add = FALSE) # name of the weather station file - usm_list$usm$fstation <- paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml") + SticsRFiles::set_usm_txt(usm_file, "fstation", paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml"), add = FALSE) # name of the first climate file - usm_list$usm$fclim1 <- paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$start.date)) - + SticsRFiles::set_usm_txt(usm_file, "fclim1", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$start.date)), add = FALSE) # name of the last climate file - usm_list$usm$fclim2 <- paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$end.date)) + SticsRFiles::set_usm_txt(usm_file, "fclim2", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$end.date)), add = FALSE) + + # number of simulation years + SticsRFiles::set_usm_txt(usm_file, "nbans", 2, add = FALSE) # hardcode for now # number of calendar years involved in the crop cycle # 1 = 1 year e.g. for spring crops, 0 = two years, e.g. for winter crops - usm_list$usm$culturean <- trait.values$timothy$crop_cycle - - # number of simulated plants (sole crop=1; intercropping=2) - usm_list$usm$nbplantes <- 1 # hardcode for now - - # Type of LAI simulation - # 0 = culture (LAI calculated by the model), 1 = feuille (LAI forced) - usm_list$usm$codesimul <- 0 # hardcode for now + SticsRFiles::set_usm_txt(usm_file, "culturean", trait.values$timothy$crop_cycle, add = FALSE) # name of the plant file for main plant - usm_list[[1]][[11]]$fplt <- paste0(defaults$pft$name, "_plt.xml") + SticsRFiles::set_usm_txt(usm_file, "fplt1", paste0(defaults$pft$name, "_plt.xml"), add = FALSE) # name of the technical file for main plant - usm_list[[1]][[11]]$ftec <- paste0(defaults$pft$name, "_tec.xml") + SticsRFiles::set_usm_txt(usm_file, "ftec1", paste0(defaults$pft$name, "_tec.xml"), add = FALSE) # name of the LAI forcing file for main plant (null if none) - usm_list[[1]][[11]]$flai <- "default.lai" # hardcode for now, doesn't matter when codesimul==0 - - # name of the plant file for associated plant (intercropping) - usm_list[[1]][[12]]$fplt <- "null" # hardcode for now + SticsRFiles::set_usm_txt(usm_file, "flai1", "default.lai", add = FALSE) # hardcode for now, doesn't matter when codesimul==0 - # name of the technical file for associated plant (intercropping) - usm_list[[1]][[12]]$ftec <- "null" # hardcode for now - - # name of the LAI forcing file for associated plant (intercropping) (null if none) - usm_list[[1]][[12]]$flai <- "null" # hardcode for now + - # write USMs - XML::saveXML(PEcAn.settings::listToXml(usm_list, "usms"), - file = file.path(rundir, "usms.xml"), - prefix = '\n') + # TODO: more than 1 USM (can use add=TRUE) and PFTs + # STICS can run 2 PFTs max: main crop + intercrop + ################################ Prepare Run ###################################### # symlink climate files met_path <- settings$run$inputs$met$path + clim_list <- list() # temporary solution + ctr <- 1 for(clim in seq(lubridate::year(settings$run$start.date), lubridate::year(settings$run$end.date))){ met_file <- gsub(paste0(lubridate::year(settings$run$start.date), ".climate"), paste0(clim, ".climate"), met_path) - clim_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", clim)) - file.symlink(met_file, clim_file) + clim_list[[ctr]] <- read.table(met_file) + ctr <- ctr + 1 } + clim_run <- do.call("rbind", clim_list) + write.table(clim_run, file.path(usmdir, "climat.txt"), col.names = FALSE, row.names = FALSE) # stics path # stics_path <- settings$model$binary @@ -501,12 +497,12 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { file.symlink(stics_path, bindir) # generate STICS input files using JavaStics - jexe <- file.path(gsub("bin","", dirname(stics_path)), "JavaSticsCmd.exe") + # jexe <- file.path(gsub("bin","", dirname(stics_path)), "JavaSticsCmd.exe") usm_name <- defaults$pft$name # if this script can already create the txts, bypass this step - cmd_generate <- paste("java -jar", jexe,"--generate-txt", rundir, usm_name) + # cmd_generate <- paste("java -jar", jexe,"--generate-txt", rundir, usm_name) # copy *.mod files mod_files <- c(file.path(gsub("bin","example", dirname(stics_path)), "var.mod"), @@ -514,9 +510,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { file.path(gsub("bin","example", dirname(stics_path)), "prof.mod")) file.copy(mod_files, rundir) - cmd_run <- paste("java -jar", jexe,"--run", rundir, usm_name) + #cmd_run <- paste("java -jar", jexe,"--run", rundir, usm_name) + + # using SticsOnR wrapper in job.sh now - SticsOnR::stics_wrapper(model_options = wrapper_options) - #----------------------------------------------------------------------- # create launch script (which will create symlink) @@ -558,8 +555,8 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { jobsh <- gsub("@MODFILE@", paste0("mod_s", usm_name, ".sti"), jobsh) - jobsh <- gsub("@CMD_GENERATE@", cmd_generate, jobsh) - jobsh <- gsub("@CMD_RUN@", cmd_run, jobsh) + #jobsh <- gsub("@CMD_GENERATE@", cmd_generate, jobsh) + jobsh <- gsub("@JAVASTICS_PATH@", javastics_path, jobsh) writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh")) Sys.chmod(file.path(settings$rundir, run.id, "job.sh")) diff --git a/models/stics/inst/template.job b/models/stics/inst/template.job index bf96e69685c..a625bd77b6a 100644 --- a/models/stics/inst/template.job +++ b/models/stics/inst/template.job @@ -15,8 +15,13 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then cd "@RUNDIR@" - @CMD_GENERATE@ - @CMD_RUN@ + # call stics_wrapper + echo " + javastics_path = '@JAVASTICS_PATH@' + sticsrun_dir = '@RUNDIR@' + wrapper_options = SticsOnR::stics_wrapper_options(javastics_path = javastics_path, data_dir = sticsrun_dir) + SticsOnR::stics_wrapper(model_options = wrapper_options) + " | R --vanilla # copy log mv "@RUNDIR@/stics_errors.log" "@OUTDIR@" diff --git a/models/stics/inst/template.usm b/models/stics/inst/template.usm new file mode 100644 index 00000000000..a2aa436effe --- /dev/null +++ b/models/stics/inst/template.usm @@ -0,0 +1,36 @@ +:codesimul +codesimul_placeholder +:codeoptim +codeoptim_placeholder +:codesuite +codesuite_placeholder +:nbplantes +nbplantes_placeholder +:nom +nom_placeholder +:datedebut +datedebut_placeholder +:datefin +datefin_placeholder +:finit +finit_placeholder +:numsol +numsol_placeholder +:nomsol +nomsol_placeholder +:fstation +fstation_placeholder +:fclim1 +fclim1_placeholder +:fclim2 +fclim2_placeholder +:nbans +nbans_placeholder +:culturean +culturean_placeholder +:fplt1 +fplt1_placeholder +:ftec1 +ftec1_placeholder +:flai1 +flai1_placeholder From fd4c8d3718c573b7f4f1f8c7acc29e6a67a92b4b Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 25 Feb 2022 11:53:00 +0200 Subject: [PATCH 13/68] move mod files to the package --- models/stics/inst/prof.mod | 4 ++++ models/stics/inst/rap.mod | 15 +++++++++++++ models/stics/inst/var.mod | 45 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 64 insertions(+) create mode 100644 models/stics/inst/prof.mod create mode 100644 models/stics/inst/rap.mod create mode 100644 models/stics/inst/var.mod diff --git a/models/stics/inst/prof.mod b/models/stics/inst/prof.mod new file mode 100644 index 00000000000..00510174f7d --- /dev/null +++ b/models/stics/inst/prof.mod @@ -0,0 +1,4 @@ +2 +tsol(iz) +10 +01 01 2000 diff --git a/models/stics/inst/rap.mod b/models/stics/inst/rap.mod new file mode 100644 index 00000000000..d64396c2bac --- /dev/null +++ b/models/stics/inst/rap.mod @@ -0,0 +1,15 @@ +1 +1 +2 +1 +rec +masec(n) +mafruit +iflos +imats +irecs +laimax +cprecip +cet +QNplante +Qles diff --git a/models/stics/inst/var.mod b/models/stics/inst/var.mod new file mode 100644 index 00000000000..34bce34862b --- /dev/null +++ b/models/stics/inst/var.mod @@ -0,0 +1,45 @@ +lai(n) +masec(n) +mafruit +HR(1) +HR(2) +HR(3) +HR(4) +HR(5) +resmes +drain +esol +et +zrac +tcult +AZnit(1) +AZnit(2) +AZnit(3) +AZnit(4) +AZnit(5) +Qles +QNplante +azomes +inn +chargefruit +AZamm(1) +AZamm(2) +AZamm(3) +AZamm(4) +AZamm(5) +CNgrain +concNO3les +drat +fapar +hauteur +Hmax +humidite +LRACH(1) +LRACH(2) +LRACH(3) +LRACH(4) +LRACH(5) +mafrais +pdsfruitfrais +Qdrain +rnet From f7278787585b897bd8713d995086a9aa1f683545 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 25 Feb 2022 12:13:30 +0200 Subject: [PATCH 14/68] set usm dir --- models/stics/R/write.config.STICS.R | 44 +++++++++++------------------ models/stics/inst/template.job | 5 ++-- 2 files changed, 19 insertions(+), 30 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index e46934c22f0..a85174c01a1 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -31,7 +31,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # find out where to write run/ouput rundir <- file.path(settings$host$rundir, run.id) - pltdir <- file.path(settings$host$rundir, run.id, "plant") + usmdir <- file.path(settings$host$rundir, run.id, settings$pfts$pft$name) # will think about multiple usm/PFT case cfgdir <- file.path(settings$host$rundir, run.id, "config") bindir <- file.path(settings$host$rundir, run.id, "bin") outdir <- file.path(settings$host$outdir, run.id) @@ -41,7 +41,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { dir.create(outdir, showWarnings = FALSE, recursive = TRUE) ## create plant, config and bin dirs - dir.create(pltdir, showWarnings = FALSE, recursive = TRUE) + dir.create(usmdir, showWarnings = FALSE, recursive = TRUE) dir.create(cfgdir, showWarnings = FALSE, recursive = TRUE) dir.create(bindir, showWarnings = FALSE, recursive = TRUE) @@ -79,7 +79,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { pft.traits <- unlist(trait.values[[pft]]) pft.names <- names(pft.traits) - plant_file <- file.path(pltdir, paste0(names(trait.values)[pft], "_plt.xml")) + plant_file <- file.path(usmdir, paste0(names(trait.values)[pft], "_plt.xml")) if(names(trait.values)[pft] != "env"){ # save the template, will be overwritten below @@ -227,13 +227,13 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## DO NOTHING FOR NOW gen_xml <- XML::xmlParse(system.file("param_gen.xml", package = "PEcAn.STICS")) - gen_file <- file.path(rundir, "param_gen.xml") + gen_file <- file.path(usmdir, "param_gen.xml") XML::saveXML(gen_xml, file = gen_file) SticsRFiles::convert_xml2txt(xml_file = gen_file, java_dir = javastics_path) # may delete the xml after this newf_xml <- XML::xmlParse(system.file("param_newform.xml", package = "PEcAn.STICS")) - newf_file <- file.path(rundir, "param_newform.xml") + newf_file <- file.path(usmdir, "param_newform.xml") XML::saveXML(newf_xml, file = newf_file) SticsRFiles::convert_xml2txt(xml_file = newf_file, java_dir = javastics_path) @@ -245,7 +245,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template ini file ini_xml <- XML::xmlParse(system.file("pecan_ini.xml", package = "PEcAn.STICS")) - ini_file <- file.path(rundir, paste0(defaults$pft$name, "_ini.xml")) + ini_file <- file.path(usmdir, paste0(defaults$pft$name, "_ini.xml")) # write the ini file XML::saveXML(ini_xml, file = ini_file) @@ -265,7 +265,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol #### TODO: revise others to have txt templates directly in the package - sols_file <- file.path(rundir, "param.sol") + sols_file <- file.path(usmdir, "param.sol") # cp template sols file (txt) file.copy(system.file("param.sol", package = "PEcAn.STICS"), sols_file) @@ -286,7 +286,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) + sta_file <- file.path(usmdir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) XML::saveXML(sta_xml, file = sta_file) @@ -313,10 +313,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## this is where we modify management practices - # probably need to do this before in the code - usmdir <- file.path(rundir, defaults$pft$name) - dir.create(usmdir, showWarnings = FALSE, recursive = TRUE) - ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases @@ -490,29 +486,23 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { clim_run <- do.call("rbind", clim_list) write.table(clim_run, file.path(usmdir, "climat.txt"), col.names = FALSE, row.names = FALSE) - # stics path - # stics_path <- settings$model$binary - # symlink to binary file.symlink(stics_path, bindir) - - # generate STICS input files using JavaStics - # jexe <- file.path(gsub("bin","", dirname(stics_path)), "JavaSticsCmd.exe") + stics_exe <- file.path(bindir, basename(stics_path)) usm_name <- defaults$pft$name - # if this script can already create the txts, bypass this step - # cmd_generate <- paste("java -jar", jexe,"--generate-txt", rundir, usm_name) - # copy *.mod files - mod_files <- c(file.path(gsub("bin","example", dirname(stics_path)), "var.mod"), - file.path(gsub("bin","example", dirname(stics_path)), "rap.mod"), - file.path(gsub("bin","example", dirname(stics_path)), "prof.mod")) - file.copy(mod_files, rundir) + file.copy(system.file("var.mod", package = "PEcAn.STICS"), usmdir) + file.copy(system.file("rap.mod", package = "PEcAn.STICS"), usmdir) + file.copy(system.file("prof.mod", package = "PEcAn.STICS"), usmdir) #cmd_run <- paste("java -jar", jexe,"--run", rundir, usm_name) # using SticsOnR wrapper in job.sh now - SticsOnR::stics_wrapper(model_options = wrapper_options) + # used to be: + # cmd_generate <- paste("java -jar", jexe,"--generate-txt", rundir, usm_name) + # cmd_run <- paste("java -jar", jexe,"--run", rundir, usm_name) #----------------------------------------------------------------------- @@ -554,9 +544,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { jobsh <- gsub("@RUNDIR@", rundir, jobsh) jobsh <- gsub("@MODFILE@", paste0("mod_s", usm_name, ".sti"), jobsh) - - #jobsh <- gsub("@CMD_GENERATE@", cmd_generate, jobsh) - jobsh <- gsub("@JAVASTICS_PATH@", javastics_path, jobsh) + jobsh <- gsub("@STICSEXE@", stics_exe, jobsh) writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh")) Sys.chmod(file.path(settings$rundir, run.id, "job.sh")) diff --git a/models/stics/inst/template.job b/models/stics/inst/template.job index a625bd77b6a..403b5d0e27e 100644 --- a/models/stics/inst/template.job +++ b/models/stics/inst/template.job @@ -17,9 +17,10 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then # call stics_wrapper echo " - javastics_path = '@JAVASTICS_PATH@' + javastics_path = '@RUNDIR@' + stics_exe = '@STICSEXE@' sticsrun_dir = '@RUNDIR@' - wrapper_options = SticsOnR::stics_wrapper_options(javastics_path = javastics_path, data_dir = sticsrun_dir) + wrapper_options = SticsOnR::stics_wrapper_options(javastics_path = javastics_path, stics_exe = stics_exe, data_dir = sticsrun_dir, force = TRUE) SticsOnR::stics_wrapper(model_options = wrapper_options) " | R --vanilla From 5e17d7998450256b27376baa987a883235d267bd Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 25 Feb 2022 12:20:26 +0200 Subject: [PATCH 15/68] fix error log dir --- models/stics/R/write.config.STICS.R | 1 + models/stics/inst/template.job | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index a85174c01a1..a1c12266f75 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -542,6 +542,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { jobsh <- gsub("@OUTDIR@", outdir, jobsh) jobsh <- gsub("@RUNDIR@", rundir, jobsh) + jobsh <- gsub("@USMDIR@", usmdir, jobsh) jobsh <- gsub("@MODFILE@", paste0("mod_s", usm_name, ".sti"), jobsh) jobsh <- gsub("@STICSEXE@", stics_exe, jobsh) diff --git a/models/stics/inst/template.job b/models/stics/inst/template.job index 403b5d0e27e..d10a1c73a96 100644 --- a/models/stics/inst/template.job +++ b/models/stics/inst/template.job @@ -25,7 +25,7 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then " | R --vanilla # copy log - mv "@RUNDIR@/stics_errors.log" "@OUTDIR@" + mv "@USMDIR@/stics_errors.log" "@OUTDIR@" STATUS=$? From 04c3893c8b3305f62b530deeded63e12f4adfb37 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 25 Feb 2022 13:01:21 +0200 Subject: [PATCH 16/68] minor revision to model2netcdf --- models/stics/R/model2netcdf.STICS.R | 27 ++++++++++++++++++++++++--- models/stics/R/write.config.STICS.R | 4 +--- models/stics/inst/template.job | 2 +- models/stics/inst/var.mod | 3 +++ 4 files changed, 29 insertions(+), 7 deletions(-) diff --git a/models/stics/R/model2netcdf.STICS.R b/models/stics/R/model2netcdf.STICS.R index 012db5e80bf..92a6419c3ce 100644 --- a/models/stics/R/model2netcdf.STICS.R +++ b/models/stics/R/model2netcdf.STICS.R @@ -22,9 +22,9 @@ ##' @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) @@ -38,7 +38,7 @@ 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 ") + PEcAn.logger::logger.info("Years selected for model run and STICS output years do not match.") } # determine time step? @@ -54,6 +54,26 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o outlist <- list() 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 <- stics_output[thisyear, "CO2sol"] + + #outlist[[2]] <- udunits2::ud.convert(HeteroResp, "ha-1 day-1", "m-2 s-1") + + # need to study NEE in STICS, growth rate is calcultaed after photosynthesis and remobilized reserves + # also check the final utput is not shoot growth only + + # dltams(n): daily growth rate of the plant (t.ha-1.d-1) + dltams <- stics_output[thisyear, "dltams.n."] * 1000 * 0.48 + + # dltaremobil: daily amount of perennial reserves remobilised (t.ha-1.d-1) + dltaremobil <- stics_output[thisyear, "dltaremobil"] * 1000 * 0.48 + + # is it? is autotrophic respiration accounted for? + NPP <- dltams - dltaremobil # kgC ha-1 d-1 + #outlist[[3]] <- udunits2::ud.convert(NPP, "ha-1 day-1", "m-2 s-1") + + outlist[[2]] <- -1 * udunits2::ud.convert(NPP - HeteroResp, "ha-1 day-1", "m-2 s-1") + # ******************** Declare netCDF dimensions and variables ********************# t <- ncdf4::ncdim_def(name = "time", units = paste0("days since ", y, "-01-01 00:00:00"), @@ -69,6 +89,7 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o nc_var <- list() nc_var[[1]] <- PEcAn.utils::to_ncvar("LAI", dims) + nc_var[[2]] <- PEcAn.utils::to_ncvar("NEE", dims) # ******************** Declare netCDF variables ********************# diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index a1c12266f75..aa32d9cda94 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -312,6 +312,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ############################ Prepare Technical File ################################## ## this is where we modify management practices + ## TODO: use ICASA compatible json file ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases @@ -333,9 +334,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { h_days <- as.matrix(utils::read.table(settings$run$inputs$harvest$path, header = TRUE, sep = ",")) - #temporary dev hack, give 2-years of intervention until multiple usms - h_days <- h_days[1:4,] - # param names h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d "hautcoupe" , # cut height for forage crops, m diff --git a/models/stics/inst/template.job b/models/stics/inst/template.job index d10a1c73a96..0320e886f76 100644 --- a/models/stics/inst/template.job +++ b/models/stics/inst/template.job @@ -36,7 +36,7 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then fi # copy output - mv @RUNDIR@/*.sti @OUTDIR@ + mv @USMDIR@/*.sti @OUTDIR@ # convert to MsTMIP echo "library (PEcAn.STICS) diff --git a/models/stics/inst/var.mod b/models/stics/inst/var.mod index 34bce34862b..5e8688406c2 100644 --- a/models/stics/inst/var.mod +++ b/models/stics/inst/var.mod @@ -43,3 +43,6 @@ mafrais pdsfruitfrais Qdrain rnet +CO2sol +dltams(n) +dltaremobil From 7becb2aebceb808d2beae99ac9b609d96eb5fca4 Mon Sep 17 00:00:00 2001 From: istfer Date: Wed, 2 Mar 2022 16:29:20 +0200 Subject: [PATCH 17/68] first pass for multiple usms --- models/stics/R/model2netcdf.STICS.R | 6 +- models/stics/R/write.config.STICS.R | 344 ++++++++++++++++------------ models/stics/inst/template.job | 9 +- 3 files changed, 210 insertions(+), 149 deletions(-) diff --git a/models/stics/R/model2netcdf.STICS.R b/models/stics/R/model2netcdf.STICS.R index 92a6419c3ce..f97ff5663f5 100644 --- a/models/stics/R/model2netcdf.STICS.R +++ b/models/stics/R/model2netcdf.STICS.R @@ -29,7 +29,11 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o 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, read.table, header = T, 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) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index aa32d9cda94..944d9d6ae61 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -31,19 +31,32 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # find out where to write run/ouput rundir <- file.path(settings$host$rundir, run.id) - usmdir <- file.path(settings$host$rundir, run.id, settings$pfts$pft$name) # will think about multiple usm/PFT case cfgdir <- file.path(settings$host$rundir, run.id, "config") bindir <- file.path(settings$host$rundir, run.id, "bin") outdir <- file.path(settings$host$outdir, run.id) + # will think about multiple PFT case, currently considering successive (>2 years) case only + usmdir_root <- paste0(file.path(settings$host$rundir, run.id, settings$pfts$pft$name), "_") + # In STICS, it is 1 UMS per crop cycle, where the cycle can be 2-years max + # If we have a consecutive monoculture for > 2 years, we still need to divide it into 2-year USMs + years_requested <- unique(lubridate::year(dseq)) + if(length(years_requested) %%2 == 1) years_requested <- c(years_requested, years_requested[length(years_requested)]) + + if(length(years_requested) > 2){ + years_indices <- rep(seq(1, length(years_requested), by=2), each=2) + usmdirs <- tapply(years_requested, years_indices, function(x) paste0(usmdir_root, paste(x, collapse = '-'))) + }else{ + usmdirs <- paste0(usmdir_root, paste(years_requested, collapse = '-')) + } + ## make sure rundir and outdir exist dir.create(rundir, showWarnings = FALSE, recursive = TRUE) dir.create(outdir, showWarnings = FALSE, recursive = TRUE) - ## create plant, config and bin dirs - dir.create(usmdir, showWarnings = FALSE, recursive = TRUE) - dir.create(cfgdir, showWarnings = FALSE, recursive = TRUE) - dir.create(bindir, showWarnings = FALSE, recursive = TRUE) + ## create usm, config and bin dirs + dir.create(cfgdir, showWarnings = FALSE, recursive = TRUE) + dir.create(bindir, showWarnings = FALSE, recursive = TRUE) + sapply(usmdirs, dir.create, showWarnings = FALSE, recursive = TRUE) # write preferences prf_xml <- XML::xmlParse(system.file("preferences.xml", package = "PEcAn.STICS")) @@ -79,7 +92,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { pft.traits <- unlist(trait.values[[pft]]) pft.names <- names(pft.traits) - plant_file <- file.path(usmdir, paste0(names(trait.values)[pft], "_plt.xml")) + plant_file <- file.path(rundir, paste0(names(trait.values)[pft], "_plt.xml")) if(names(trait.values)[pft] != "env"){ # save the template, will be overwritten below @@ -222,20 +235,23 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { } # pft-loop ends + file.copy(file.path(rundir, "ficplt1.txt"), file.path(usmdirs, "ficplt1.txt")) ############################## Param gen / newform #################################### ## DO NOTHING FOR NOW gen_xml <- XML::xmlParse(system.file("param_gen.xml", package = "PEcAn.STICS")) - gen_file <- file.path(usmdir, "param_gen.xml") + gen_file <- file.path(rundir, "param_gen.xml") XML::saveXML(gen_xml, file = gen_file) SticsRFiles::convert_xml2txt(xml_file = gen_file, java_dir = javastics_path) # may delete the xml after this + file.copy(file.path(rundir, "tempopar.sti"), file.path(usmdirs, "tempopar.sti")) newf_xml <- XML::xmlParse(system.file("param_newform.xml", package = "PEcAn.STICS")) - newf_file <- file.path(usmdir, "param_newform.xml") + newf_file <- file.path(rundir, "param_newform.xml") XML::saveXML(newf_xml, file = newf_file) SticsRFiles::convert_xml2txt(xml_file = newf_file, java_dir = javastics_path) + file.copy(file.path(rundir, "tempoparv6.sti"), file.path(usmdirs, "tempoparv6.sti")) @@ -245,7 +261,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template ini file ini_xml <- XML::xmlParse(system.file("pecan_ini.xml", package = "PEcAn.STICS")) - ini_file <- file.path(usmdir, paste0(defaults$pft$name, "_ini.xml")) + ini_file <- file.path(rundir, paste0(defaults$pft$name, "_ini.xml")) # write the ini file XML::saveXML(ini_xml, file = ini_file) @@ -255,7 +271,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # SticsRFiles::set_param_xml(xml_file = ini_file, param_name = "lai0", param_value = 1, select = "plante", value = "1", overwrite = TRUE) SticsRFiles::convert_xml2txt(xml_file = ini_file, java_dir = javastics_path) - + file.copy(file.path(rundir, "ficini.txt"), file.path(usmdirs, "ficini.txt")) ############################ Prepare Soils ################################## @@ -265,7 +281,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol #### TODO: revise others to have txt templates directly in the package - sols_file <- file.path(usmdir, "param.sol") + sols_file <- file.path(rundir, "param.sol") # cp template sols file (txt) file.copy(system.file("param.sol", package = "PEcAn.STICS"), sols_file) @@ -274,6 +290,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # sols_vals <- SticsRFiles::get_soil_txt(sols_file) SticsRFiles::set_soil_txt(filepath = sols_file, param="typsol", value=paste0("sol", defaults$pft$name)) + file.copy(sols_file, file.path(usmdirs, "param.sol")) # DO NOTHING ELSE FOR NOW @@ -286,7 +303,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # read in template sta file sta_xml <- XML::xmlParse(system.file("pecan_sta.xml", package = "PEcAn.STICS")) - sta_file <- file.path(usmdir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) + sta_file <- file.path(rundir, paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml")) XML::saveXML(sta_xml, file = sta_file) @@ -294,6 +311,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { SticsRFiles::set_param_xml(sta_file, "latitude", settings$run$site$lat, overwrite = TRUE) SticsRFiles::convert_xml2txt(xml_file = sta_file, java_dir = javastics_path) + file.copy(file.path(rundir, "station.txt"), file.path(usmdirs, "station.txt")) # another way to change latitute # sta_txt <- file.path(rundir, "station.txt") @@ -317,6 +335,17 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## instead of using a template, this could be easier if we prepare a dataframe and use SticsRFiles::gen_tec_xml tec_df <- data.frame(Tec_name = paste0(defaults$pft$name, "_tec.xml")) # note more than one PFT cases + # these shouldn't be empty even if we don't use them (values from timothy example in STICS) + tec_df$iplt0 <- 999 # date of sowing + tec_df$profsem <- 2 # depth of sowing + tec_df$densitesem <- 100 # plant sowing density + tec_df$variete <- 1 # cultivar number corresponding to the cultivar name in the plant file (could be passed via a field activity file) + tec_df$irecbutoir <- 999 #latest date of harvest (imposed if the crop cycle is not finished at this date) + tec_df$profmes <- 120 # depth of measurement of the soil water reserve (cm) + tec_df$engrais <- 1 # fertilizer type + tec_df$concirr <- 0.11 # concentration of mineral N in irrigation water (kg ha-1 mm-1) + tec_df$ressuite <- 'straw+roots' # type of crop residue + tec_df$h2ograinmax <- 0.32 # maximal water content of fruits at harvest # the following formalisms exist in the tec file: ## supply of organic residus @@ -333,7 +362,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { if(!is.null(settings$run$inputs$harvest)){ h_days <- as.matrix(utils::read.table(settings$run$inputs$harvest$path, header = TRUE, sep = ",")) - + # param names h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d "hautcoupe" , # cut height for forage crops, m @@ -341,132 +370,145 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { "msresiduel" , # residual aerial biomass after a cut of a forage crop, t.ha-1 "anitcoupe") # amount of mineral N added by fertiliser application at each cut of a forage crop, kg.ha-1 - - harvest_list <- list() - for(hrow in seq_len(nrow(h_days))){ - - # empty - harvest_df <- data.frame(julfauche = NA, hautcoupe = NA, lairesiduel = NA, msresiduel = NA, anitcoupe = NA) + # loop for each USM + for(usmi in seq_along(usmdirs)){ + usm_years <- years_requested[(usmi*2-1):(usmi*2)] + h_sub <- h_days[h_days[,1] %in% usm_years, ] + dseq_sub <- dseq[lubridate::year(dseq) %in% usm_years] - # If given harvest date is within simulation days - # probably need to break down >2 years into multiple usms - if(as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01")) %in% dseq){ + harvest_list <- list() + for(hrow in seq_len(nrow(h_sub))){ + + # empty + harvest_df <- data.frame(julfauche = NA, hautcoupe = NA, lairesiduel = NA, msresiduel = NA, anitcoupe = NA) - # STICS needs cutting days in cumulative julian days - # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) - # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 - # the following code should give that - harvest_df$julfauche <- which(dseq == as.Date(h_days[hrow, 2], origin = paste0(h_days[hrow, 1], "-01-01"))) + lubridate::yday(dseq[1]) - 2 - harvest_df$lairesiduel <- h_days[hrow, 3] - # also pass cutting height + + # If given harvest date is within simulation days + # probably need to break down >2 years into multiple usms + if(as.Date(h_sub[hrow, 2], origin = paste0(h_sub[hrow, 1], "-01-01")) %in% dseq_sub){ + + # STICS needs cutting days in cumulative julian days + # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) + # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 + # the following code should give that + harvest_df$julfauche <- which(dseq_sub == as.Date(h_sub[hrow, 2], origin = paste0(h_sub[hrow, 1], "-01-01"))) + lubridate::yday(dseq_sub[1]) - 2 + harvest_df$lairesiduel <- h_sub[hrow, 3] + # also pass cutting height + } + + colnames(harvest_df) <- paste0(h_param_names, "_", hrow) + harvest_list[[hrow]] <- harvest_df } + harvest_tec <- do.call("cbind", harvest_list) - colnames(harvest_df) <- paste0(h_param_names, "_", hrow) - harvest_list[[hrow]] <- harvest_df - } - harvest_tec <- do.call("cbind", harvest_list) - - harvest_tec$codefauche <- 2 # use calendar days - } - # TODO: if no harvest file is given modify other harvest parameters, e.g. harvest decision - - # DO NOTHING ELSE FOR NOW - # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION - - # same usm -> continue columns - tec_df <- cbind(tec_df, harvest_tec) - - # these shouldn't be empty even if we don't use them (values from timothy example in STICS) - tec_df$iplt0 <- 999 # date of sowing - tec_df$profsem <- 2 # depth of sowing - tec_df$densitesem <- 100 # plant sowing density - tec_df$variete <- 1 # cultivar number corresponding to the cultivar name in the plant file (could be passed via a field activity file) - tec_df$irecbutoir <- 999 #latest date of harvest (imposed if the crop cycle is not finished at this date) - tec_df$profmes <- 120 # depth of measurement of the soil water reserve (cm) - tec_df$engrais <- 1 # fertilizer type - tec_df$concirr <- 0.11 # concentration of mineral N in irrigation water (kg ha-1 mm-1) - tec_df$ressuite <- 'straw+roots' # type of crop residue - tec_df$h2ograinmax <- 0.32 # maximal water content of fruits at harvest - - SticsRFiles::gen_tec_xml(out_path = usmdir, param_table = tec_df) - - # TODO: more than 1 USM, rbind - - SticsRFiles::convert_xml2txt(xml_file = file.path(usmdir, paste0(defaults$pft$name, "_tec.xml")), - java_dir = javastics_path) - + harvest_tec$codefauche <- 2 # use calendar days + + # DO NOTHING ELSE FOR NOW + # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION + + # same usm -> continue columns + usm_tec_df <- cbind(tec_df, harvest_tec) + + SticsRFiles::gen_tec_xml(out_path = usmdirs[usmi], param_table = usm_tec_df) + + # TODO: more than 1 USM, rbind + + SticsRFiles::convert_xml2txt(xml_file = file.path(usmdirs[usmi], paste0(defaults$pft$name, "_tec.xml")), + java_dir = javastics_path) + + } # end-loop over usms + } # TODO: if no harvest file is given modify other harvest parameters, e.g. harvest decision ################################ Prepare USM file ###################################### - # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. - usm_file <- file.path(usmdir, "new_travail.usm") - - # cp template usm file - file.copy(system.file("template.usm", package = "PEcAn.STICS"), usm_file) + # loop for each USM + #ncodesuite <- ifelse(length(usmdirs) > 1, 1,0) - # Type of LAI simulation - # 0 = culture (LAI calculated by the model), 1 = feuille (LAI forced) - SticsRFiles::set_usm_txt(usm_file, "codesimul", 0, add = FALSE) # hardcode for now - - # use optimization - # 0 = no; 1 = yes main plant; 2 = yes associated plant - SticsRFiles::set_usm_txt(usm_file, "codeoptim", 0, add = FALSE) - - # option to simulate several - # successive USM (0 = no, 1 = yes) - SticsRFiles::set_usm_txt(usm_file, "codesuite", 0, add = FALSE) # for now - - # number of simulated plants (sole crop=1; intercropping=2) - SticsRFiles::set_usm_txt(usm_file, "nbplantes", 1, add = FALSE) # hardcode for now - - # pft name - SticsRFiles::set_usm_txt(usm_file, "nom", defaults$pft$name, add = FALSE) + for(usmi in seq_along(usmdirs)){ - # beginning day of the simulation (julian.d) - SticsRFiles::set_usm_txt(usm_file, "datedebut", lubridate::yday(settings$run$start.date), add = FALSE) - - # end day of the simulation (julian.d) (at the end of consecutive years, i.e. can be greater than 366) - SticsRFiles::set_usm_txt(usm_file, "datefin", (lubridate::yday(settings$run$start.date) + length(dseq) - 1), add = FALSE) - - # name of the initialization file - SticsRFiles::set_usm_txt(usm_file, "finit", paste0(defaults$pft$name, "_ini.xml"), add = FALSE) - - # soil number - SticsRFiles::set_usm_txt(usm_file, "numsol", 1, add = FALSE) - - # name of the soil in the sols.xml file - SticsRFiles::set_usm_txt(usm_file, "nomsol", paste0("sol", defaults$pft$name), add = FALSE) - - # name of the weather station file - SticsRFiles::set_usm_txt(usm_file, "fstation", paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml"), add = FALSE) - - # name of the first climate file - SticsRFiles::set_usm_txt(usm_file, "fclim1", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$start.date)), add = FALSE) - - # name of the last climate file - SticsRFiles::set_usm_txt(usm_file, "fclim2", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", lubridate::year(settings$run$end.date)), add = FALSE) - - # number of simulation years - SticsRFiles::set_usm_txt(usm_file, "nbans", 2, add = FALSE) # hardcode for now - - # number of calendar years involved in the crop cycle - # 1 = 1 year e.g. for spring crops, 0 = two years, e.g. for winter crops - SticsRFiles::set_usm_txt(usm_file, "culturean", trait.values$timothy$crop_cycle, add = FALSE) - - # name of the plant file for main plant - SticsRFiles::set_usm_txt(usm_file, "fplt1", paste0(defaults$pft$name, "_plt.xml"), add = FALSE) - - # name of the technical file for main plant - SticsRFiles::set_usm_txt(usm_file, "ftec1", paste0(defaults$pft$name, "_tec.xml"), add = FALSE) - - # name of the LAI forcing file for main plant (null if none) - SticsRFiles::set_usm_txt(usm_file, "flai1", "default.lai", add = FALSE) # hardcode for now, doesn't matter when codesimul==0 - + usm_years <- years_requested[(usmi*2-1):(usmi*2)] + dseq_sub <- dseq[lubridate::year(dseq) %in% usm_years] + + # read in template USM (Unit of SiMulation) file, has the master settings, file names etc. + usm_file <- file.path(usmdirs[usmi], "new_travail.usm") + + # cp template usm file + file.copy(system.file("template.usm", package = "PEcAn.STICS"), usm_file) + + # Type of LAI simulation + # 0 = culture (LAI calculated by the model), 1 = feuille (LAI forced) + SticsRFiles::set_usm_txt(usm_file, "codesimul", 0, add = FALSE) # hardcode for now + + # use optimization + # 0 = no; 1 = yes main plant; 2 = yes associated plant + SticsRFiles::set_usm_txt(usm_file, "codeoptim", 0, add = FALSE) + + # option to simulate several + # successive USM (0 = no, 1 = yes) + if(usmi == 1){ + SticsRFiles::set_usm_txt(usm_file, "codesuite", 0, add = FALSE) + }else{ + SticsRFiles::set_usm_txt(usm_file, "codesuite", 1, add = FALSE) + } + + + # number of simulated plants (sole crop=1; intercropping=2) + SticsRFiles::set_usm_txt(usm_file, "nbplantes", 1, add = FALSE) # hardcode for now + + # pft name + SticsRFiles::set_usm_txt(usm_file, "nom", basename(usmdirs[usmi]), add = FALSE) + + + ## handle dates, also for partial year(s) + if(usmi == 1){ + # beginning day of the simulation (julian.d) + # end day of the simulation (julian.d) (at the end of consecutive years, i.e. can be greater than 366) + SticsRFiles::set_usm_txt(usm_file, "datedebut", lubridate::yday(settings$run$start.date), add = FALSE) + SticsRFiles::set_usm_txt(usm_file, "datefin", (lubridate::yday(settings$run$start.date) + length(dseq_sub) - 1), add = FALSE) + }else{ + SticsRFiles::set_usm_txt(usm_file, "datedebut", 1, add = FALSE) + SticsRFiles::set_usm_txt(usm_file, "datefin", length(dseq_sub), add = FALSE) + } + + # name of the initialization file + SticsRFiles::set_usm_txt(usm_file, "finit", paste0(defaults$pft$name, "_ini.xml"), add = FALSE) + + # soil number + SticsRFiles::set_usm_txt(usm_file, "numsol", 1, add = FALSE) + + # name of the soil in the sols.xml file + SticsRFiles::set_usm_txt(usm_file, "nomsol", paste0("sol", defaults$pft$name), add = FALSE) + + # name of the weather station file + SticsRFiles::set_usm_txt(usm_file, "fstation", paste0(tolower(sub(" .*", "", settings$run$site$name)), "_sta.xml"), add = FALSE) + + # name of the first climate file + SticsRFiles::set_usm_txt(usm_file, "fclim1", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", usm_years[1]), add = FALSE) + + # name of the last climate file + SticsRFiles::set_usm_txt(usm_file, "fclim2", paste0(tolower(sub(" .*", "", settings$run$site$name)), ".", usm_years[2]), add = FALSE) + + # number of simulation years + SticsRFiles::set_usm_txt(usm_file, "nbans", 2, add = FALSE) # hardcode for now + + # number of calendar years involved in the crop cycle + # 1 = 1 year e.g. for spring crops, 0 = two years, e.g. for winter crops + SticsRFiles::set_usm_txt(usm_file, "culturean", trait.values$timothy$crop_cycle, add = FALSE) + + # name of the plant file for main plant + SticsRFiles::set_usm_txt(usm_file, "fplt1", paste0(defaults$pft$name, "_plt.xml"), add = FALSE) + + # name of the technical file for main plant + SticsRFiles::set_usm_txt(usm_file, "ftec1", paste0(defaults$pft$name, "_tec.xml"), add = FALSE) + + # name of the LAI forcing file for main plant (null if none) + SticsRFiles::set_usm_txt(usm_file, "flai1", "default.lai", add = FALSE) # hardcode for now, doesn't matter when codesimul==0 + + # TODO: more than 1 PFTs + # STICS can run 2 PFTs max: main crop + intercrop + } - - # TODO: more than 1 USM (can use add=TRUE) and PFTs - # STICS can run 2 PFTs max: main crop + intercrop @@ -474,26 +516,31 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # symlink climate files met_path <- settings$run$inputs$met$path - clim_list <- list() # temporary solution - ctr <- 1 - for(clim in seq(lubridate::year(settings$run$start.date), lubridate::year(settings$run$end.date))){ - met_file <- gsub(paste0(lubridate::year(settings$run$start.date), ".climate"), paste0(clim, ".climate"), met_path) - clim_list[[ctr]] <- read.table(met_file) - ctr <- ctr + 1 + + for(usmi in seq_along(usmdirs)){ + + usm_years <- unique(years_requested[(usmi*2-1):(usmi*2)]) + dseq_sub <- dseq[lubridate::year(dseq) %in% usm_years] + + clim_list <- list() # temporary solution + for(clim in seq_along(usm_years)){ + # currently assuming only first year file has been passed to the settings, modify met2model if changing the logic + met_file <- gsub(paste0(lubridate::year(settings$run$start.date), ".climate"), paste0(usm_years[clim], ".climate"), met_path) + clim_list[[clim]] <- read.table(met_file) + } + clim_run <- do.call("rbind", clim_list) + write.table(clim_run, file.path(usmdirs[usmi], "climat.txt"), col.names = FALSE, row.names = FALSE) + } - clim_run <- do.call("rbind", clim_list) - write.table(clim_run, file.path(usmdir, "climat.txt"), col.names = FALSE, row.names = FALSE) - + # symlink to binary file.symlink(stics_path, bindir) stics_exe <- file.path(bindir, basename(stics_path)) - usm_name <- defaults$pft$name - - # copy *.mod files - file.copy(system.file("var.mod", package = "PEcAn.STICS"), usmdir) - file.copy(system.file("rap.mod", package = "PEcAn.STICS"), usmdir) - file.copy(system.file("prof.mod", package = "PEcAn.STICS"), usmdir) + # symlink *.mod files + file.symlink(system.file("var.mod", package = "PEcAn.STICS"), file.path(usmdirs, "var.mod")) + file.symlink(system.file("rap.mod", package = "PEcAn.STICS"), file.path(usmdirs, "rap.mod")) + file.symlink(system.file("prof.mod", package = "PEcAn.STICS"), file.path(usmdirs, "prof.mod")) #cmd_run <- paste("java -jar", jexe,"--run", rundir, usm_name) @@ -540,9 +587,16 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { jobsh <- gsub("@OUTDIR@", outdir, jobsh) jobsh <- gsub("@RUNDIR@", rundir, jobsh) - jobsh <- gsub("@USMDIR@", usmdir, jobsh) - jobsh <- gsub("@MODFILE@", paste0("mod_s", usm_name, ".sti"), jobsh) + if(length(usmdirs)>1){ + jobsh <- gsub("@SUCCESSIVE_USMS@", paste0("list(c('", paste(basename(usmdirs), collapse="','"), "'))"), jobsh) + }else{ + jobsh <- gsub("@SUCCESSIVE_USMS@", 'NULL', jobsh) + } + + jobsh <- gsub("@USMDIR@", usmdirs[1], jobsh) # for now + + jobsh <- gsub("@MODFILE@", paste0("mod_s", basename(usmdirs[1]), ".sti"), jobsh) jobsh <- gsub("@STICSEXE@", stics_exe, jobsh) writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh")) diff --git a/models/stics/inst/template.job b/models/stics/inst/template.job index 0320e886f76..02d0312e995 100644 --- a/models/stics/inst/template.job +++ b/models/stics/inst/template.job @@ -20,12 +20,13 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then javastics_path = '@RUNDIR@' stics_exe = '@STICSEXE@' sticsrun_dir = '@RUNDIR@' - wrapper_options = SticsOnR::stics_wrapper_options(javastics_path = javastics_path, stics_exe = stics_exe, data_dir = sticsrun_dir, force = TRUE) + successive_usms = @SUCCESSIVE_USMS@ + wrapper_options = SticsOnR::stics_wrapper_options(javastics_path = javastics_path, stics_exe = stics_exe, data_dir = sticsrun_dir, force = TRUE, successive_usms = successive_usms) SticsOnR::stics_wrapper(model_options = wrapper_options) " | R --vanilla # copy log - mv "@USMDIR@/stics_errors.log" "@OUTDIR@" + cp "@USMDIR@/stics_errors.log" "@OUTDIR@" STATUS=$? @@ -36,7 +37,9 @@ if [ ! -e "@OUTDIR@/@MODFILE@" ]; then fi # copy output - mv @USMDIR@/*.sti @OUTDIR@ + mv @RUNDIR@/**/mod_b* @OUTDIR@ + mv @RUNDIR@/**/mod_s* @OUTDIR@ + # convert to MsTMIP echo "library (PEcAn.STICS) From dc6c0c1ded765e36a71fceb22c5cd07cb2070e70 Mon Sep 17 00:00:00 2001 From: istfer Date: Thu, 3 Mar 2022 16:48:35 +0200 Subject: [PATCH 18/68] changes to fix met input workflow --- base/db/R/dbfiles.R | 2 +- modules/data.atmosphere/R/met.process.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/base/db/R/dbfiles.R b/base/db/R/dbfiles.R index 28509daef8f..c79501667a4 100644 --- a/base/db/R/dbfiles.R +++ b/base/db/R/dbfiles.R @@ -51,7 +51,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, formatid <- get.id( table = "formats", colnames = c("mimetype_id", "name"), - values = c(mimetypeid, formatname), + values = c(format(mimetypeid), formatname), con = con, create = TRUE, dates = TRUE diff --git a/modules/data.atmosphere/R/met.process.R b/modules/data.atmosphere/R/met.process.R index 3ce070eb42d..538f96fbfed 100644 --- a/modules/data.atmosphere/R/met.process.R +++ b/modules/data.atmosphere/R/met.process.R @@ -275,7 +275,7 @@ met.process <- function(site, input_met, start_date, end_date, model, format.vars = format.vars, bety = con) } else { - if (! met %in% c("ERA5")) cf.id = input_met$id + if (! met %in% c("ERA5", "FieldObservatory")) cf.id = input_met$id } #--------------------------------------------------------------------------------------------------# From 96f5296ed4aa73a4976e72669045488a47663ef2 Mon Sep 17 00:00:00 2001 From: istfer Date: Sat, 5 Mar 2022 12:09:31 +0200 Subject: [PATCH 19/68] switch using json for management --- models/stics/R/write.config.STICS.R | 111 +++++++++++++++------------- 1 file changed, 61 insertions(+), 50 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 944d9d6ae61..0867448b54a 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -358,67 +358,78 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## special techniques ## soil modification by techniques (compaction-fragmentation) - # if a harvest file is given, most (all?) of our harvest cases are actually fall under special techniques - cut crop - if(!is.null(settings$run$inputs$harvest)){ + # if a field activity file is given, most (all?) of our harvest cases are actually fall under special techniques - cut crop + if(!is.null(settings$run$inputs$field.activity)){ - h_days <- as.matrix(utils::read.table(settings$run$inputs$harvest$path, header = TRUE, sep = ",")) - - # param names - h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d - "hautcoupe" , # cut height for forage crops, m - "lairesiduel", # residual LAI after each cut of forage crop, m2 m-2 - "msresiduel" , # residual aerial biomass after a cut of a forage crop, t.ha-1 - "anitcoupe") # amount of mineral N added by fertiliser application at each cut of a forage crop, kg.ha-1 - + events_file <- jsonlite::read_json(settings$run$inputs$field.activity$path, simplifyVector = TRUE)[[1]] # loop for each USM for(usmi in seq_along(usmdirs)){ - usm_years <- years_requested[(usmi*2-1):(usmi*2)] - h_sub <- h_days[h_days[,1] %in% usm_years, ] + usm_years <- years_requested[(usmi*2-1):(usmi*2)] dseq_sub <- dseq[lubridate::year(dseq) %in% usm_years] - harvest_list <- list() - for(hrow in seq_len(nrow(h_sub))){ - - # empty - harvest_df <- data.frame(julfauche = NA, hautcoupe = NA, lairesiduel = NA, msresiduel = NA, anitcoupe = NA) - - - # If given harvest date is within simulation days - # probably need to break down >2 years into multiple usms - if(as.Date(h_sub[hrow, 2], origin = paste0(h_sub[hrow, 1], "-01-01")) %in% dseq_sub){ - - # STICS needs cutting days in cumulative julian days - # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) - # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 - # the following code should give that - harvest_df$julfauche <- which(dseq_sub == as.Date(h_sub[hrow, 2], origin = paste0(h_sub[hrow, 1], "-01-01"))) + lubridate::yday(dseq_sub[1]) - 2 - harvest_df$lairesiduel <- h_sub[hrow, 3] - # also pass cutting height - } + events_sub <- events_file$events[lubridate::year(events_file$events$date) %in% usm_years, ] + + if("planting" %in% events_sub$mgmt_operations_event){ - colnames(harvest_df) <- paste0(h_param_names, "_", hrow) - harvest_list[[hrow]] <- harvest_df + pl_date <- events_sub$date[events_sub$mgmt_operations_event == "planting"] + tec_df$iplt0 <- lubridate::yday(as.Date(pl_date)) } - harvest_tec <- do.call("cbind", harvest_list) - - harvest_tec$codefauche <- 2 # use calendar days - - # DO NOTHING ELSE FOR NOW - # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION - - # same usm -> continue columns - usm_tec_df <- cbind(tec_df, harvest_tec) - - SticsRFiles::gen_tec_xml(out_path = usmdirs[usmi], param_table = usm_tec_df) - - # TODO: more than 1 USM, rbind - SticsRFiles::convert_xml2txt(xml_file = file.path(usmdirs[usmi], paste0(defaults$pft$name, "_tec.xml")), - java_dir = javastics_path) + if("harvest" %in% events_sub$mgmt_operations_event){ + # param names + h_param_names <- c("julfauche" , # date of each cut for forage crops, julian.d + "hautcoupe" , # cut height for forage crops, m + "lairesiduel", # residual LAI after each cut of forage crop, m2 m-2 + "msresiduel" , # residual aerial biomass after a cut of a forage crop, t.ha-1 + "anitcoupe") # amount of mineral N added by fertiliser application at each cut of a forage crop, kg.ha-1 + + harvest_sub <- events_sub[events_sub$mgmt_operations_event == "harvest",] + + harvest_list <- list() + for(hrow in seq_len(nrow(harvest_sub))){ + + # empty + harvest_df <- data.frame(julfauche = NA, hautcoupe = NA, lairesiduel = NA, msresiduel = NA, anitcoupe = NA) + + + # If given harvest date is within simulation days + # probably need to break down >2 years into multiple usms + if(as.Date(harvest_sub$date[hrow]) %in% dseq_sub){ + + # STICS needs cutting days in cumulative julian days + # e.g. first cutting day of the first simulation year can be 163 (2018-06-13) + # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 + # the following code should give that + harvest_df$julfauche <- which(dseq_sub == as.Date(harvest_sub$date[hrow])) + lubridate::yday(dseq_sub[1]) - 1 + harvest_df$lairesiduel <- 0.1 # hardcode for now + # also pass cutting height + } + + colnames(harvest_df) <- paste0(h_param_names, "_", hrow) + harvest_list[[hrow]] <- harvest_df + } + harvest_tec <- do.call("cbind", harvest_list) + + harvest_tec$codefauche <- 2 # use calendar days + + # DO NOTHING ELSE FOR NOW + # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION + + # same usm -> continue columns + usm_tec_df <- cbind(tec_df, harvest_tec) + + SticsRFiles::gen_tec_xml(out_path = usmdirs[usmi], param_table = usm_tec_df) + + # TODO: more than 1 USM, rbind + + SticsRFiles::convert_xml2txt(xml_file = file.path(usmdirs[usmi], paste0(defaults$pft$name, "_tec.xml")), + java_dir = javastics_path) + + } #harvest-if end } # end-loop over usms - } # TODO: if no harvest file is given modify other harvest parameters, e.g. harvest decision + } # TODO: if no events file is given modify other harvest parameters, e.g. harvest decision ################################ Prepare USM file ###################################### From 25cb6f8e8eed0b530f2f15d7d6b2177a2f6f587c Mon Sep 17 00:00:00 2001 From: istfer Date: Sat, 5 Mar 2022 13:20:08 +0200 Subject: [PATCH 20/68] change tag --- models/stics/R/write.config.STICS.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 0867448b54a..089db6d3c5d 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -359,9 +359,9 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ## soil modification by techniques (compaction-fragmentation) # if a field activity file is given, most (all?) of our harvest cases are actually fall under special techniques - cut crop - if(!is.null(settings$run$inputs$field.activity)){ + if(!is.null(settings$run$inputs$fielddata)){ - events_file <- jsonlite::read_json(settings$run$inputs$field.activity$path, simplifyVector = TRUE)[[1]] + events_file <- jsonlite::read_json(settings$run$inputs$fielddata$path, simplifyVector = TRUE)[[1]] # loop for each USM for(usmi in seq_along(usmdirs)){ From 18ee0735e78161516721981ef919ca7e0ff62d26 Mon Sep 17 00:00:00 2001 From: istfer Date: Sun, 6 Mar 2022 11:11:17 +0200 Subject: [PATCH 21/68] better params --- models/stics/R/write.config.STICS.R | 13 ++++++++++--- models/stics/inst/crop_plt.xml | 14 +++++++------- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 089db6d3c5d..8a049dc8341 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -403,8 +403,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # in following years it should be cumulative, meaning a cutting day on 2019-06-12 is 527, not 162 # the following code should give that harvest_df$julfauche <- which(dseq_sub == as.Date(harvest_sub$date[hrow])) + lubridate::yday(dseq_sub[1]) - 1 - harvest_df$lairesiduel <- 0.1 # hardcode for now - # also pass cutting height + harvest_df$hautcoupe <- 0.1 # hardcode for now + harvest_df$lairesiduel <- 0.2 # hardcode for now + harvest_df$msresiduel <- 0 + harvest_df$anitcoupe <- 0 } colnames(harvest_df) <- paste0(h_param_names, "_", hrow) @@ -412,7 +414,11 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { } harvest_tec <- do.call("cbind", harvest_list) - harvest_tec$codefauche <- 2 # use calendar days + harvest_tec$codefauche <- 1 # cut crop - 1:yes, 2:no + harvest_tec$mscoupemini <- 0 # min val of aerial biomass to make a cut + harvest_tec$codemodfauche <- 2 # use calendar days + harvest_tec$hautcoupedefaut <- 0.05 # cut height for forage crops + harvest_tec$stadecoupedf <- "rec" # DO NOTHING ELSE FOR NOW # TODO: ADD OTHER MANAGEMENT, E.G. FERTILIZATION @@ -540,6 +546,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { clim_list[[clim]] <- read.table(met_file) } clim_run <- do.call("rbind", clim_list) + clim_run[,8] <- clim_run[,8]/2 write.table(clim_run, file.path(usmdirs[usmi], "climat.txt"), col.names = FALSE, row.names = FALSE) } diff --git a/models/stics/inst/crop_plt.xml b/models/stics/inst/crop_plt.xml index a0596782441..c460cdbf48e 100644 --- a/models/stics/inst/crop_plt.xml +++ b/models/stics/inst/crop_plt.xml @@ -51,8 +51,8 @@ 7.00000 - 274 - 6.50000 + 284 + 1.50000 10.00000 @@ -260,7 +260,7 @@ 12.440 - 480.000 + 1480.000 0.05000 @longsperac@ @@ -78,7 +86,6 @@ 0.00000 - 0 -999 + -999 -999 -999 @@ -111,10 +119,7 @@ 200.00000 - 140.00000 - 0.50000 - 0.02000 - 0.70000 + 0.00000 0.00000 25.00000 30.00000 @@ -131,24 +136,20 @@ 0.0000 @@ -161,11 +162,8 @@ - 1 20 - 10 15 2.6 2.7 @@ -186,25 +182,57 @@ 25.00000 - 377.000 246.000 - 0.35500 0.00000 0.00000 + + + + 0.879 @@ -237,21 +259,23 @@ - 0.00000 - 0.00000 - 0.01757 - 0.01100 - 0.55000 0.00000 rec 0.00000 0.30000 - + @@ -294,21 +344,18 @@ 2 - -6.00000 -20.00000 @@ -321,14 +368,11 @@ - 13.9 - 2.2 0.90000 0.15000 0.60000 0.70000 0.40000 - 0.00800 0.00500 - @@ -68,7 +68,6 @@ 0.08150 0.20000 0.70000 - 0.70000 0.70000 @@ -78,9 +77,7 @@ 1.40000 0.50000 - 0.02000 0.04500 - 0.02000 10.00000 20.00000 0.00000 @@ -103,10 +100,13 @@ 0.10300 12.00000 15.00000 - 0.65000 - 0.00060 - 0.02720 - 0.01670 + 0.0007 + 0.02519 + 0.015 + 0.11200 + 8.50000 + 0.06000 + 11.00000 0.10500 5.50000 8.50000 @@ -136,7 +136,7 @@ 47 25 148 - 3.0 + 3.0 - + + - - + + 0.00100 + 0.10000 + + + + + - - - - - - - - - - - - - - - + + + + + + + + + + + + 1 + + + - - - - - - - + + + + + diff --git a/models/stics/inst/pecan_sta.xml b/models/stics/inst/pecan_sta.xml index a10b4f2773c..01cd3f9c2a5 100644 --- a/models/stics/inst/pecan_sta.xml +++ b/models/stics/inst/pecan_sta.xml @@ -1,65 +1,86 @@ - - - 2.50000 - 0.00000 - 46.50000 - 1000.00000 - 20.000000 - - - + + + 50.00000 + 0.23000 + 0.18000 + 0.62000 + 1.00000 + + + + + 0.70000 + 6.00000 + 0.50000 + 0.16000 + 0.00400 + 0.59000 + + + + + diff --git a/models/stics/inst/pecan_tec.xml b/models/stics/inst/pecan_tec.xml index c333c76f133..c9bfc287a25 100644 --- a/models/stics/inst/pecan_tec.xml +++ b/models/stics/inst/pecan_tec.xml @@ -1,319 +1,361 @@ - - + + - - - - - - - - - - - - - - - - - - - - - - 999 - 2.00000 - 100.00000 - 1 - - - - + + + + + + + + + + + + + + + + + + + + + 0 + 0.00 + 6.00 + + + + + 0 + 5.50000 + 450 + 1 + + + + - 999 - - + + 999 + 442 + 999 + 999 + 999 + 999 + 999 + 999 + 999 + + + + + 1.00000 - - - - 1 - 0.11000 - - + + + + 20.00000 + + + 0.00000 + + + + + + + + + + 180 + 80 + 1 + + + + + + 245 + straw+roots + + + - - - + + + + + + + + - - - + + 1.10000 + 1.30000 + 0.00500 + 0.05000 + + + + + + + + diff --git a/models/stics/inst/sols.xml b/models/stics/inst/sols.xml index fa606e236e3..0a0e3d5d0d2 100644 --- a/models/stics/inst/sols.xml +++ b/models/stics/inst/sols.xml @@ -1,121 +1,122 @@ - - - - 17.0 - 0.1500 - 15.0000 - 1.0000 - 6.5000 - 0.0100 - 0.2000 - 4.0000 - 0.0000 - 200.0000 - 50.0000 - 0.5000 - 60.0000 - 5.0000 - 0.0100 - 0.0000 - 0.3300 - - - - - - - - - - - - - - - - - - - 20.00 - 16.00 - 4.00 - 1.20 - 0.00 - 1 - 50.00 - 5 - - - 20.00 - 15.30 - 9.30 - 1.30 - 0.00 - 1 - 50.00 - 2 - - - 60.00 - 15.30 - 9.30 - 1.30 - 0.00 - 1 - 50.00 - 1 - - - 0.00 - 0.00 - 0.00 - 0.00 - 0.00 - 1 - 50.00 - 1 - - - 0.00 - 0.00 - 0.00 - 0.00 - 0.00 - 1 - 50.00 - 1 - + + + + 30.2 + 0.2700 + 40.0000 + 0.0000 + 7.0000 + 0.2000 + 0.1700 + 12.0000 + 0.0000 + 200.0000 + 50.0000 + 0.5000 + 60.0000 + 5.0000 + 0.0100 + 0.0000 + 0.65000 + 0.3300 + + + + + + + + + + + + + + + + + + + 20.00 + 46.80 + 26.20 + 1.08 + 0.00 + 1 + 50.00 + 10 + + + 20.00 + 46.40 + 27.40 + 1.09 + 0.00 + 1 + 50.00 + 10 + + + 20.00 + 48.50 + 29.10 + 1.02 + 0.00 + 1 + 50.00 + 10 + + + 20.00 + 50.10 + 25.50 + 0.99 + 0.00 + 1 + 50.00 + 10 + + + 20.00 + 50.10 + 25.50 + 0.99 + 0.00 + 1 + 50.00 + 10 + diff --git a/models/stics/inst/usms.xml b/models/stics/inst/usms.xml index 74792d30b5c..5b927fadb77 100644 --- a/models/stics/inst/usms.xml +++ b/models/stics/inst/usms.xml @@ -1,16 +1,16 @@ - - - + + + @DATEBUT@ @DATEFIN@ - @FINIT@ + @FINIT@ @NOMSOL@ @FSTATION@ @FCLIM1@ - @FCLIM2@ + @FCLIM2@ @CULTUREAN@ @NBPLANTES@ - @CODESIMUL@ + @CODESIMUL@ @FPLT1@ @FTEC1@ @@ -20,6 +20,6 @@ @FPLT2@ @FTEC2@ @FLAI2@ - - + + From c1d9a49f29ebc323819f11d46e9e5c956dfc2289 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 14:59:05 +0200 Subject: [PATCH 54/68] further updates for newer model version --- models/stics/R/write.config.STICS.R | 85 ++++++++++++++--------------- models/stics/inst/prof.mod | 2 +- 2 files changed, 42 insertions(+), 45 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 1ba0f11b666..b4995b12d4b 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -217,7 +217,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { codeplante <- 'fou' codeperenne <- 2 }else{ - codeplante <- substr(names(trait.values)[pft],1,3) + codeplante <- base::substr(names(trait.values)[pft],1,3) codeperenne <- 1 } codebfroid <- 2 # vernalization requirement, hardcoding for now, 2==yes @@ -732,6 +732,21 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { SticsRFiles::set_param_xml(plant_file, "croirac", pft.traits[which(pft.names == "croirac")], overwrite = TRUE) } + # extinction coefficient connecting LAI to crop height + if ("LAI2height" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "khaut", pft.traits[which(pft.names == "LAI2height")], overwrite = TRUE) + } + + # average root radius + if ("rayon" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "rayon", pft.traits[which(pft.names == "rayon")], overwrite = TRUE) + } + + # minimal value for drought stress index + if ("swfacmin" %in% pft.names) { + SticsRFiles::set_param_xml(plant_file, "swfacmin", pft.traits[which(pft.names == "swfacmin")], overwrite = TRUE) + } + # convert xml2txt if(names(trait.values)[pft] != "env"){ SticsRFiles::convert_xml2txt(file = plant_file, javastics = javastics_path) @@ -791,12 +806,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { SticsRFiles::set_param_xml(gen_file, "y0msrac", pft.traits[which(pft.names == "rootmin_harvest")], overwrite = TRUE) } - # move to plt!!! file - # extinction coefficient connecting LAI to crop height - if ("LAI2height" %in% pft.names) { - SticsRFiles::set_param_xml(gen_file, "khaut", pft.traits[which(pft.names == "LAI2height")], overwrite = TRUE) - } - ### Root growth # bulk density of soil below which root growth is reduced due to a lack of soil cohesion (g.cm-3) @@ -825,12 +834,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { if ("lvopt" %in% pft.names) { SticsRFiles::set_param_xml(gen_file, "lvopt", pft.traits[which(pft.names == "lvopt")], overwrite = TRUE) } - - # move to plt!!! file - # average root radius - if ("rayon" %in% pft.names) { - SticsRFiles::set_param_xml(gen_file, "rayon", pft.traits[which(pft.names == "rayon")], overwrite = TRUE) - } # diffusion coefficient of nitrate N in soil at field capacity if ("difN_FC" %in% soil.names) { @@ -891,29 +894,24 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { SticsRFiles::set_param_xml(gen_file, "TREFr", soil_params[which(soil.names == "T_r_ORdecomp")], overwrite = TRUE) } - # not used anymore, or at least not with this name!!! - # initial fraction of soil organic N inactive for mineralisation (= stable SON/ total SON) - if ("FINERT" %in% soil.names) { - SticsRFiles::set_param_xml(gen_file, "FINERT", soil_params[which(soil.names == "FINERT")], overwrite = TRUE) - } - - # not used anymore, or at least not with this name!!! - # relative potential mineralization rate: K2 = fmin1 * exp(- fmin2*argi) / (1+fmin3*calc) - if ("FMIN1" %in% soil.names) { - SticsRFiles::set_param_xml(gen_file, "FMIN1", soil_params[which(soil.names == "FMIN1")], overwrite = TRUE) - } - - # not used anymore, or at least not with this name!!! - # parameter defining the effect of clay on the potential mineralization rate: K2 = fmin1 * exp(-fmin2*argi) / (1+fmin3*calc) - if ("FMIN2" %in% soil.names) { - SticsRFiles::set_param_xml(gen_file, "FMIN2", soil_params[which(soil.names == "FMIN2")], overwrite = TRUE) - } - - # not used anymore, or at least not with this name!!! - # parameter defining the effect of CaCO3 on the potential mineralization rate: K2 = fmin1 * exp(-fmin2*argi) / (1+fmin3*calc) - if ("FMIN3" %in% soil.names) { - SticsRFiles::set_param_xml(gen_file, "FMIN3", soil_params[which(soil.names == "FMIN3")], overwrite = TRUE) - } + # TODO: come back to these + # # not used anymore, or at least not with this name!!! + # # relative potential mineralization rate: K2 = fmin1 * exp(- fmin2*argi) / (1+fmin3*calc) + # if ("FMIN1" %in% soil.names) { + # SticsRFiles::set_param_xml(gen_file, "FMIN1", soil_params[which(soil.names == "FMIN1")], overwrite = TRUE) + # } + # + # # not used anymore, or at least not with this name!!! + # # parameter defining the effect of clay on the potential mineralization rate: K2 = fmin1 * exp(-fmin2*argi) / (1+fmin3*calc) + # if ("FMIN2" %in% soil.names) { + # SticsRFiles::set_param_xml(gen_file, "FMIN2", soil_params[which(soil.names == "FMIN2")], overwrite = TRUE) + # } + # + # # not used anymore, or at least not with this name!!! + # # parameter defining the effect of CaCO3 on the potential mineralization rate: K2 = fmin1 * exp(-fmin2*argi) / (1+fmin3*calc) + # if ("FMIN3" %in% soil.names) { + # SticsRFiles::set_param_xml(gen_file, "FMIN3", soil_params[which(soil.names == "FMIN3")], overwrite = TRUE) + # } # N/C ratio of soil humus if ("Wh" %in% soil.names) { @@ -1165,13 +1163,6 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { }) ### new formulations - - # move to plt!!! file - # minimal value for drought stress index - if ("swfacmin" %in% pft.names) { - SticsRFiles::set_param_xml(newf_file, "swfacmin", pft.traits[which(pft.names == "swfacmin")], overwrite = TRUE) - } - # DO NOTHING ELSE FOR NOW SticsRFiles::convert_xml2txt(file = newf_file, javastics = javastics_path) @@ -1268,6 +1259,12 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { #### I NOW PUT TXT VERSION TO THE MODEL PACKAGE: param.sol #### TODO: revise others to have txt templates directly in the package + # # changed from FINERT to finert and moved to the sols.xml + # # initial fraction of soil organic N inactive for mineralisation (= stable SON/ total SON) + # if ("FINERT" %in% soil.names) { + # SticsRFiles::set_param_xml(gen_file, "finert", soil_params[which(soil.names == "FINERT")], overwrite = TRUE) + # } + sols_file <- file.path(rundir, "param.sol") # cp template sols file (txt) @@ -1577,7 +1574,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # Type of LAI simulation # 0 = culture (LAI calculated by the model), 1 = feuille (LAI forced) - SticsRFiles::set_usm_txt(usm_file, "codesimul", 0, append = FALSE) # hardcode for now + SticsRFiles::set_usm_txt(usm_file, "codesimul", "culture", append = FALSE) # hardcode for now # use optimization # 0 = no; 1 = yes main plant; 2 = yes associated plant diff --git a/models/stics/inst/prof.mod b/models/stics/inst/prof.mod index 00510174f7d..e0b00c6fe59 100644 --- a/models/stics/inst/prof.mod +++ b/models/stics/inst/prof.mod @@ -1,4 +1,4 @@ 2 -tsol(iz) +Chum 10 01 01 2000 From 0fc4b841344338dab7407a127eccd95827663e0d Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 15:11:22 +0200 Subject: [PATCH 55/68] dont change nonstics scripts --- base/db/R/dbfiles.R | 234 +++++++++--------- .../assim.batch/R/pda.generate.externals.R | 17 +- .../assim.batch/man/pda.generate.externals.Rd | 16 +- 3 files changed, 132 insertions(+), 135 deletions(-) diff --git a/base/db/R/dbfiles.R b/base/db/R/dbfiles.R index 4621bfb35fd..d151b655783 100644 --- a/base/db/R/dbfiles.R +++ b/base/db/R/dbfiles.R @@ -43,28 +43,28 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, parentid = NA, con, hostname = PEcAn.remote::fqdn(), allow.conflicting.dates = FALSE, ens = FALSE) { name <- basename(in.path) hostname <- default_hostname(hostname) - + # find mimetype, if it does not exist, it will create one mimetypeid <- get.id("mimetypes", "type_string", mimetype, con, create = TRUE) - + # find appropriate format, create if it does not exist formatid <- get.id( table = "formats", colnames = c("mimetype_id", "name"), - values = c(format(mimetypeid), formatname), + values = c(mimetypeid, formatname), con = con, create = TRUE, dates = TRUE ) - - + + # setup parent part of query if specified if (is.na(parentid)) { parent <- "" } else { parent <- paste0(" AND parent_id=", parentid) } - + # find appropriate input, if not in database, insert new input existing.input <- db.query( query = paste0( @@ -75,7 +75,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, ), con = con ) - + inputid <- NULL if (nrow(existing.input) > 0) { # Convert dates to Date objects and strip all time zones (DB values are timezone-free) @@ -87,7 +87,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } existing.input$start_date <- lubridate::force_tz(time = lubridate::as_date(existing.input$start_date), tzone = "UTC") existing.input$end_date <- lubridate::force_tz(time = lubridate::as_date(existing.input$end_date), tzone = "UTC") - + for (i in seq_len(nrow(existing.input))) { existing.input.i <- existing.input[i, ] if (is.na(existing.input.i$start_date) && is.null(startdate)) { @@ -96,7 +96,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, inputid <- existing.input.i[["id"]] } } - + if (is.null(inputid) && !allow.conflicting.dates) { print(existing.input, digits = 10) PEcAn.logger::logger.error(paste0( @@ -110,7 +110,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, return(NULL) } } - + if (is.null(inputid)) { # Either there was no existing input, or there was but the dates don't match and # allow.conflicting.dates==TRUE. So, insert new input record. @@ -145,7 +145,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, # This is the id that we just registered inserted.id <- db.query(query = cmd, con = con) name.s <- name - + if (is.null(startdate)) { inputid <- db.query( query = paste0( @@ -169,7 +169,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } else { inserted.id <- data.frame(id = inputid) # in the case that inputid is not null then this means that there was an exsiting input } - + if (length(inputid) > 1 && !ens) { PEcAn.logger::logger.warn(paste0( "Multiple input files found matching parameters format_id = ", formatid, @@ -181,17 +181,17 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } else if (ens) { inputid <- inserted.id$id } - + # find appropriate dbfile, if not in database, insert new dbfile dbfile <- dbfile.check(type = "Input", container.id = inputid, con = con, hostname = hostname) - + if (nrow(dbfile) > 0 & !ens) { if (nrow(dbfile) > 1) { print(dbfile) PEcAn.logger::logger.warn("Multiple dbfiles found. Using last.") dbfile <- dbfile[nrow(dbfile), ] } - + if (dbfile$file_name != in.prefix || dbfile$file_path != in.path && !ens) { print(dbfile, digits = 10) PEcAn.logger::logger.error(paste0( @@ -210,7 +210,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, con = con, reuse = TRUE, hostname = hostname ) } - + invisible(list(input.id = inputid, dbfile.id = dbfileid)) } @@ -241,26 +241,26 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetype, formatname, parentid = NA, con, hostname = PEcAn.remote::fqdn(), exact.dates = FALSE, pattern = NULL, return.all = FALSE) { hostname <- default_hostname(hostname) - + mimetypeid <- get.id(table = "mimetypes", colnames = "type_string", values = mimetype, con = con) if (is.null(mimetypeid)) { return(invisible(data.frame())) } - + # find appropriate format formatid <- get.id(table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con) - + if (is.null(formatid)) { invisible(data.frame()) } - + # setup parent part of query if specified if (is.na(parentid)) { parent <- "" } else { parent <- paste0(" AND parent_id=", parentid) } - + # find appropriate input if (exact.dates) { if (!is.null(enddate)) { @@ -287,15 +287,15 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp } } else { # not exact dates inputs <- db.query( - query = paste0( - "SELECT * FROM inputs WHERE site_id=", siteid, - " AND format_id=", formatid, - parent - ), - con = con + query = paste0( + "SELECT * FROM inputs WHERE site_id=", siteid, + " AND format_id=", formatid, + parent + ), + con = con ) } - + if (is.null(inputs) | length(inputs$id) == 0) { return(data.frame()) } else { @@ -303,12 +303,12 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp ## Case where pattern is not NULL inputs <- inputs[grepl(pattern, inputs$name), ] } - + ## parent check when NA # if (is.na(parentid)) { # inputs <- inputs[is.na(inputs$parent_id),] # } - + if (length(inputs$id) > 1) { PEcAn.logger::logger.warn("Found multiple matching inputs. Checking for one with associate files on host machine") print(inputs) @@ -317,7 +317,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp # for(i in seq_len(ni)){ # dbfile[[i]] <- dbfile.check(type = 'Input', container.id = inputs$id[i], con = con, hostname = hostname, machine.check = TRUE) # } - + dbfile <- dbfile.check( type = "Input", @@ -327,8 +327,8 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp machine.check = TRUE, return.all = return.all ) - - + + if (nrow(dbfile) == 0) { ## With the possibility of dbfile.check returning nothing, ## as.data.frame ensures a empty data.frame is returned @@ -336,10 +336,10 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp PEcAn.logger::logger.info("File not found on host machine. Returning Valid input with file associated on different machine if possible") return(as.data.frame(dbfile.check(type = "Input", container.id = inputs$id, con = con, hostname = hostname, machine.check = FALSE))) } - + return(dbfile) } else if (length(inputs$id) == 0) { - + # need this third case here because prent check above can return an empty inputs return(data.frame()) } else { @@ -353,7 +353,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp machine.check = TRUE, return.all = return.all ) - + if (nrow(dbfile) == 0) { ## With the possibility of dbfile.check returning nothing, ## as.data.frame ensures an empty data.frame is returned @@ -361,7 +361,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp PEcAn.logger::logger.info("File not found on host machine. Returning Valid input with file associated on different machine if possible") return(as.data.frame(dbfile.check(type = "Input", container.id = inputs$id, con = con, hostname = hostname, machine.check = FALSE))) } - + return(dbfile) } } @@ -388,28 +388,28 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp ##' } dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate pft pftid <- get.id("pfts", "name", pft, con) if (is.null(pftid)) { PEcAn.logger::logger.severe("Could not find pft, could not store file", filename) } - + mimetypeid <- get.id( table = "mimetypes", colnames = "type_string", values = mimetype, con = con, create = TRUE ) - + # find appropriate format formatid <- get.id( table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con, create = TRUE, dates = TRUE ) - + # find appropriate posterior # NOTE: This is defined but not used # posterior_ids <- get.id("posteriors", "pft_id", pftid, con) - + posteriorid_query <- paste0( "SELECT id FROM posteriors WHERE pft_id=", pftid, " AND format_id=", formatid @@ -426,7 +426,7 @@ dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, ho ) posteriorid <- db.query(posteriorid_query, con)[["id"]] } - + # NOTE: Modified by Alexey Shiklomanov. # I'm not sure how this is supposed to work, but I think it's like this invisible(dbfile.insert( @@ -455,24 +455,24 @@ dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, ho ##' } dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate pft pftid <- get.id(table = "pfts", values = "name", colnames = pft, con = con) if (is.null(pftid)) { invisible(data.frame()) } - + # find appropriate format mimetypeid <- get.id(table = "mimetypes", values = "type_string", colnames = mimetype, con = con) if (is.null(mimetypeid)) { PEcAn.logger::logger.error("mimetype ", mimetype, "does not exist") } formatid <- get.id(table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con) - + if (is.null(formatid)) { invisible(data.frame()) } - + # find appropriate posterior posteriorid <- db.query( query = paste0( @@ -484,7 +484,7 @@ dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PE if (is.null(posteriorid)) { invisible(data.frame()) } - + invisible(dbfile.check(type = "Posterior", container.id = posteriorid, con = con, hostname = hostname)) } @@ -509,14 +509,14 @@ dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PE ##' } dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + if (substr(in.path, 1, 1) != "/") { PEcAn.logger::logger.error("path to dbfiles:", in.path, " is not a valid full path") } - + # find appropriate host hostid <- get.id(table = "machines", colnames = "hostname", values = hostname, con = con, create = TRUE, dates = TRUE) - + # Query for existing dbfile record with same file_name, file_path, machine_id , # container_type, and container_id. dbfile <- invisible(db.query( @@ -528,10 +528,10 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn ), con = con )) - + if (nrow(dbfile) == 0) { # If no existing record, insert one - + insert_result <- db.query( query = paste0( "INSERT INTO dbfiles ", @@ -541,7 +541,7 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn ), con = con ) - + file.id <- insert_result[["id"]] } else if (!reuse) { # If there is an existing record but reuse==FALSE, return NA. @@ -559,7 +559,7 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn file.id <- dbfile[["id"]] } } - + # Return the new dbfile ID, or the one that existed already (reuse==T), or NA (reuse==F) return(file.id) } @@ -591,34 +591,34 @@ dbfile.check <- function(type, container.id, con, machine.check = TRUE, return.all = FALSE) { type <- match.arg(type, c("Input", "Posterior", "Model")) - + hostname <- default_hostname(hostname) - + # find appropriate host hostid <- get.id(table = "machines", colnames = "hostname", values = hostname, con = con) if (is.null(hostid)) { return(data.frame()) } - + dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::filter( .data$container_type == !!type, .data$container_id %in% !!container.id ) - + if (machine.check) { dbfiles <- dbfiles %>% dplyr::filter(.data$machine_id == !!hostid) } - + dbfiles <- dplyr::collect(dbfiles) - + if (nrow(dbfiles) > 1 && !return.all) { PEcAn.logger::logger.warn("Multiple Valid Files found on host machine. Returning last updated record.") dbfiles <- dbfiles %>% dplyr::filter(.data$updated_at == max(.data$updated_at)) } - + dbfiles } @@ -643,9 +643,9 @@ dbfile.check <- function(type, container.id, con, ##' } dbfile.file <- function(type, id, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + files <- dbfile.check(type = type, container.id = id, con = con, hostname = hostname) - + if (nrow(files) > 1) { PEcAn.logger::logger.warn("multiple files found for", id, "returned; using the first one found") invisible(file.path(files[1, "file_path"], files[1, "file_name"])) @@ -667,13 +667,13 @@ dbfile.file <- function(type, id, con, hostname = PEcAn.remote::fqdn()) { ##' } dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate host hostid <- db.query(query = paste0("SELECT id FROM machines WHERE hostname='", hostname, "'"), con = con)[["id"]] if (is.null(hostid)) { invisible(NA) } - + # find file file_name <- basename(file) file_path <- dirname(file) @@ -686,7 +686,7 @@ dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { ), con = con ) - + if (nrow(ids) > 1) { PEcAn.logger::logger.warn("multiple ids found for", file, "returned; using the first one found") invisible(ids[1, "container_id"]) @@ -727,58 +727,58 @@ dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = FALSE) { - - + + # create nulls for file movement and error info error <- 0 files.sym <- 0 files.changed <- 0 files.reg <- 0 files.indb <- 0 - + # check for file type and update to make it *.file type if (file.type != "clim" | file.type != "nc") { PEcAn.logger::logger.error("File type not supported by move at this time. Currently only supports NC and CLIM files") error <- 1 } file.pattern <- paste0("*.", file.type) - - - + + + # create new directory if it doesn't exist if (!dir.exists(new.dir)) { dir.create(new.dir) } - - + + # check to make sure both directories exist if (!dir.exists(old.dir)) { PEcAn.logger::logger.error("Old File directory does not exist. Please enter valid file path") error <- 1 } - + if (!dir.exists(new.dir)) { PEcAn.logger::logger.error("New File directory does not exist. Please enter valid file path") error <- 1 } - + if (basename(new.dir) != basename(old.dir)) { PEcAn.logger::logger.error("Basenames of files do not match") } - + # list files in the old directory old.files <- list.files(path = old.dir, pattern = file.pattern) - + # check to make sure there are files if (length(old.files) == 0) { PEcAn.logger::logger.warn("No files found") error <- 1 } - + # create full file path full.old.file <- file.path(old.dir, old.files) - - + + ### Get BETY information ### con <- db.open( params = list( @@ -788,46 +788,46 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F user = "bety", password = "bety") ) - + # get matching dbfiles from BETY dbfile.path <- dirname(full.old.file) dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::collect() %>% dplyr::filter(.data$file_name %in% basename(full.old.file)) %>% dplyr::filter(.data$file_path %in% dbfile.path) - - + + # if there are matching db files if (dim(dbfiles)[1] > 0) { - + # Check to make sure files line up if (dim(dbfiles)[1] != length(full.old.file)) { PEcAn.logger::logger.warn("Files to be moved don't match up with BETY files, only moving the files that match") - + # IF DB FILES AND FULL FILES DONT MATCH, remove those not in BETY - will take care of the rest below index <- which(basename(full.old.file) %in% dbfiles$file_name) index1 <- seq(1, length(full.old.file)) check <- index1[-which(index1 %in% index)] full.old.file <- full.old.file[-check] - + # record the number of files that are being moved files.changed <- length(full.old.file) } - + # Check to make sure the files line up if (dim(dbfiles)[1] != length(full.old.file)) { PEcAn.logger::logger.error("Files to be moved don't match up with BETY files, canceling move") error <- 1 } - - + + # Make sure the files line up dbfiles <- dbfiles[order(dbfiles$file_name), ] full.old.file <- sort(full.old.file) - + # Record number of files moved and changed in BETY files.indb <- dim(dbfiles)[1] - + # Move files and update BETY if (error == 0) { for (i in 1:length(full.old.file)) { @@ -836,40 +836,40 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error if statement } # end dbfile loop - - + + # if there are files that are in the folder but not in BETY, we can either register them or not if (dim(dbfiles)[1] == 0 | files.changed > 0) { - + # Recheck what files are in the directory since others may have been moved above old.files <- list.files(path = old.dir, pattern = file.pattern) - + # Recreate full file path full.old.file <- file.path(old.dir, old.files) - - + + # Error check again to make sure there aren't any matching dbfiles dbfile.path <- dirname(full.old.file) dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::collect() %>% dplyr::filter(.data$file_name %in% basename(full.old.file)) %>% dplyr::filter(.data$file_path %in% dbfile.path) - + if (dim(dbfiles)[1] > 0) { PEcAn.logger::logger.error("There are still dbfiles matching these files! Canceling link or registration") error <- 1 } - - + + if (error == 0 & register == TRUE) { - + # Record how many files are being registered to BETY files.reg <- length(full.old.file) - + for (i in 1:length(full.old.file)) { file_path <- dirname(full.old.file[i]) file_name <- basename(full.old.file[i]) - + if (file.type == "nc") { mimetype <- "application/x-netcdf" formatname <- "CF Meteorology application" @@ -881,8 +881,8 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F else { PEcAn.logger::logger.error("File Type is currently not supported") } - - + + dbfile.input.insert( in.path = file_path, in.prefix = file_name, @@ -900,26 +900,26 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error loop } # end register == TRUE - + if (error == 0 & register == FALSE) { # Create file path for symbolic link full.new.file <- file.path(new.dir, old.files) - + # Record number of files that will have a symbolic link made files.sym <- length(full.new.file) - + # Line up files full.new.file <- sort(full.new.file) full.old.file <- sort(full.old.file) - + # Check to make sure the files are the same length if (length(full.new.file) != length(full.old.file)) { PEcAn.logger::logger.error("Files to be moved don't match up with BETY. Canceling Move") error <- 1 } - + # Move file and create symbolic link if there are no errors - + if (error == 0) { for (i in 1:length(full.old.file)) { fs::file_move(full.old.file[i], new.dir) @@ -927,13 +927,13 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error loop } # end Register == FALSE - - + + if (error > 0) { PEcAn.logger::logger.error("There was an error, files were not moved or linked") } - + if (error == 0) { PEcAn.logger::logger.info(paste0(files.changed + files.indb, " files were moved and updated on BETY, ", files.sym, " were moved and had a symbolic link created, and ", files.reg, " files were moved and then registered in BETY")) } -} # end dbfile.move() +} # end dbfile.move() \ No newline at end of file diff --git a/modules/assim.batch/R/pda.generate.externals.R b/modules/assim.batch/R/pda.generate.externals.R index 6d65d19d7c8..a1e6163d2fd 100644 --- a/modules/assim.batch/R/pda.generate.externals.R +++ b/modules/assim.batch/R/pda.generate.externals.R @@ -52,15 +52,14 @@ ##' @param ind.list a named list of vectors (one per pft), where each vector indicates the indices of the parameters on the prior.list targeted in the PDA ##' e.g. ind.list <- list(temperate.deciduous = c(2), temperate.conifer = c(1,2)) ##' @param nknots number of knots you want to train the emulator on -##' -##' example -##' +##' @export +##' @examples +##' \dontrun{ ##' pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, external.formats = TRUE, ##' model_data_diag = TRUE, model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", ##' start_date = "2017-01-01", end_date = "2018-12-31") -##' -##' @author Istem Fer -##' @export +##' } + pda.generate.externals <- function(external.data = FALSE, obs = NULL, varn = NULL, varid = NULL, n_eff = NULL, align_method = "match_timestep", par = NULL, model_data_diag = FALSE, model.out = NULL, start_date = NULL, end_date = NULL, external.formats = FALSE, @@ -158,8 +157,8 @@ pda.generate.externals <- function(external.data = FALSE, obs = NULL, varn = } pda.externals$external.knots <- external.knots - - + + ##################### model & data alignment diagnostics ##################### if(model_data_diag){ if(is.null(obs) & is.null(model.out) & is.null(varn) & is.null(start_date) & is.null(end_date)){ @@ -188,4 +187,4 @@ pda.generate.externals <- function(external.data = FALSE, obs = NULL, varn = pda.externals$model_data_diag <- model_data_diag return(pda.externals) -} +} \ No newline at end of file diff --git a/modules/assim.batch/man/pda.generate.externals.Rd b/modules/assim.batch/man/pda.generate.externals.Rd index f0917ad29f3..4a649bd4a16 100644 --- a/modules/assim.batch/man/pda.generate.externals.Rd +++ b/modules/assim.batch/man/pda.generate.externals.Rd @@ -95,18 +95,16 @@ If not NULL these are used, otherwise knots will be generated using prior.list, \item{ind.list}{a named list of vectors (one per pft), where each vector indicates the indices of the parameters on the prior.list targeted in the PDA e.g. ind.list <- list(temperate.deciduous = c(2), temperate.conifer = c(1,2))} -\item{nknots}{number of knots you want to train the emulator on - -example - -pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, external.formats = TRUE, -model_data_diag = TRUE, model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", -start_date = "2017-01-01", end_date = "2018-12-31")} +\item{nknots}{number of knots you want to train the emulator on} } \description{ This is a helper function for preparing PDA external objects, but it doesn't cover all the cases yet, use it with care You can use this function just to generate either one of the external.* PDA objects, but note that some args cannot be blank depending on what you aim to generate } -\author{ -Istem Fer +\examples{ +\dontrun{ +pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, external.formats = TRUE, + model_data_diag = TRUE, model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", + start_date = "2017-01-01", end_date = "2018-12-31") +} } From a855426ffa30b786a1325b1355eb77881adb8211 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 15:16:35 +0200 Subject: [PATCH 56/68] revert one more --- modules/assim.batch/R/minimize.GP.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/assim.batch/R/minimize.GP.R b/modules/assim.batch/R/minimize.GP.R index db2a8fced44..1c54b3e3aa6 100644 --- a/modules/assim.batch/R/minimize.GP.R +++ b/modules/assim.batch/R/minimize.GP.R @@ -343,7 +343,7 @@ mcmc.GP <- function(gp, x0, nmcmc, rng, format = "lin", mix = "joint", splinefun newllp <- pda.calc.llik.par(settings, n.of.obs, newSS, hyper.pars) - ynew <- get_y(newSS, xnew, llik.fn, priors, currllp) + ynew <- get_y(newSS, xnew, llik.fn, priors, newllp) if (is.accepted(ycurr, ynew)) { xcurr <- xnew currSS <- newSS From dde53c1a2f4c83742718fcfc01619db1088aa334 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 15:34:22 +0200 Subject: [PATCH 57/68] update some IC names --- models/stics/R/write.config.STICS.R | 18 +++++++++--------- modules/assim.batch/R/pda.generate.externals.R | 5 +++-- .../assim.batch/man/pda.generate.externals.Rd | 5 +++-- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index b4995b12d4b..64731c61f84 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -1289,28 +1289,28 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { sapply(1:5, function(x) SticsRFiles::set_soil_txt(file = sols_file, param="epc", value=20, layer = x)) - # water_concentration_at_field_capacity - hccf <- ncdf4::ncvar_get(ic_nc, "water_concentration_at_field_capacity") + # volume_fraction_of_water_in_soil_at_field_capacity + hccf <- ncdf4::ncvar_get(ic_nc, "volume_fraction_of_water_in_soil_at_field_capacity") hccf <- round(hccf*100, digits = 2) sapply(seq_along(hccf), function(x) SticsRFiles::set_soil_txt(file = sols_file, param="hccf", value=hccf[x], layer = x)) - # water_concentration_at_wilting_point - hminf <- ncdf4::ncvar_get(ic_nc, "water_concentration_at_wilting_point") + # volume_fraction_of_condensed_water_in_soil_at_wilting_point + hminf <- ncdf4::ncvar_get(ic_nc, "volume_fraction_of_condensed_water_in_soil_at_wilting_point") hminf <- round(hminf*100, digits = 2) sapply(seq_along(hminf), function(x) SticsRFiles::set_soil_txt(file = sols_file, param="hminf", value=hminf[x], layer = x)) # soil_organic_nitrogen_content - Norg <- ncdf4::ncvar_get(ic_nc, "soil_nitrogen_content") + Norg <- ncdf4::ncvar_get(ic_nc, "soil_organic_nitrogen_content") Norg <- round(Norg[1]*100, digits = 2) # STICS uses 1 Norg value SticsRFiles::set_soil_txt(file = sols_file, param="Norg", value=Norg) - # soil_clay_content - argi <- ncdf4::ncvar_get(ic_nc, "soil_clay_content") + # mass_fraction_of_clay_in_soil + argi <- ncdf4::ncvar_get(ic_nc, "mass_fraction_of_clay_in_soil") argi <- round(argi[1]*100, digits = 0) # STICS uses 1 argi value SticsRFiles::set_soil_txt(file = sols_file, param="argi", value=argi) - # soil_bulk_density (kg m-3 --> g cm-3) - DAF <- ncdf4::ncvar_get(ic_nc, "soil_bulk_density") + # soil_density (kg m-3 --> g cm-3) + DAF <- ncdf4::ncvar_get(ic_nc, "soil_density") DAF <- round(udunits2::ud.convert(DAF, "kg m-3", "g cm-3"), digits = 1) sapply(seq_along(DAF), function(x) SticsRFiles::set_soil_txt(file = sols_file, param="DAF", value=DAF[x], layer = x)) diff --git a/modules/assim.batch/R/pda.generate.externals.R b/modules/assim.batch/R/pda.generate.externals.R index a1e6163d2fd..46de30818df 100644 --- a/modules/assim.batch/R/pda.generate.externals.R +++ b/modules/assim.batch/R/pda.generate.externals.R @@ -55,8 +55,9 @@ ##' @export ##' @examples ##' \dontrun{ -##' pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, external.formats = TRUE, -##' model_data_diag = TRUE, model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", +##' pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, +##' external.formats = TRUE, model_data_diag = TRUE, +##' model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", ##' start_date = "2017-01-01", end_date = "2018-12-31") ##' } diff --git a/modules/assim.batch/man/pda.generate.externals.Rd b/modules/assim.batch/man/pda.generate.externals.Rd index 4a649bd4a16..20d06c5c322 100644 --- a/modules/assim.batch/man/pda.generate.externals.Rd +++ b/modules/assim.batch/man/pda.generate.externals.Rd @@ -103,8 +103,9 @@ You can use this function just to generate either one of the external.* PDA obje } \examples{ \dontrun{ -pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, external.formats = TRUE, - model_data_diag = TRUE, model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", +pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, + external.formats = TRUE, model_data_diag = TRUE, + model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", start_date = "2017-01-01", end_date = "2018-12-31") } } From c643a62474274c442ca9522105db55c31dd329cb Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 15:38:06 +0200 Subject: [PATCH 58/68] remove udunits calls --- models/stics/R/model2netcdf.STICS.R | 8 ++++---- models/stics/R/write.config.STICS.R | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/models/stics/R/model2netcdf.STICS.R b/models/stics/R/model2netcdf.STICS.R index c1974d03c04..e296d6a450c 100644 --- a/models/stics/R/model2netcdf.STICS.R +++ b/models/stics/R/model2netcdf.STICS.R @@ -60,15 +60,15 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o 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 <- udunits2::ud.convert(stics_output[thisyear, "CO2sol"], "ha-1 day-1", "m-2 s-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 <- udunits2::ud.convert(stics_output[thisyear, "dltams.n."], "ton", "kg") * 0.48 # ton to kgC + 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 <- udunits2::ud.convert(stics_output[thisyear, "dltaremobil"], "ton", "kg") * 0.48 # ton to kgC + 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 @@ -79,7 +79,7 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o ## should be roughly equal to this: #diff(stics_output[thisyear, "masec.n."])+ diff(stics_output[thisyear, "msrac.n."]) # t.ha-1 - NPP <- udunits2::ud.convert(NPP, "ha-1 day-1", "m-2 s-1") # kg C m-2 s-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) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index 64731c61f84..aa0e15c4f3c 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -542,14 +542,14 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # maximum SLA (specific leaf area) of green leaves (cm2 g-1) if ("SLAMAX" %in% pft.names) { slamax <- pft.traits[which(pft.names == "SLAMAX")] - slamax <- udunits2::ud.convert(udunits2::ud.convert(slamax, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 + slamax <- PEcAn.utils::ud_convert(PEcAn.utils::ud_convert(slamax, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 SticsRFiles::set_param_xml(plant_file, "slamax", slamax, overwrite = TRUE) } # minimum SLA (specific leaf area) of green leaves (cm2 g-1) if ("SLAMIN" %in% pft.names) { slamin <- pft.traits[which(pft.names == "SLAMIN")] - slamin <- udunits2::ud.convert(udunits2::ud.convert(slamin, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 + slamin <- PEcAn.utils::ud_convert(PEcAn.utils::ud_convert(slamin, "m2", "cm2"), "kg-1", "g-1") # m2 kg-1 to cm2 g-1 SticsRFiles::set_param_xml(plant_file, "slamin", slamin, overwrite = TRUE) } @@ -598,7 +598,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # longsperac - specific root length (cm g-1) if ("SRL" %in% pft.names) { - srl_val <- udunits2::ud.convert(pft.traits[which(pft.names == "SRL")], "m", "cm") + srl_val <- PEcAn.utils::ud_convert(pft.traits[which(pft.names == "SRL")], "m", "cm") SticsRFiles::set_param_xml(plant_file, "longsperac", srl_val, overwrite = TRUE) } @@ -1207,24 +1207,24 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # initial aerial biomass (kg m-2 --> t ha-1) masec0 <- ncdf4::ncvar_get(ic_nc, "AGB") - SticsRFiles::set_param_xml(file = ini_file, param = "masec0", values = udunits2::ud.convert(masec0, "kg m-2", "t ha-1"), select = "plante", select_value = "1", overwrite = TRUE) + SticsRFiles::set_param_xml(file = ini_file, param = "masec0", values = PEcAn.utils::ud_convert(masec0, "kg m-2", "t ha-1"), select = "plante", select_value = "1", overwrite = TRUE) # initial depth of root apex of the crop (m --> cm) zrac0 <- ncdf4::ncvar_get(ic_nc, "rooting_depth") if(zrac0 < 0.2) zrac0 <- 0.2 - SticsRFiles::set_param_xml(file = ini_file, param = "zrac0", values = udunits2::ud.convert(zrac0, "m", "cm"), select = "plante", select_value = "1", overwrite = TRUE) + SticsRFiles::set_param_xml(file = ini_file, param = "zrac0", values = PEcAn.utils::ud_convert(zrac0, "m", "cm"), select = "plante", select_value = "1", overwrite = TRUE) # initial grain dry weight - haven't started any simulations from this stage yet # SticsRFiles::set_param_xml(file = ini_file, param = "magrain0", values = 0, select = "plante", select_value = "1", overwrite = TRUE) # initial N amount in the plant (kg m-2 --> kg ha-1) QNplante0 <- ncdf4::ncvar_get(ic_nc, "plant_nitrogen_content") - SticsRFiles::set_param_xml(file = ini_file, param = "QNplante0", values = udunits2::ud.convert(QNplante0, "kg m-2", "kg ha-1"), select = "plante", select_value = "1", overwrite = TRUE) + SticsRFiles::set_param_xml(file = ini_file, param = "QNplante0", values = PEcAn.utils::ud_convert(QNplante0, "kg m-2", "kg ha-1"), select = "plante", select_value = "1", overwrite = TRUE) # Not anymore # initial reserve of biomass (kg m-2 --> t ha-1) #resperenne0 <- ncdf4::ncvar_get(ic_nc, "reserve_biomass") - #SticsRFiles::set_param_xml(file = ini_file, param = "resperenne0", values = udunits2::ud.convert(resperenne0, "kg m-2", "t ha-1"), select = "plante", select_value = "1", overwrite = TRUE) + #SticsRFiles::set_param_xml(file = ini_file, param = "resperenne0", values = PEcAn.utils::ud_convert(resperenne0, "kg m-2", "t ha-1"), select = "plante", select_value = "1", overwrite = TRUE) # initial root density in each of the five soil layers densinitial <- ncdf4::ncvar_get(ic_nc, "root_density") @@ -1311,7 +1311,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { # soil_density (kg m-3 --> g cm-3) DAF <- ncdf4::ncvar_get(ic_nc, "soil_density") - DAF <- round(udunits2::ud.convert(DAF, "kg m-3", "g cm-3"), digits = 1) + DAF <- round(PEcAn.utils::ud_convert(DAF, "kg m-3", "g cm-3"), digits = 1) sapply(seq_along(DAF), function(x) SticsRFiles::set_soil_txt(file = sols_file, param="DAF", value=DAF[x], layer = x)) # c2n_humus From 28dc8acfb87a6fc88b3acc04b15c0f30a7bc82ff Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 16:08:01 +0200 Subject: [PATCH 59/68] working with namespace calls --- models/stics/DESCRIPTION | 5 +++-- models/stics/R/met2model.STICS.R | 2 +- models/stics/R/model2netcdf.STICS.R | 2 +- models/stics/R/write.config.STICS.R | 4 ++-- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/models/stics/DESCRIPTION b/models/stics/DESCRIPTION index f37764650f0..bb2db15e4c0 100644 --- a/models/stics/DESCRIPTION +++ b/models/stics/DESCRIPTION @@ -10,17 +10,18 @@ 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: testthat (>= 1.0.2) +Remotes: + github::SticsRPacks/SticsRPacks SystemRequirements: STICS OS_type: unix License: BSD_3_clause + file LICENSE diff --git a/models/stics/R/met2model.STICS.R b/models/stics/R/met2model.STICS.R index c4586fbc4af..127aba01a80 100644 --- a/models/stics/R/met2model.STICS.R +++ b/models/stics/R/met2model.STICS.R @@ -131,7 +131,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date, # column 6: minimum temperature (°C) Tair <- ncdf4::ncvar_get(nc, "air_temperature") ## in Kelvin Tair <- Tair[ydays %in% simdays] - Tair_C <- udunits2::ud.convert(Tair, "K", "degC") + 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 diff --git a/models/stics/R/model2netcdf.STICS.R b/models/stics/R/model2netcdf.STICS.R index e296d6a450c..e7e8409232e 100644 --- a/models/stics/R/model2netcdf.STICS.R +++ b/models/stics/R/model2netcdf.STICS.R @@ -29,7 +29,7 @@ model2netcdf.STICS <- function(outdir, sitelat, sitelon, start_date, end_date, o out_files <- list.files(outdir) stics_out_file <- file.path(outdir, out_files[grepl("mod_s.*", out_files)]) - stics_output <- lapply(stics_out_file, read.table, header = TRUE, 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]), ] diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index aa0e15c4f3c..c3bedf41cde 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -1678,10 +1678,10 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { for(clim in seq_along(usm_years)){ # currently assuming only first year file has been passed to the settings, modify met2model if changing the logic met_file <- gsub(paste0(lubridate::year(settings$run$start.date), ".climate"), paste0(usm_years[clim], ".climate"), met_path) - clim_list[[clim]] <- read.table(met_file) + clim_list[[clim]] <- utils::read.table(met_file) } clim_run <- do.call("rbind", clim_list) - write.table(clim_run, file.path(usmdirs[usmi], "climat.txt"), col.names = FALSE, row.names = FALSE) + utils::write.table(clim_run, file.path(usmdirs[usmi], "climat.txt"), col.names = FALSE, row.names = FALSE) } From d3fe9a6667bbd38fec1dbf76ad8bcfd238364408 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 16:15:38 +0200 Subject: [PATCH 60/68] install separately --- models/stics/DESCRIPTION | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/models/stics/DESCRIPTION b/models/stics/DESCRIPTION index bb2db15e4c0..6b4ef911fb4 100644 --- a/models/stics/DESCRIPTION +++ b/models/stics/DESCRIPTION @@ -17,11 +17,12 @@ Imports: lubridate, ncdf4, XML, - dplyr + dplyr, Suggests: testthat (>= 1.0.2) Remotes: - github::SticsRPacks/SticsRPacks + github::SticsRPacks/SticsRFiles, + github::SticsRPacks/SticsOnR SystemRequirements: STICS OS_type: unix License: BSD_3_clause + file LICENSE From 4d592e1127a7976b29ae4d88af895c649580865e Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 17:06:44 +0200 Subject: [PATCH 61/68] update depends --- docker/depends/pecan_deps_from_github.txt | 2 ++ docker/depends/pecan_package_dependencies.csv | 3 +-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/docker/depends/pecan_deps_from_github.txt b/docker/depends/pecan_deps_from_github.txt index 058170f37f0..2cbeb0a854f 100644 --- a/docker/depends/pecan_deps_from_github.txt +++ b/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 diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv index 2ffb26841a7..e5e9568c811 100644 --- a/docker/depends/pecan_package_dependencies.csv +++ b/docker/depends/pecan_package_dependencies.csv @@ -113,6 +113,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 @@ -270,7 +271,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 @@ -428,7 +428,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 From 43b556deba84333ab5a7f7f511e9a95d22b69a8b Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 17:17:22 +0200 Subject: [PATCH 62/68] working wit CI --- models/stics/DESCRIPTION | 1 + 1 file changed, 1 insertion(+) diff --git a/models/stics/DESCRIPTION b/models/stics/DESCRIPTION index 6b4ef911fb4..5e2cfbd1255 100644 --- a/models/stics/DESCRIPTION +++ b/models/stics/DESCRIPTION @@ -19,6 +19,7 @@ Imports: XML, dplyr, Suggests: + SticsRFiles, testthat (>= 1.0.2) Remotes: github::SticsRPacks/SticsRFiles, From 611c34a40c0c4c03ff0a015697c956767774effe Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 17:22:46 +0200 Subject: [PATCH 63/68] forgot to add these --- Makefile.depends | 2 +- modules/assim.batch/R/pda.generate.externals.R | 9 +++++---- modules/assim.batch/man/pda.generate.externals.Rd | 9 +++++---- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/Makefile.depends b/Makefile.depends index bc5143e0671..563b64c4107 100644 --- a/Makefile.depends +++ b/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 diff --git a/modules/assim.batch/R/pda.generate.externals.R b/modules/assim.batch/R/pda.generate.externals.R index 46de30818df..b20633eedf1 100644 --- a/modules/assim.batch/R/pda.generate.externals.R +++ b/modules/assim.batch/R/pda.generate.externals.R @@ -55,10 +55,11 @@ ##' @export ##' @examples ##' \dontrun{ -##' pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, -##' external.formats = TRUE, model_data_diag = TRUE, -##' model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", -##' start_date = "2017-01-01", end_date = "2018-12-31") +##' pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, +##' varn = "NEE", varid = 297, n_eff = 106.9386, +##' external.formats = TRUE, model_data_diag = TRUE, +##' model.out = "/tmp/out/outdir", +##' start_date = "2017-01-01", end_date = "2018-12-31") ##' } pda.generate.externals <- function(external.data = FALSE, obs = NULL, varn = NULL, varid = NULL, n_eff = NULL, align_method = "match_timestep", par = NULL, diff --git a/modules/assim.batch/man/pda.generate.externals.Rd b/modules/assim.batch/man/pda.generate.externals.Rd index 20d06c5c322..fc5b3ce409c 100644 --- a/modules/assim.batch/man/pda.generate.externals.Rd +++ b/modules/assim.batch/man/pda.generate.externals.Rd @@ -103,9 +103,10 @@ You can use this function just to generate either one of the external.* PDA obje } \examples{ \dontrun{ -pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, varn = "NEE", varid = 297, n_eff = 106.9386, - external.formats = TRUE, model_data_diag = TRUE, - model.out = "/data/workflows/PEcAn_15000000111/out/15000186876", - start_date = "2017-01-01", end_date = "2018-12-31") +pda.externals <- pda.generate.externals(external.data = TRUE, obs = obs, +varn = "NEE", varid = 297, n_eff = 106.9386, +external.formats = TRUE, model_data_diag = TRUE, +model.out = "/tmp/out/outdir", +start_date = "2017-01-01", end_date = "2018-12-31") } } From 58b01f72ef406df574f02eb7f28ea60f1ccd48bb Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 17:31:47 +0200 Subject: [PATCH 64/68] still one more dependency --- docker/depends/pecan_package_dependencies.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv index e5e9568c811..9b76e20f75c 100644 --- a/docker/depends/pecan_package_dependencies.csv +++ b/docker/depends/pecan_package_dependencies.csv @@ -544,6 +544,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 From 60f77c8498fd3bd87f3d07ecc8dd1801c087cda0 Mon Sep 17 00:00:00 2001 From: istfer Date: Fri, 8 Mar 2024 17:49:49 +0200 Subject: [PATCH 65/68] test actually caught a bug --- models/stics/R/met2model.STICS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/stics/R/met2model.STICS.R b/models/stics/R/met2model.STICS.R index 127aba01a80..43ba5edc16b 100644 --- a/models/stics/R/met2model.STICS.R +++ b/models/stics/R/met2model.STICS.R @@ -157,7 +157,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date, # 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) & !all(is.nan(U)) & all(is.nan(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) From e09c3cada014cc67421ba7f94b0b942e63a2d956 Mon Sep 17 00:00:00 2001 From: Istem Fer Date: Fri, 8 Mar 2024 19:31:01 +0200 Subject: [PATCH 66/68] untouch dbfiles.R --- base/db/R/dbfiles.R | 232 ++++++++++++++++++++++---------------------- 1 file changed, 116 insertions(+), 116 deletions(-) diff --git a/base/db/R/dbfiles.R b/base/db/R/dbfiles.R index d151b655783..e0f9556779b 100644 --- a/base/db/R/dbfiles.R +++ b/base/db/R/dbfiles.R @@ -43,10 +43,10 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, parentid = NA, con, hostname = PEcAn.remote::fqdn(), allow.conflicting.dates = FALSE, ens = FALSE) { name <- basename(in.path) hostname <- default_hostname(hostname) - + # find mimetype, if it does not exist, it will create one mimetypeid <- get.id("mimetypes", "type_string", mimetype, con, create = TRUE) - + # find appropriate format, create if it does not exist formatid <- get.id( table = "formats", @@ -56,15 +56,15 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, create = TRUE, dates = TRUE ) - - + + # setup parent part of query if specified if (is.na(parentid)) { parent <- "" } else { parent <- paste0(" AND parent_id=", parentid) } - + # find appropriate input, if not in database, insert new input existing.input <- db.query( query = paste0( @@ -75,7 +75,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, ), con = con ) - + inputid <- NULL if (nrow(existing.input) > 0) { # Convert dates to Date objects and strip all time zones (DB values are timezone-free) @@ -87,7 +87,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } existing.input$start_date <- lubridate::force_tz(time = lubridate::as_date(existing.input$start_date), tzone = "UTC") existing.input$end_date <- lubridate::force_tz(time = lubridate::as_date(existing.input$end_date), tzone = "UTC") - + for (i in seq_len(nrow(existing.input))) { existing.input.i <- existing.input[i, ] if (is.na(existing.input.i$start_date) && is.null(startdate)) { @@ -96,7 +96,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, inputid <- existing.input.i[["id"]] } } - + if (is.null(inputid) && !allow.conflicting.dates) { print(existing.input, digits = 10) PEcAn.logger::logger.error(paste0( @@ -110,7 +110,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, return(NULL) } } - + if (is.null(inputid)) { # Either there was no existing input, or there was but the dates don't match and # allow.conflicting.dates==TRUE. So, insert new input record. @@ -145,7 +145,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, # This is the id that we just registered inserted.id <- db.query(query = cmd, con = con) name.s <- name - + if (is.null(startdate)) { inputid <- db.query( query = paste0( @@ -169,7 +169,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } else { inserted.id <- data.frame(id = inputid) # in the case that inputid is not null then this means that there was an exsiting input } - + if (length(inputid) > 1 && !ens) { PEcAn.logger::logger.warn(paste0( "Multiple input files found matching parameters format_id = ", formatid, @@ -181,17 +181,17 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, } else if (ens) { inputid <- inserted.id$id } - + # find appropriate dbfile, if not in database, insert new dbfile dbfile <- dbfile.check(type = "Input", container.id = inputid, con = con, hostname = hostname) - + if (nrow(dbfile) > 0 & !ens) { if (nrow(dbfile) > 1) { print(dbfile) PEcAn.logger::logger.warn("Multiple dbfiles found. Using last.") dbfile <- dbfile[nrow(dbfile), ] } - + if (dbfile$file_name != in.prefix || dbfile$file_path != in.path && !ens) { print(dbfile, digits = 10) PEcAn.logger::logger.error(paste0( @@ -210,7 +210,7 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, con = con, reuse = TRUE, hostname = hostname ) } - + invisible(list(input.id = inputid, dbfile.id = dbfileid)) } @@ -241,26 +241,26 @@ dbfile.input.insert <- function(in.path, in.prefix, siteid, startdate, enddate, dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetype, formatname, parentid = NA, con, hostname = PEcAn.remote::fqdn(), exact.dates = FALSE, pattern = NULL, return.all = FALSE) { hostname <- default_hostname(hostname) - + mimetypeid <- get.id(table = "mimetypes", colnames = "type_string", values = mimetype, con = con) if (is.null(mimetypeid)) { return(invisible(data.frame())) } - + # find appropriate format formatid <- get.id(table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con) - + if (is.null(formatid)) { invisible(data.frame()) } - + # setup parent part of query if specified if (is.na(parentid)) { parent <- "" } else { parent <- paste0(" AND parent_id=", parentid) } - + # find appropriate input if (exact.dates) { if (!is.null(enddate)) { @@ -287,15 +287,15 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp } } else { # not exact dates inputs <- db.query( - query = paste0( - "SELECT * FROM inputs WHERE site_id=", siteid, - " AND format_id=", formatid, - parent - ), - con = con + query = paste0( + "SELECT * FROM inputs WHERE site_id=", siteid, + " AND format_id=", formatid, + parent + ), + con = con ) } - + if (is.null(inputs) | length(inputs$id) == 0) { return(data.frame()) } else { @@ -303,12 +303,12 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp ## Case where pattern is not NULL inputs <- inputs[grepl(pattern, inputs$name), ] } - + ## parent check when NA # if (is.na(parentid)) { # inputs <- inputs[is.na(inputs$parent_id),] # } - + if (length(inputs$id) > 1) { PEcAn.logger::logger.warn("Found multiple matching inputs. Checking for one with associate files on host machine") print(inputs) @@ -317,7 +317,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp # for(i in seq_len(ni)){ # dbfile[[i]] <- dbfile.check(type = 'Input', container.id = inputs$id[i], con = con, hostname = hostname, machine.check = TRUE) # } - + dbfile <- dbfile.check( type = "Input", @@ -327,8 +327,8 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp machine.check = TRUE, return.all = return.all ) - - + + if (nrow(dbfile) == 0) { ## With the possibility of dbfile.check returning nothing, ## as.data.frame ensures a empty data.frame is returned @@ -336,10 +336,10 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp PEcAn.logger::logger.info("File not found on host machine. Returning Valid input with file associated on different machine if possible") return(as.data.frame(dbfile.check(type = "Input", container.id = inputs$id, con = con, hostname = hostname, machine.check = FALSE))) } - + return(dbfile) } else if (length(inputs$id) == 0) { - + # need this third case here because prent check above can return an empty inputs return(data.frame()) } else { @@ -353,7 +353,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp machine.check = TRUE, return.all = return.all ) - + if (nrow(dbfile) == 0) { ## With the possibility of dbfile.check returning nothing, ## as.data.frame ensures an empty data.frame is returned @@ -361,7 +361,7 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp PEcAn.logger::logger.info("File not found on host machine. Returning Valid input with file associated on different machine if possible") return(as.data.frame(dbfile.check(type = "Input", container.id = inputs$id, con = con, hostname = hostname, machine.check = FALSE))) } - + return(dbfile) } } @@ -388,28 +388,28 @@ dbfile.input.check <- function(siteid, startdate = NULL, enddate = NULL, mimetyp ##' } dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate pft pftid <- get.id("pfts", "name", pft, con) if (is.null(pftid)) { PEcAn.logger::logger.severe("Could not find pft, could not store file", filename) } - + mimetypeid <- get.id( table = "mimetypes", colnames = "type_string", values = mimetype, con = con, create = TRUE ) - + # find appropriate format formatid <- get.id( table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con, create = TRUE, dates = TRUE ) - + # find appropriate posterior # NOTE: This is defined but not used # posterior_ids <- get.id("posteriors", "pft_id", pftid, con) - + posteriorid_query <- paste0( "SELECT id FROM posteriors WHERE pft_id=", pftid, " AND format_id=", formatid @@ -426,7 +426,7 @@ dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, ho ) posteriorid <- db.query(posteriorid_query, con)[["id"]] } - + # NOTE: Modified by Alexey Shiklomanov. # I'm not sure how this is supposed to work, but I think it's like this invisible(dbfile.insert( @@ -455,24 +455,24 @@ dbfile.posterior.insert <- function(filename, pft, mimetype, formatname, con, ho ##' } dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate pft pftid <- get.id(table = "pfts", values = "name", colnames = pft, con = con) if (is.null(pftid)) { invisible(data.frame()) } - + # find appropriate format mimetypeid <- get.id(table = "mimetypes", values = "type_string", colnames = mimetype, con = con) if (is.null(mimetypeid)) { PEcAn.logger::logger.error("mimetype ", mimetype, "does not exist") } formatid <- get.id(table = "formats", colnames = c("mimetype_id", "name"), values = c(mimetypeid, formatname), con = con) - + if (is.null(formatid)) { invisible(data.frame()) } - + # find appropriate posterior posteriorid <- db.query( query = paste0( @@ -484,7 +484,7 @@ dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PE if (is.null(posteriorid)) { invisible(data.frame()) } - + invisible(dbfile.check(type = "Posterior", container.id = posteriorid, con = con, hostname = hostname)) } @@ -509,14 +509,14 @@ dbfile.posterior.check <- function(pft, mimetype, formatname, con, hostname = PE ##' } dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + if (substr(in.path, 1, 1) != "/") { PEcAn.logger::logger.error("path to dbfiles:", in.path, " is not a valid full path") } - + # find appropriate host hostid <- get.id(table = "machines", colnames = "hostname", values = hostname, con = con, create = TRUE, dates = TRUE) - + # Query for existing dbfile record with same file_name, file_path, machine_id , # container_type, and container_id. dbfile <- invisible(db.query( @@ -528,10 +528,10 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn ), con = con )) - + if (nrow(dbfile) == 0) { # If no existing record, insert one - + insert_result <- db.query( query = paste0( "INSERT INTO dbfiles ", @@ -541,7 +541,7 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn ), con = con ) - + file.id <- insert_result[["id"]] } else if (!reuse) { # If there is an existing record but reuse==FALSE, return NA. @@ -559,7 +559,7 @@ dbfile.insert <- function(in.path, in.prefix, type, id, con, reuse = TRUE, hostn file.id <- dbfile[["id"]] } } - + # Return the new dbfile ID, or the one that existed already (reuse==T), or NA (reuse==F) return(file.id) } @@ -591,34 +591,34 @@ dbfile.check <- function(type, container.id, con, machine.check = TRUE, return.all = FALSE) { type <- match.arg(type, c("Input", "Posterior", "Model")) - + hostname <- default_hostname(hostname) - + # find appropriate host hostid <- get.id(table = "machines", colnames = "hostname", values = hostname, con = con) if (is.null(hostid)) { return(data.frame()) } - + dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::filter( .data$container_type == !!type, .data$container_id %in% !!container.id ) - + if (machine.check) { dbfiles <- dbfiles %>% dplyr::filter(.data$machine_id == !!hostid) } - + dbfiles <- dplyr::collect(dbfiles) - + if (nrow(dbfiles) > 1 && !return.all) { PEcAn.logger::logger.warn("Multiple Valid Files found on host machine. Returning last updated record.") dbfiles <- dbfiles %>% dplyr::filter(.data$updated_at == max(.data$updated_at)) } - + dbfiles } @@ -643,9 +643,9 @@ dbfile.check <- function(type, container.id, con, ##' } dbfile.file <- function(type, id, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + files <- dbfile.check(type = type, container.id = id, con = con, hostname = hostname) - + if (nrow(files) > 1) { PEcAn.logger::logger.warn("multiple files found for", id, "returned; using the first one found") invisible(file.path(files[1, "file_path"], files[1, "file_name"])) @@ -667,13 +667,13 @@ dbfile.file <- function(type, id, con, hostname = PEcAn.remote::fqdn()) { ##' } dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { hostname <- default_hostname(hostname) - + # find appropriate host hostid <- db.query(query = paste0("SELECT id FROM machines WHERE hostname='", hostname, "'"), con = con)[["id"]] if (is.null(hostid)) { invisible(NA) } - + # find file file_name <- basename(file) file_path <- dirname(file) @@ -686,7 +686,7 @@ dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { ), con = con ) - + if (nrow(ids) > 1) { PEcAn.logger::logger.warn("multiple ids found for", file, "returned; using the first one found") invisible(ids[1, "container_id"]) @@ -727,58 +727,58 @@ dbfile.id <- function(type, file, con, hostname = PEcAn.remote::fqdn()) { dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = FALSE) { - - + + # create nulls for file movement and error info error <- 0 files.sym <- 0 files.changed <- 0 files.reg <- 0 files.indb <- 0 - + # check for file type and update to make it *.file type if (file.type != "clim" | file.type != "nc") { PEcAn.logger::logger.error("File type not supported by move at this time. Currently only supports NC and CLIM files") error <- 1 } file.pattern <- paste0("*.", file.type) - - - + + + # create new directory if it doesn't exist if (!dir.exists(new.dir)) { dir.create(new.dir) } - - + + # check to make sure both directories exist if (!dir.exists(old.dir)) { PEcAn.logger::logger.error("Old File directory does not exist. Please enter valid file path") error <- 1 } - + if (!dir.exists(new.dir)) { PEcAn.logger::logger.error("New File directory does not exist. Please enter valid file path") error <- 1 } - + if (basename(new.dir) != basename(old.dir)) { PEcAn.logger::logger.error("Basenames of files do not match") } - + # list files in the old directory old.files <- list.files(path = old.dir, pattern = file.pattern) - + # check to make sure there are files if (length(old.files) == 0) { PEcAn.logger::logger.warn("No files found") error <- 1 } - + # create full file path full.old.file <- file.path(old.dir, old.files) - - + + ### Get BETY information ### con <- db.open( params = list( @@ -788,46 +788,46 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F user = "bety", password = "bety") ) - + # get matching dbfiles from BETY dbfile.path <- dirname(full.old.file) dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::collect() %>% dplyr::filter(.data$file_name %in% basename(full.old.file)) %>% dplyr::filter(.data$file_path %in% dbfile.path) - - + + # if there are matching db files if (dim(dbfiles)[1] > 0) { - + # Check to make sure files line up if (dim(dbfiles)[1] != length(full.old.file)) { PEcAn.logger::logger.warn("Files to be moved don't match up with BETY files, only moving the files that match") - + # IF DB FILES AND FULL FILES DONT MATCH, remove those not in BETY - will take care of the rest below index <- which(basename(full.old.file) %in% dbfiles$file_name) index1 <- seq(1, length(full.old.file)) check <- index1[-which(index1 %in% index)] full.old.file <- full.old.file[-check] - + # record the number of files that are being moved files.changed <- length(full.old.file) } - + # Check to make sure the files line up if (dim(dbfiles)[1] != length(full.old.file)) { PEcAn.logger::logger.error("Files to be moved don't match up with BETY files, canceling move") error <- 1 } - - + + # Make sure the files line up dbfiles <- dbfiles[order(dbfiles$file_name), ] full.old.file <- sort(full.old.file) - + # Record number of files moved and changed in BETY files.indb <- dim(dbfiles)[1] - + # Move files and update BETY if (error == 0) { for (i in 1:length(full.old.file)) { @@ -836,40 +836,40 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error if statement } # end dbfile loop - - + + # if there are files that are in the folder but not in BETY, we can either register them or not if (dim(dbfiles)[1] == 0 | files.changed > 0) { - + # Recheck what files are in the directory since others may have been moved above old.files <- list.files(path = old.dir, pattern = file.pattern) - + # Recreate full file path full.old.file <- file.path(old.dir, old.files) - - + + # Error check again to make sure there aren't any matching dbfiles dbfile.path <- dirname(full.old.file) dbfiles <- dplyr::tbl(con, "dbfiles") %>% dplyr::collect() %>% dplyr::filter(.data$file_name %in% basename(full.old.file)) %>% dplyr::filter(.data$file_path %in% dbfile.path) - + if (dim(dbfiles)[1] > 0) { PEcAn.logger::logger.error("There are still dbfiles matching these files! Canceling link or registration") error <- 1 } - - + + if (error == 0 & register == TRUE) { - + # Record how many files are being registered to BETY files.reg <- length(full.old.file) - + for (i in 1:length(full.old.file)) { file_path <- dirname(full.old.file[i]) file_name <- basename(full.old.file[i]) - + if (file.type == "nc") { mimetype <- "application/x-netcdf" formatname <- "CF Meteorology application" @@ -881,8 +881,8 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F else { PEcAn.logger::logger.error("File Type is currently not supported") } - - + + dbfile.input.insert( in.path = file_path, in.prefix = file_name, @@ -900,26 +900,26 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error loop } # end register == TRUE - + if (error == 0 & register == FALSE) { # Create file path for symbolic link full.new.file <- file.path(new.dir, old.files) - + # Record number of files that will have a symbolic link made files.sym <- length(full.new.file) - + # Line up files full.new.file <- sort(full.new.file) full.old.file <- sort(full.old.file) - + # Check to make sure the files are the same length if (length(full.new.file) != length(full.old.file)) { PEcAn.logger::logger.error("Files to be moved don't match up with BETY. Canceling Move") error <- 1 } - + # Move file and create symbolic link if there are no errors - + if (error == 0) { for (i in 1:length(full.old.file)) { fs::file_move(full.old.file[i], new.dir) @@ -927,13 +927,13 @@ dbfile.move <- function(old.dir, new.dir, file.type, siteid = NULL, register = F } # end i loop } # end error loop } # end Register == FALSE - - + + if (error > 0) { PEcAn.logger::logger.error("There was an error, files were not moved or linked") } - + if (error == 0) { PEcAn.logger::logger.info(paste0(files.changed + files.indb, " files were moved and updated on BETY, ", files.sym, " were moved and had a symbolic link created, and ", files.reg, " files were moved and then registered in BETY")) } -} # end dbfile.move() \ No newline at end of file +} # end dbfile.move() From 2cc45ec4149b87a06de18687b67c0a7cbb303c13 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Mon, 11 Mar 2024 21:47:42 -0700 Subject: [PATCH 67/68] Update models/stics/DESCRIPTION --- models/stics/DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/stics/DESCRIPTION b/models/stics/DESCRIPTION index 5e2cfbd1255..81624629cdb 100644 --- a/models/stics/DESCRIPTION +++ b/models/stics/DESCRIPTION @@ -17,7 +17,7 @@ Imports: lubridate, ncdf4, XML, - dplyr, + dplyr Suggests: SticsRFiles, testthat (>= 1.0.2) From 06c888747994c86ec02706e052af9521633ab403 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Mon, 11 Mar 2024 21:49:02 -0700 Subject: [PATCH 68/68] Update models/stics/R/write.config.STICS.R --- models/stics/R/write.config.STICS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/stics/R/write.config.STICS.R b/models/stics/R/write.config.STICS.R index c3bedf41cde..87cfbd5fa7c 100644 --- a/models/stics/R/write.config.STICS.R +++ b/models/stics/R/write.config.STICS.R @@ -48,7 +48,7 @@ write.config.STICS <- function(defaults, trait.values, settings, run.id) { ########## Determining number of USMs (could be made its own function) - # In STICS, it is 1 UMS per crop cycle, where each cycle can be 2-years max + # In STICS, it is 1 USM per crop cycle, where each cycle can be 2-years max # If we have a consecutive monoculture for > 2 years, we still need to divide it into 2-year USMs # If there are multiple pfts, this is a strong clue that there are multiple crop cycles # but it can also be the case that there is one cycle with intercropping