From b35d85d3a6eef9e01e76dabf9fdbb571da8cb311 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 10 Apr 2026 13:40:49 -0700 Subject: [PATCH 1/2] Avoid 4x slowdown from ud_convert --- models/sipnet/NEWS.md | 4 + models/sipnet/R/model2netcdf.SIPNET.R | 221 +++++++++++------- models/sipnet/man/model2netcdf.SIPNET.Rd | 5 +- .../tests/testthat/test-model2netcdf.SIPNET.R | 13 +- 4 files changed, 157 insertions(+), 86 deletions(-) diff --git a/models/sipnet/NEWS.md b/models/sipnet/NEWS.md index ef34ba30b58..7384238d53d 100644 --- a/models/sipnet/NEWS.md +++ b/models/sipnet/NEWS.md @@ -1,5 +1,9 @@ # PEcAn.SIPNET 1.10.0.9000 +* Fixed a unit error in model2netcdf.SIPNET's calculation of `AGBI` (kgC/m2/sec) + from `woodCreation` (actually gC/m2/timestep, was being treated as gC/m2/day). +* `model2netcdf.SIPNET` now takes NPP directly from sipnet.out rather than repeat + Sipnet's internal calculation (as GPP - Ra) with identical results. * Updated README with a more complete model description and instructions for installing SIPNET (#3705) * Removed `tests/Rcheck_reference.log`, which was used to ignore historic check messages that have now been fixed. * Initial support for SIPNET v2.0, whose features include simplified input files, diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R index 7cda9231a28..95031a82c87 100644 --- a/models/sipnet/R/model2netcdf.SIPNET.R +++ b/models/sipnet/R/model2netcdf.SIPNET.R @@ -39,24 +39,25 @@ mergeNC <- function( } #--------------------------------------------------------------------------------------------------# -##' Convert SIPNET output to netCDF -##' -##' Converts all output contained in a folder to netCDF. -##' -##' @param outdir Location of SIPNET model output -##' @param sitelat Latitude of the site -##' @param sitelon Longitude of the site -##' @param start_date Start time of the simulation -##' @param end_date End time of the simulation -##' @param revision model revision -##' @param overwrite Flag for overwriting nc files or not -##' @param conflict Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next. -##' @param prefix prefix to read the output files -##' @param delete.raw logical: remove sipnet.out files after converting? -##' -##' @export -##' @author Shawn Serbin, Michael Dietze -model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out", +#' Convert SIPNET output to netCDF +#' +#' Converts all output contained in a folder to netCDF. +#' +#' @param outdir Location of SIPNET model output +#' @param sitelat Latitude of the site +#' @param sitelon Longitude of the site +#' @param start_date Start time of the simulation +#' @param end_date End time of the simulation +#' @param revision model revision. +#' Ignored: PEcAn detects all relevant version differences from the format of the output file. +#' @param overwrite Flag for overwriting nc files or not +#' @param conflict Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next. +#' @param prefix prefix to read the output files +#' @param delete.raw logical: remove sipnet.out files after converting? +#' +#' @export +#' @author Shawn Serbin, Michael Dietze +model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision = NULL, prefix = "sipnet.out", overwrite = FALSE, conflict = FALSE) { ### Read in model output in SIPNET format sipnet_out_file <- file.path(outdir, prefix) @@ -118,7 +119,75 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, timestep.s <- 86400 / out_day - + + + ## Unit conversions + # + # CKB 20260407: Not using ud_convert here is intentional! + # This step is a consistent bottleneck to whole-run speed, and tests using + # ud_convert show a surprisingly large slowdown: + # In a test batch with ~500 rundirs run in parallel on a 2022-era SSD Macbook, + # the model stage took ~4.5x(!) longer with ud_convert than with simple scalars. + g_to_kg <- function(x) x / 1000 + g_step_to_kg_sec <- function(x) x / 1000 / timestep.s + cm_to_mm <- function(x) x * 10 + cm_step_to_mm_sec <- function(x) x * 10 / timestep.s + sipnet_output <- sipnet_output |> + dplyr::mutate( + + # C and N pools + dplyr::across( + .cols = c( + # C pools are mandatory + all_of(c("plantWoodC", "plantLeafC", "coarseRootC", "fineRootC", "soil", "litter")), + # N only present when turned on + any_of(c("minN", "soilOrgN", "litterN")) + ), + .fns = g_to_kg + ), + + # C and N fluxes + dplyr::across( + .cols = c( + all_of(c("gpp", "nee", "npp", "rAboveground", "rRoot", "rtot", "rSoil")), + any_of(c("woodCreation", "n2o", "nLeaching", "nFixation", "nUptake", "ch4")) + ), + .fns = g_step_to_kg_sec + ), + + # Water pools + dplyr::across( + .cols = c( + all_of(c("soilWater", "soilWetnessFrac", "snow")), + any_of("litterWater") # Only present in V1 output + ), + .fns = cm_to_mm + ), + + # Water fluxes + dplyr::across( + .cols = all_of("evapotranspiration"), + .fns = cm_step_to_mm_sec + ), + # Water flux special case: + # Sipnet reports transpiration, and no other variables, in cm/day not cm/timestep. + fluxestranspiration = cm_to_mm(.data$fluxestranspiration) / 86400, # cm/day -> mm/sec + + # Date and time + datetime = sipnet2datetime(.data$year, .data$day, .data$time) + ) + + + # calculate LAI for standard output + # LAI = plantLeafC / leafCSpWt + # both operands are in carbon units (gC/m2 and gC/m2_leaf), + # so no carbon fraction conversion (e.g. cFracLeaf) is needed. + param <- utils::read.table(file.path(gsub(pattern = "/out/", + replacement = "/run/", x = outdir), + "sipnet.param"), stringsAsFactors = FALSE) + leafCSpWt <- param[param[, 1] == "leafCSpWt", 2] + SLA <- 1000 / leafCSpWt # m2 leaf / kg C + ### Loop over years in SIPNET output to create separate netCDF outputs for (y in year_seq) { @@ -133,18 +202,12 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, file.rename(file.path(outdir, paste(y, "nc", sep = ".")), file.path(outdir, "previous.nc")) } print(paste("---- Processing year: ", y)) # turn on for debugging - + ## Subset data for processing sub.sipnet.output <- subset(sipnet_output, sipnet_output$year == y) - sub_dates <- sipnet2datetime( - y, - sub.sipnet.output[["day"]], - sub.sipnet.output[["time"]] - ) - sub_dates_cf <- PEcAn.utils::datetime2cf( - sub_dates, + sub.sipnet.output$datetime, paste0("days since ", y, "-01-01"), tz = "UTC" ) @@ -161,81 +224,71 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, ## Setup outputs for netCDF file in appropriate units output <- list( - "GPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "g/m2", "kg/m2") / timestep.s, - "NPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp - (sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot), "g/m2", "kg/m2") / timestep.s, - "TotalResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rtot, "g/m2", "kg/m2") / timestep.s, - "AutoResp" = (PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, "g/m2", "kg/m2")) / timestep.s, - "HeteroResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil - sub.sipnet.output$rRoot, "g/m2", "kg/m2") / timestep.s, - "SoilResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil, "g/m2", "kg/m2") / timestep.s, - "NEE" = PEcAn.utils::ud_convert(sub.sipnet.output$nee, "g/m2", "kg/m2") / timestep.s, - "AbvGrndWood" = PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC, "g/m2", "kg/m2"), - "leaf_carbon_content" = PEcAn.utils::ud_convert(sub.sipnet.output$plantLeafC, "g/m2", "kg/m2"), - "TotLivBiom" = (PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC + - sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC, "g/m2", "kg/m2")), - "TotSoilCarb" = PEcAn.utils::ud_convert(sub.sipnet.output$soil + sub.sipnet.output$litter, "g/m2", "kg/m2") + "GPP" = sub.sipnet.output$gpp, + "NPP" = sub.sipnet.output$npp, + "TotalResp" = sub.sipnet.output$rtot, + "AutoResp" = sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, + "HeteroResp" = sub.sipnet.output$rSoil - sub.sipnet.output$rRoot, + "SoilResp" = sub.sipnet.output$rSoil, + "NEE" = sub.sipnet.output$nee, + "AbvGrndWood" = sub.sipnet.output$plantWoodC, + "leaf_carbon_content" = sub.sipnet.output$plantLeafC, + "litter_carbon_content" = sub.sipnet.output$litter, + "fine_root_carbon_content" = sub.sipnet.output$fineRootC, + "coarse_root_carbon_content" = sub.sipnet.output$coarseRootC, + "LAI" = sub.sipnet.output$plantLeafC * SLA, + "TotLivBiom" = sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC + + sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC, + "TotSoilCarb" = sub.sipnet.output$soil + sub.sipnet.output$litter, + "AGB" = sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, + + # Water variables: + # Liquid water units are cm in Sipnet; in PEcAn they're kg water m-2 + # (which is equivalent to mm: (water density = 1000 kg m-3) * (1 m/ 1000 mm) = (1 kg m-2)/mm + # Evapotranspiration in SIPNET is cm^3 water per cm^2 of area, + # already converted above to mm sec-1. + # To convert it to latent heat units W/m2 multiply by latent heat of vaporization (J kg-1) + # Latent heat of vaporization is not constant and it varies slightly with temperature, get.lv() returns 2.5e6 J kg-1 by default + "Qle" = sub.sipnet.output$evapotranspiration * PEcAn.data.atmosphere::get.lv(), # Qle W/m2/sec + "Transp" = sub.sipnet.output$fluxestranspiration, + "SoilMoist" = sub.sipnet.output$soilWater, + "SoilMoistFrac" = sub.sipnet.output$soilWetnessFrac, + "SWE" = sub.sipnet.output$snow # Snow Water Equivalent ) - # Water variables: - # Liquid water units are cm in Sipnet; in PEcAn they're kg water m-2 - # (which is equivalent to mm, but ud_convert doesn't know that) - # Evapotranspiration in SIPNET is cm^3 water per cm^2 of area, to convert it to latent heat units W/m2 multiply with : - # 0.01 (cm2m) * 1000 (water density, kg m-3) * latent heat of vaporization (J kg-1) - # Latent heat of vaporization is not constant and it varies slightly with temperature, get.lv() returns 2.5e6 J kg-1 by default - output[["Qle"]] <- (PEcAn.utils::ud_convert(sub.sipnet.output$evapotranspiration, "cm", "mm") * PEcAn.data.atmosphere::get.lv()) / timestep.s # Qle W/m2 - # Note that Sipnet reports transpiration, and no other variables, in cm/day not cm/timestep. - output[["Transp"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$fluxestranspiration, "cm/day", "mm/sec") # Transpiration - output[["SoilMoist"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$soilWater, "cm", "mm") # Soil moisture kgW/m2 - output[["SoilMoistFrac"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$soilWetnessFrac, "cm", "mm") # Fractional soil wetness - output[["SWE"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$snow, "cm", "mm") # Snow Water Equivalent - output[["litter_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litter, "g/m2", "kg/m2") - # litterWater was removed in SIPNET v2; only extract if present - if ("litterWater" %in% names(sub.sipnet.output)) { - output[["litter_mass_content_of_water"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litterWater, "cm", "mm") + if ("litterWater" %in% names(sub.sipnet.output)) { # Removed in SIPNET v2; only extract if present + output[["litter_mass_content_of_water"]] <- sub.sipnet.output$litterWater } - - #calculate LAI for standard output - # LAI = plantLeafC / leafCSpWt - # both operands are in carbon units (gC/m2 and gC/m2_leaf), - # so no carbon fraction conversion (e.g. cFracLeaf) is needed. - param <- utils::read.table(file.path(gsub(pattern = "/out/", - replacement = "/run/", x = outdir), - "sipnet.param"), stringsAsFactors = FALSE) - leafCSpWt <- param[param[, 1] == "leafCSpWt", 2] - SLA <- 1000 / leafCSpWt # m2 leaf / kg C - output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA - output[["fine_root_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$fineRootC, "g/m2", "kg/m2") - output[["coarse_root_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$coarseRootC, "g/m2", "kg/m2") - if ("woodCreation" %in% names(sub.sipnet.output)) { - output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "g/m2/day", "kg/m2/s") + if ("woodCreation" %in% names(sub.sipnet.output)) { # Added in SIPNET v2; only extract if present + output[["GWBI"]] <- sub.sipnet.output$woodCreation } - output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "g/m2", "kg/m2") + # columns only present in sipnet >= v2 with N and methane turned on if ("minN" %in% names(sub.sipnet.output)) { - output[["mineral_N"]] <- sub.sipnet.output$minN * 0.001 # gN/m2 -> kgN/m2 + output[["mineral_N"]] <- sub.sipnet.output$minN } if ("soilOrgN" %in% names(sub.sipnet.output)) { - output[["soil_organic_N"]] <- sub.sipnet.output$soilOrgN * 0.001 # kgN/m2 + output[["soil_organic_N"]] <- sub.sipnet.output$soilOrgN } if ("litterN" %in% names(sub.sipnet.output)) { - output[["litter_N"]] <- sub.sipnet.output$litterN * 0.001 # kgN/m2 + output[["litter_N"]] <- sub.sipnet.output$litterN } if ("n2o" %in% names(sub.sipnet.output)) { - output[["N2O_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$n2o, "g/m2", "kg/m2") / timestep.s - # convert g N m-2 per timestep -> kg N m-2 s-1 + output[["N2O_flux"]] <- sub.sipnet.output$n2o } if ("nLeaching" %in% names(sub.sipnet.output)) { - output[["N_leaching"]] <- (sub.sipnet.output$nLeaching * 0.001) / timestep.s # kgN/m2/s + output[["N_leaching"]] <- sub.sipnet.output$nLeaching } if ("nFixation" %in% names(sub.sipnet.output)) { - output[["N_fixation"]] <- (sub.sipnet.output$nFixation * 0.001) / timestep.s # kgN/m2/s + output[["N_fixation"]] <- sub.sipnet.output$nFixation } if ("nUptake" %in% names(sub.sipnet.output)) { - output[["N_uptake"]] <- (sub.sipnet.output$nUptake * 0.001) / timestep.s # kgN/m2/s + output[["N_uptake"]] <- sub.sipnet.output$nUptake } if ("ch4" %in% names(sub.sipnet.output)) { - output[["CH4_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$ch4, "g/m2", "kg/m2") / timestep.s - # convert g C m-2 per timestep -> kg C m-2 s-1 + output[["CH4_flux"]] <- sub.sipnet.output$ch4 } + output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2])) # ******************** Declare netCDF variables ********************# @@ -376,10 +429,12 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, sipnet2datetime <- function(year, doy, hour){ hr <- floor(hour) - minsec <- PEcAn.utils::ud_convert(hour - hr, "hour", "min") + # minsec <- PEcAn.utils::ud_convert(hour - hr, "hour", "min") + minsec <- (hour - hr) * 60 minute <- floor(minsec) - sec <- PEcAn.utils::ud_convert(minsec - minute, "minute", "second") + # sec <- PEcAn.utils::ud_convert(minsec - minute, "minute", "second") + sec <- (minsec - minute) * 60 minute <- ifelse(sec == 60, minute + 1, minute) sec <- ifelse(sec == 60, 0, sec) diff --git a/models/sipnet/man/model2netcdf.SIPNET.Rd b/models/sipnet/man/model2netcdf.SIPNET.Rd index b88a96198cf..5042dbc7343 100644 --- a/models/sipnet/man/model2netcdf.SIPNET.Rd +++ b/models/sipnet/man/model2netcdf.SIPNET.Rd @@ -11,7 +11,7 @@ model2netcdf.SIPNET( start_date, end_date, delete.raw = FALSE, - revision, + revision = NULL, prefix = "sipnet.out", overwrite = FALSE, conflict = FALSE @@ -30,7 +30,8 @@ model2netcdf.SIPNET( \item{delete.raw}{logical: remove sipnet.out files after converting?} -\item{revision}{model revision} +\item{revision}{model revision. +Ignored: PEcAn detects all relevant version differences from the format of the output file.} \item{prefix}{prefix to read the output files} diff --git a/models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R b/models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R index 849d250ef04..6af515f42e8 100644 --- a/models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R +++ b/models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R @@ -44,7 +44,7 @@ make_base_sipnet <- function(n = 4L) { year = 2002, day = rep(c(1, 2), each = n / 2, length.out = n), time = rep(c(6, 18), length.out = n), - plantWoodC = 5000, plantLeafC = 200, woodCreation = 0.5, + plantWoodC = 5000, plantLeafC = 200, soil = 10000, microbeC = 8, coarseRootC = 1200, fineRootC = 800, litter = 400, soilWater = 14, soilWetnessFrac = 0.85, snow = 0, npp = 0.05, nee = 0.10, cumNEE = cumsum(rep(0.1, n)), @@ -226,6 +226,17 @@ test_that("fluxes are converted from gC/m2/timestep to kg/m2/sec", { expect_equal(pec$GPP, sip$gpp / 1000 / ts) expect_equal(pec$Transp, sip$fluxestranspiration * 10 / ts, tolerance = 1e-6) + + sip2 <- make_v2_sipnet() + out2 <- out_dir <- setup_sipnet_test(sip2)$outdir + pec2 <- PEcAn.utils::read.output( + ncfiles = file.path(out2, "2002.nc"), + variables = c("GPP", "GWBI", "Transp"), + dataframe = TRUE, + verbose = FALSE, + print_summary = FALSE + ) + expect_equal(pec2$GWBI, sip2$woodCreation / 1000 / ts, tolerance = 1e-6) }) From 89577f1600569f4ec9ec6e0815465da4a1ca58a3 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 10 Apr 2026 19:37:18 -0700 Subject: [PATCH 2/2] namespace --- models/sipnet/R/model2netcdf.SIPNET.R | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R index 95031a82c87..5d96f615fee 100644 --- a/models/sipnet/R/model2netcdf.SIPNET.R +++ b/models/sipnet/R/model2netcdf.SIPNET.R @@ -139,9 +139,9 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, dplyr::across( .cols = c( # C pools are mandatory - all_of(c("plantWoodC", "plantLeafC", "coarseRootC", "fineRootC", "soil", "litter")), + dplyr::all_of(c("plantWoodC", "plantLeafC", "coarseRootC", "fineRootC", "soil", "litter")), # N only present when turned on - any_of(c("minN", "soilOrgN", "litterN")) + dplyr::any_of(c("minN", "soilOrgN", "litterN")) ), .fns = g_to_kg ), @@ -149,8 +149,8 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, # C and N fluxes dplyr::across( .cols = c( - all_of(c("gpp", "nee", "npp", "rAboveground", "rRoot", "rtot", "rSoil")), - any_of(c("woodCreation", "n2o", "nLeaching", "nFixation", "nUptake", "ch4")) + dplyr::all_of(c("gpp", "nee", "npp", "rAboveground", "rRoot", "rtot", "rSoil")), + dplyr::any_of(c("woodCreation", "n2o", "nLeaching", "nFixation", "nUptake", "ch4")) ), .fns = g_step_to_kg_sec ), @@ -158,15 +158,15 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, # Water pools dplyr::across( .cols = c( - all_of(c("soilWater", "soilWetnessFrac", "snow")), - any_of("litterWater") # Only present in V1 output + dplyr::all_of(c("soilWater", "soilWetnessFrac", "snow")), + dplyr::any_of("litterWater") # Only present in V1 output ), .fns = cm_to_mm ), # Water fluxes dplyr::across( - .cols = all_of("evapotranspiration"), + .cols = dplyr::all_of("evapotranspiration"), .fns = cm_step_to_mm_sec ), # Water flux special case: