diff --git a/DESCRIPTION b/DESCRIPTION index fba3448..e26e721 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: jNSMR Title: Wrapper for Java Newhall Simulation Model (jNSM) "A Traditional Soil Climate Simulation Model" -Version: 0.2.0 +Version: 0.2.0.9001 Authors@R: person(given = "Andrew", family = "Brown", diff --git a/LICENSE b/LICENSE index 246a979..8e2de53 100644 --- a/LICENSE +++ b/LICENSE @@ -1,11 +1,4 @@ -Newhall 1.6.1, Copyright (C) 2010-2011 -United States Department of Agriculture - Natural Resources Conservation Service, -Penn State University Center for Environmental Informatics -All rights reserved. - -This product includes software developed by the JDOM Project (http://www.jdom.org/). - -JDOM 1.1, Copyright (C) 2000-2004 Jason Hunter & Brett McLaughlin +Copyright 2010-2023 United States Department of Agriculture - Natural Resources Conservation Service, Soil and Plant Science Division Staff, Penn State University Center for Environmental Informatics All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: diff --git a/R/AAAA.R b/R/AAAA.R index 31a0d97..b99545d 100644 --- a/R/AAAA.R +++ b/R/AAAA.R @@ -53,4 +53,4 @@ newhall_version <- function() { )) } -.get_default_suffix <- function() "-1.6.4-3" \ No newline at end of file +.get_default_suffix <- function() "-1.6.5" \ No newline at end of file diff --git a/R/NewhallBatch.R b/R/NewhallBatch.R index 15ff460..80a3ba4 100644 --- a/R/NewhallBatch.R +++ b/R/NewhallBatch.R @@ -7,8 +7,10 @@ #' #' @param .data a _data.frame_ or _character_ vector of paths to CSV files; or a SpatRaster or RasterStack containin the same data elements and names as included in the batch `data.frame`/CSV format #' @param unitSystem Default: `"metric"` OR `"mm"` OR `"cm"` use _millimeters_ of rainfall (default for the BASIC model); set to `unitSystem="english"` OR `unitSystem="in"` to transform English (inches of precipitation; degrees Fahrenheit) inputs to metric (millimeters of precipitation; degrees Celsius) before running simulation -#' @param soilAirOffset air-soil temperature offset. Conventionally for jNSM: `2.5` for metric units (default); `4.5` for english units. -#' @param amplitude difference in amplitude between soil and air temperature sine waves. Default `0.66` +#' @param soilAirOffset air-soil temperature offset. Conventionally for jNSM: `2.5` for metric units (default); `4.5` for english units. Can optionally be specified as a layer in a raster input. +#' @param amplitude difference in amplitude between soil and air temperature sine waves. Default `0.66`. Can optionally be specified as a layer in a raster input. +#' @param hasOHorizon Used for cryic soil temperature regime criteria. Default: `FALSE`. Can optionally be specified as a layer in a raster input. +#' @param isSaturated Used for cryic soil temperature regime and aquic soil moisture regime mask. Default: `FALSE`. Can optionally be specified as a layer in a raster input. #' @param verbose print message about number of simulations and elapsed time #' @param toString call `toString()` method on each _NewhallResults_ object and store in `output` column of result? #' @param checkargs _logical_; check argument length and data types for each run? Default: `TRUE` @@ -84,6 +86,8 @@ newhall_batch.default <- function(.data = NULL, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in","english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -94,7 +98,18 @@ newhall_batch.default <- function(.data = NULL, overwrite = NULL) { # if newer JAR is available, use the fastest batch method - if (newhall_version() >= "1.6.3") { + if (newhall_version() >= "1.6.5") { + batch3( + .data, + unitSystem = unitSystem, + soilAirOffset = soilAirOffset, + amplitude = amplitude, + verbose = verbose, + toString = toString, + checkargs = checkargs + ) + } + if (newhall_version() >= "1.6.3" && newhall_version() < "1.6.5") { batch2( .data, unitSystem = unitSystem, @@ -330,7 +345,9 @@ batch2 <- function(.data, rJava::.jarray(rep(unitSystem == "cm", nrow(.data))), rJava::.jarray(as.double(.data$awc)), rJava::.jarray(rep(soilAirOffset, nrow(.data))), - rJava::.jarray(rep(amplitude, nrow(.data))) + rJava::.jarray(rep(amplitude, nrow(.data))), + rJava::.jarray(rep(FALSE, nrow(.data))), # O horizon + rJava::.jarray(rep(FALSE, nrow(.data))) # saturation ) b <- rJava::.jcast(res, "Lorg/psu/newhall/sim/NewhallBatchResults") @@ -374,6 +391,144 @@ batch2 <- function(.data, type.convert(res, as.is = TRUE) } +# data.frame -> data.frame +#' +#' @importFrom utils type.convert +batch3 <- function(.data, + unitSystem = "metric", + soilAirOffset = ifelse(unitSystem %in% + c("in", "english"), + 4.5, 2.5), + amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, + verbose = TRUE, + toString = TRUE, + checkargs = TRUE) { + + t1 <- Sys.time() + x <- BASICSimulationModel() + + unitSystem <- match.arg(tolower(unitSystem), + choices = c("metric","mm","cm","in","english")) + + if (unitSystem %in% c("metric","mm")) { + # "cm" is the internal convention in the NewhallDatasetMetadata for _millimeters_ of rainfall, degrees Celsius + unitSystem <- "cm" + } else if (unitSystem == "english") { + # "in" is the internal convention in the NewhallDatasetMetadata for inches of rainfall, degrees Fahrenheit + unitSystem <- "in" + } + + # minimum dataset includes all of the codes specified in colnames of batch file template + mincols <- !(.colnamesNewhallBatch() %in% colnames(.data)) + + if (sum(mincols) > 0) { + stop(sprintf( + "columns %s are required in the Newhall batch CSV input format", + paste0(.colnamesNewhallBatch()[mincols], collapse = ", ") + ), call. = FALSE) + } + + # convert deg F to deg C + .doUnitsTemp <- function(x) if (unitSystem == "in") return((x - 32) * 5 / 9) else x + + # convert inches to _millimeters_ + .doUnitsLength <- function(x) if (unitSystem == "in") return(x * 25.4) else x + + .data <- data.frame(.data) + cnd <- colnames(.data) + .SD <- NULL + + # calculate hemisphere + hem <- rep(rJava::.jchar(strtoi(charToRaw('N'), 16L)), nrow(.data)) + hem[.data$latDD < 0] <- rJava::.jchar(strtoi(charToRaw('S'), 16L)) + + # handle arguments or gridded inputs for + # - soil-air temperature offset + # - soil temperature amplitude + # - O horizon presence/absence + # - saturation (50cm depth) presence/absence + + sao <- as.double(.data$soilAirOffset) + if (length(sao) == 0) { + sao <- rep(soilAirOffset, nrow(.data)) + } + + amp <- as.double(.data$amplitude) + if (length(amp) == 0) { + amp <- rep(amplitude, nrow(.data)) + } + + ohz <- as.logical(.data$hasOHorizon) + if (length(ohz) == 0) { + ohz <- rep(hasOHorizon, nrow(.data)) + } + + isa <- as.logical(.data$isSaturated) + if (length(isa) == 0) { + isa <- rep(isSaturated, nrow(.data)) + } + + res <- rJava::.jcall(x, "Lorg/psu/newhall/sim/NewhallBatchResults;", "runBatch2", + # res <- x$runBatch2( + rJava::.jarray(cbind(0.0, as.matrix(.data[, cnd[grep("^p[A-Z][a-z]{2}$", cnd)]])), dispatch = TRUE), + rJava::.jarray(cbind(0.0, as.matrix(.data[, cnd[grep("^t[A-Z][a-z]{2}$", cnd)]])), dispatch = TRUE), + rJava::.jarray(as.double(.data$latDD)), + rJava::.jarray(as.double(.data$lonDD)), + rJava::.jarray(rJava::.jchar(hem)), + rJava::.jarray(as.double(.data$elev)), + # rJava::.jarray(rJava::.jchar(c(rJava::.jchar(strtoi(charToRaw('N'), 16L)), rJava::.jchar(strtoi(charToRaw('S'), 16L))))[as.integer(.data$latDD < 0) + 1]), + # rJava::.jarray(rJava::.jchar(c(rJava::.jchar(strtoi(charToRaw('E'), 16L)), rJava::.jchar(strtoi(charToRaw('W'), 16L))))[as.integer(.data$latDD > 0) + 1]), + rJava::.jarray(rep(unitSystem == "cm", nrow(.data))), + rJava::.jarray(as.double(.data$awc)), + rJava::.jarray(sao), + rJava::.jarray(amp), + rJava::.jarray(ohz), + rJava::.jarray(isa) + ) + + b <- rJava::.jcast(res, "Lorg/psu/newhall/sim/NewhallBatchResults") + + # store arrays of values as public fields of NewhallBatchResults + fields <- c( + "annualRainfall", + "waterHoldingCapacity", + "annualWaterBalance", + "annualPotentialEvapotranspiration", + "summerWaterBalance", + "dryDaysAfterSummerSolstice", + "moistDaysAfterWinterSolstice", + "numCumulativeDaysDry", + "numCumulativeDaysMoistDry", + "numCumulativeDaysMoist", + "numCumulativeDaysDryOver5C", + "numCumulativeDaysMoistDryOver5C", + "numCumulativeDaysMoistOver5C", + "numConsecutiveDaysMoistInSomeParts", + "numConsecutiveDaysMoistInSomePartsOver8C", + "temperatureRegime", + "moistureRegime", + "regimeSubdivision1", + "regimeSubdivision2" + ) + fieldsmatrix <- c("meanPotentialEvapotranspiration","temperatureCalendar", "moistureCalendar") + + # convert to data frame + res <- lapply(fields, function(n) rJava::.jfield(b, name = n)) + res <- lapply(res, function(x) {if (length(x) == 0) return(rep(NA, length(res[[1]]))); x}) + res <- as.data.frame(res) + if (verbose) { + deltat <- signif(difftime(Sys.time(), t1, units = "auto"), digits = 2) + message(sprintf( + "newhall_batch: ran n=%s simulations in %s %s", + nrow(res), deltat, attr(deltat, 'units') + )) + } + colnames(res) <- fields #, fieldsmatrix + type.convert(res, as.is = TRUE) +} + #' @export #' @examples #' @@ -426,6 +581,8 @@ newhall_batch <- function(.data, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -443,6 +600,8 @@ newhall_batch.character <- function(.data, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -512,6 +671,8 @@ newhall_batch.SpatRaster <- function(.data, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in","english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = FALSE, checkargs = TRUE, @@ -727,6 +888,8 @@ newhall_batch.RasterBrick <- function(.data, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in","english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -758,6 +921,8 @@ newhall_batch.RasterStack <- function(.data, unitSystem = "metric", soilAirOffset = ifelse(unitSystem %in% c("in","english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, diff --git a/README.Rmd b/README.Rmd index 8bf5429..0ad4a7c 100644 --- a/README.Rmd +++ b/README.Rmd @@ -79,12 +79,7 @@ head(res) ## License information -**This package uses a modified version of the Newhall model v1.6.1 (released 2016/02/10) of the jNSM (official download here: https://www.nrcs.usda.gov/wps/portal/nrcs/detail/?cid=nrcs142p2_053559)**. The compiled JAR and source code are distributed in this R package under the "New" (3-Clause) BSD License. See _LICENSE_ for more information. Modifications to the JAR relative to legacy version are only to facilitate higher throughput batching and access to additional data elements. - -> Newhall 1.6.1, Copyright (C) 2010-2011 -> United States Department of Agriculture - Natural Resources Conservation Service, -> Penn State University Center for Environmental Informatics -> All rights reserved. +**This package uses a modified version of the Newhall model v1.6.1 (released 2016/02/10) of the jNSM (official download here: https://www.nrcs.usda.gov/resources/education-and-teaching-materials/java-newhall-simulation-model-jnsm)**. The compiled JAR and source code are distributed in this R package under the "New" (3-Clause) BSD License. See _LICENSE_ for more information. Modifications to the JAR relative to legacy version facilitate higher throughput and access to additional data elements. ## System requirements diff --git a/README.md b/README.md index 40492c7..ac7fb25 100644 --- a/README.md +++ b/README.md @@ -36,8 +36,8 @@ paths, a data.frame, a SpatRaster or RasterStack/Brick object as input. ### GeoTIFF/SpatRaster Input library(jNSMR) - #> jNSMR (0.1.2.9000) -- R interface to the classic Newhall Simulation Model - #> Added JAR file (newhall-1.6.4.jar) to Java class path. + #> jNSMR (0.1.2.9001) -- R interface to the classic Newhall Simulation Model + #> Added JAR file (newhall-1.6.5.jar) to Java class path. library(terra) #> terra 1.7.39 @@ -46,7 +46,7 @@ paths, a data.frame, a SpatRaster or RasterStack/Brick object as input. x$elev <- 0 # elevation is not currently used by the model directly y <- newhall_batch(x) ## full resolution - #> newhall_batch: ran n=18790 simulations in 24 secs + #> newhall_batch: ran n=18790 simulations in 26 secs par(mfrow = c(2, 1)) terra::plot(y$annualWaterBalance, main = "Annual Water Balance (P-PET)") @@ -136,16 +136,11 @@ directory of this package. **This package uses a modified version of the Newhall model v1.6.1 (released 2016/02/10) of the jNSM (official download here: -)**. +)**. The compiled JAR and source code are distributed in this R package under the “New” (3-Clause) BSD License. See *LICENSE* for more information. -Modifications to the JAR relative to legacy version are only to -facilitate higher throughput batching and access to additional data -elements. - -> Newhall 1.6.1, Copyright (C) 2010-2011 United States Department of -> Agriculture - Natural Resources Conservation Service, Penn State -> University Center for Environmental Informatics All rights reserved. +Modifications to the JAR relative to legacy version facilitate higher +throughput and access to additional data elements. ## System requirements diff --git a/inst/java/newhall-1.6.5.jar b/inst/java/newhall-1.6.5.jar new file mode 100644 index 0000000..bce49d1 Binary files /dev/null and b/inst/java/newhall-1.6.5.jar differ diff --git a/java/org/psu/newhall/sim/BASICSimulationModel.java b/java/org/psu/newhall/sim/BASICSimulationModel.java index 7205a5a..2e84809 100644 --- a/java/org/psu/newhall/sim/BASICSimulationModel.java +++ b/java/org/psu/newhall/sim/BASICSimulationModel.java @@ -31,13 +31,14 @@ public static NewhallBatchResults runBatch(String[] name, String[] country, doub } return new NewhallBatchResults(r); } - public static NewhallBatchResults runBatch2(double[][] precipitation, double[][] temperature, double latitude[], double longitude[], char nsHemisphere[], double elevation[], boolean[] isMetric, double[] waterholdingCapacity, double[] fc, double[] fcd) { + public static NewhallBatchResults runBatch2(double[][] precipitation, double[][] temperature, double latitude[], double longitude[], char nsHemisphere[], double elevation[], boolean[] isMetric, double[] waterholdingCapacity, double[] fc, double[] fcd, boolean[] hasOHorizon, boolean[] isSaturated) { NewhallResults[] r = new NewhallResults[latitude.length]; for (int i = 0; i <= (latitude.length - 1); i++) { - r[i] = runSimulation(temperature[i], precipitation[i], latitude[i], longitude[i], elevation[i], nsHemisphere[i], isMetric[i], waterholdingCapacity[i], fc[i], fcd[i]); + r[i] = runSimulation(temperature[i], precipitation[i], latitude[i], longitude[i], elevation[i], nsHemisphere[i], isMetric[i], waterholdingCapacity[i], fc[i], fcd[i], hasOHorizon[i], isSaturated[i]); } return new NewhallBatchResults(r); } + public static NewhallResults runSimulation(double[] temperature, double[] precip, double latDD, @@ -47,7 +48,9 @@ public static NewhallResults runSimulation(double[] temperature, boolean isMetric, double waterHoldingCapacity, double fc, - double fcd) { + double fcd, + boolean hasOHorizon, + boolean isSaturated) { // This is a hack to prevent a problematic region of code from running // multiple times. Problem datasets that test that this hack works @@ -230,13 +233,25 @@ public static NewhallResults runSimulation(double[] temperature, * that is true indicates the temp regime. */ + // properly handle cryic criteria + int cryic_ht = 15; + if (hasOHorizon) { + if(isSaturated) { + cryic_ht = 6; + } else { + cryic_ht = 8; + } + } else if (isSaturated) { + cryic_ht = 13; + } + boolean[] cr = new boolean[13]; boolean[] reg = new boolean[13]; cr[1] = tma < 0; // Mean annual air temp (MAAT) < 0C. cr[2] = 0 <= tma && tma < 8; // 0C <= MAAT <= 8C. // cr[3] = (st - cs) < 15; // Summer temp ave minus (summer/winter difference * (1 - SOIL_AIR_REL) * 0.5) < 15C. // TODO: where did latter part ^^ (... - SWD...) of this come from? Misread cryic crit 1? - cr[3] = st >= 0 && st < 15; // "non-saturated, organic surface, mean _summer_ soil temperature between 0 and 8C/15C + cr[3] = st >= 0 && st < cryic_ht; // "non-saturated, organic surface, mean _summer_ soil temperature between 0 and 8C/15C // TODO: st upper limit depends on saturation, O horizon cr[7] = (dif * fcd) >= 6; // Summer/winter difference * SOIL_AIR_REL >= 6 // NOTE: Taxonomy clearly states difference greater/equal than 6, not 5 diff --git a/man/figures/README-spatraster-ex-2.png b/man/figures/README-spatraster-ex-2.png index 9efd1cf..e1bbdb4 100644 Binary files a/man/figures/README-spatraster-ex-2.png and b/man/figures/README-spatraster-ex-2.png differ diff --git a/man/newhall_batch.Rd b/man/newhall_batch.Rd index 1c8b3b8..44abcfb 100644 --- a/man/newhall_batch.Rd +++ b/man/newhall_batch.Rd @@ -14,6 +14,8 @@ unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -29,6 +31,8 @@ newhall_batch( unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -44,6 +48,8 @@ newhall_batch( unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -59,6 +65,8 @@ newhall_batch( unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = FALSE, checkargs = TRUE, @@ -75,6 +83,8 @@ newhall_batch( unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -91,6 +101,8 @@ newhall_batch( unitSystem = "metric", soilAirOffset = ifelse(unitSystem \%in\% c("in", "english"), 4.5, 2.5), amplitude = 0.66, + hasOHorizon = FALSE, + isSaturated = FALSE, verbose = TRUE, toString = TRUE, checkargs = TRUE, @@ -107,9 +119,13 @@ newhall_batch( \item{unitSystem}{Default: \code{"metric"} OR \code{"mm"} OR \code{"cm"} use \emph{millimeters} of rainfall (default for the BASIC model); set to \code{unitSystem="english"} OR \code{unitSystem="in"} to transform English (inches of precipitation; degrees Fahrenheit) inputs to metric (millimeters of precipitation; degrees Celsius) before running simulation} -\item{soilAirOffset}{air-soil temperature offset. Conventionally for jNSM: \code{2.5} for metric units (default); \code{4.5} for english units.} +\item{soilAirOffset}{air-soil temperature offset. Conventionally for jNSM: \code{2.5} for metric units (default); \code{4.5} for english units. Can optionally be specified as a layer in a raster input.} -\item{amplitude}{difference in amplitude between soil and air temperature sine waves. Default \code{0.66}} +\item{amplitude}{difference in amplitude between soil and air temperature sine waves. Default \code{0.66}. Can optionally be specified as a layer in a raster input.} + +\item{hasOHorizon}{Used for cryic soil temperature regime criteria. Default: \code{FALSE}. Can optionally be specified as a layer in a raster input.} + +\item{isSaturated}{Used for cryic soil temperature regime and aquic soil moisture regime mask. Default: \code{FALSE}. Can optionally be specified as a layer in a raster input.} \item{verbose}{print message about number of simulations and elapsed time}