diff --git a/DESCRIPTION b/DESCRIPTION index 4733adc..f4882dd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: rosettaPTF Title: R Frontend for Rosetta Pedotransfer Functions -Version: 0.1.3 +Version: 0.1.4 Author: Soil and Plant Science Division Staff Maintainer: Andrew G. Brown Description: Access Python rosetta-soil pedotransfer functions in an R environment. Rosetta is a neural network-based model for predicting unsaturated soil hydraulic parameters from basic soil characterization data. The model predicts parameters for the van Genuchten unsaturated soil hydraulic properties model, using sand, silt, and clay, bulk density and water content. The codebase is now maintained by Dr. Todd Skaggs and other U.S. Department of Agriculture employees. This R package is intended to provide for use cases that involve many thousands of calls to the pedotransfer function. Less demanding use cases are encouraged to use the web interface or API endpoint. There are additional wrappers of the API endpoints provided by the soilDB R package `ROSETTA()` method. diff --git a/R/find_python.R b/R/find_python.R index f021fe1..62ed4a5 100644 --- a/R/find_python.R +++ b/R/find_python.R @@ -69,7 +69,7 @@ find_python <- function(envname = NULL, if (length(pypath) == 0) { # find newest python installation with rosetta installed - x <- try(reticulate::py_discover_config("rosetta")) + x <- try(reticulate::py_discover_config("rosetta"), silent = TRUE) if (!is.null(x$python_versions)) { xxx <- lapply(x$python_versions, function(x) { y <- gsub("Python ", "", system(paste(shQuote(x), "--version"), intern = TRUE, ignore.stdout = TRUE, ignore.stderr = TRUE)) @@ -82,7 +82,7 @@ find_python <- function(envname = NULL, # otherwise find newest python installation if (is.null(res) || inherits(res, 'try-error')) { - x <- try(reticulate::py_discover_config()) + x <- try(reticulate::py_discover_config(), silent = TRUE) if (!is.null(x$python_versions)) { xxx <- lapply(x$python_versions, function(x) { y <- gsub("Python ", "", system(paste(shQuote(x), "--version"), intern = TRUE, ignore.stdout = TRUE, ignore.stderr = TRUE, show.output.on.console = FALSE)) diff --git a/R/install_rosetta.R b/R/install_rosetta.R index fe8d107..bcd92c2 100644 --- a/R/install_rosetta.R +++ b/R/install_rosetta.R @@ -20,12 +20,13 @@ install_rosetta <- function(envname = NULL, arcpy_path = getOption("rosettaPTF.arcpy_path")) { # use heuristics to find python executable - if (!is.null(arcpy_path) && dir.exists(arcpy_path)){ + if (!is.null(arcpy_path) && dir.exists(arcpy_path)) { find_python(envname = envname, arcpy_path = arcpy_path) } # get rosetta-soil (and numpy if needed) - try(reticulate::py_install("rosetta-soil", envname = envname, method = method, conda = conda, pip = pip)) + try(reticulate::py_install("rosetta-soil", envname = envname, method = method, conda = conda, pip = pip), + silent = TRUE) # load modules globally in package (prevents having to reload rosettaPTF library in session) .loadModules() diff --git a/R/run_rosetta.R b/R/run_rosetta.R index 885204f..6a87a2a 100644 --- a/R/run_rosetta.R +++ b/R/run_rosetta.R @@ -54,6 +54,7 @@ run_rosetta <- function(soildata, vars = NULL, rosetta_version = 3, cores = 1, + core_thresh = NULL, file = NULL, nrows = NULL, overwrite = NULL) @@ -122,8 +123,9 @@ run_rosetta.RasterStack <- function(soildata, vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { ## for in memory only, can just convert to data.frame and use that method # res <- run_rosetta(raster::as.data.frame(soildata), @@ -153,8 +155,9 @@ run_rosetta.RasterBrick <- function(soildata, vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { run_rosetta(terra::rast(soildata), vars = vars, @@ -166,8 +169,9 @@ run_rosetta.RasterBrick <- function(soildata, ) } #' @param cores number of cores; used only for processing _SpatRaster_ or _Raster*_ input +#' @param core_thresh Magic number for determining processing chunk size. Default `20000L`. Used to calculate default `nrows` #' @param file path to write incremental raster processing output for large inputs that do not fit in memory; passed to `terra::writeStart()` and used only for processing _SpatRaster_ or _Raster*_ input; defaults to a temporary file created by `tempfile()` if needed -#' @param nrows number of rows to use per block; passed to `terra::readValues()` `terra::writeValues()`; used only for processing _SpatRaster_ or _Raster*_ input; defaults to number of rows in dataset if needed +#' @param nrows number of rows to use per block chunk; passed to `terra::readValues()` and `terra::writeValues()`; used only for processing _SpatRaster_ or _Raster*_ inputs. Defaults to the total number of rows divided by the number of cells divided by `core_thresh`. #' @param overwrite logical; overwrite `file`? passed to `terra::writeStart()`; defaults to `TRUE` if needed #' @export #' @rdname run_rosetta @@ -177,10 +181,19 @@ run_rosetta.SpatRaster <- function(soildata, vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { - terra::readStart(soildata) + + if (!terra::inMemory(soildata)) { + terra::readStart(soildata) + on.exit({ + try({ + terra::readStop(soildata) + }, silent = TRUE) + }, add = TRUE) + } # create template brick out <- terra::rast(soildata) @@ -191,15 +204,21 @@ run_rosetta.SpatRaster <- function(soildata, names(out) <- cnm out_info <- terra::writeStart(out, filename = file, overwrite = overwrite) + on.exit({ + try({ + out <- terra::writeStop(out) + }, silent = TRUE) + }, add = TRUE) + start_row <- seq(1, out_info$nrows, nrows) n_row <- diff(c(start_row, out_info$nrows + 1)) - if (cores > 1 && out_info$nrows*ncol(soildata) > 20000) { + if (cores > 1 && out_info$nrows*ncol(soildata) > core_thresh) { cls <- parallel::makeCluster(cores) on.exit(parallel::stopCluster(cls)) # TODO: can blocks be parallelized? - for(i in seq_along(start_row)) { + for (i in seq_along(start_row)) { if (n_row[i] > 0) { blockdata <- terra::readValues(soildata, row = start_row[i], nrows = n_row[i], dataframe = TRUE) ids <- 1:nrow(blockdata) @@ -207,7 +226,8 @@ run_rosetta.SpatRaster <- function(soildata, # run_rosetta is a "costly" function and not particularly fast, so in theory parallel would help # parallel within-block processing - X <- split(blockdata, rep(seq(from = 1, to = floor(length(ids) / 20000) + 1), each = 20000)[1:length(ids)]) + n <- floor(length(ids) / core_thresh / cores) + 1 + X <- split(blockdata, rep(seq(from = 1, to = n, each = n)))[1:length(ids)] r <- do.call('rbind', parallel::clusterApply(cls, X, function(x) rosettaPTF::run_rosetta(x, vars = vars, rosetta_version = rosetta_version))) @@ -216,7 +236,7 @@ run_rosetta.SpatRaster <- function(soildata, } } } else { - for(i in seq_along(start_row)) { + for (i in seq_along(start_row)) { if (n_row[i] > 0) { foo <- rosettaPTF::run_rosetta(terra::readValues(soildata, row = start_row[i], nrows = n_row[i], dataframe = TRUE), vars = vars, @@ -226,11 +246,9 @@ run_rosetta.SpatRaster <- function(soildata, } } - out <- terra::writeStop(out) - terra::readStop(soildata) + # replace NaN with NA_real_ (note: not compatible with calling writeStop() on.exit()) + # out <- terra::classify(out, cbind(NaN, NA_real_)) #terra::values(out)[is.nan(terra::values(out))] <- NA_real_ - # replace NaN with NA_real_ - terra::values(out)[is.nan(terra::values(out))] <- NA_real_ out } diff --git a/man/run_rosetta.Rd b/man/run_rosetta.Rd index a053df0..e5a7c98 100644 --- a/man/run_rosetta.Rd +++ b/man/run_rosetta.Rd @@ -21,8 +21,9 @@ vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -31,8 +32,9 @@ vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -41,8 +43,9 @@ vars = NULL, rosetta_version = 3, cores = 1, - file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + core_thresh = 20000L, + file = paste0(tempfile(), ".tif"), + nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) } @@ -57,9 +60,11 @@ \item{cores}{number of cores; used only for processing \emph{SpatRaster} or \emph{Raster*} input} +\item{core_thresh}{Magic number for determining processing chunk size. Default \code{20000L}. Used to calculate default \code{nrows}} + \item{file}{path to write incremental raster processing output for large inputs that do not fit in memory; passed to \code{terra::writeStart()} and used only for processing \emph{SpatRaster} or \emph{Raster*} input; defaults to a temporary file created by \code{tempfile()} if needed} -\item{nrows}{number of rows to use per block; passed to \code{terra::readValues()} \code{terra::writeValues()}; used only for processing \emph{SpatRaster} or \emph{Raster*} input; defaults to number of rows in dataset if needed} +\item{nrows}{number of rows to use per block chunk; passed to \code{terra::readValues()} and \code{terra::writeValues()}; used only for processing \emph{SpatRaster} or \emph{Raster*} inputs. Defaults to the total number of rows divided by the number of cells divided by \code{core_thresh}.} \item{overwrite}{logical; overwrite \code{file}? passed to \code{terra::writeStart()}; defaults to \code{TRUE} if needed} } diff --git a/tests/testthat/test-rosetta.R b/tests/testthat/test-rosetta.R index 12265ad..0cdcd3b 100644 --- a/tests/testthat/test-rosetta.R +++ b/tests/testthat/test-rosetta.R @@ -50,11 +50,11 @@ test_that("run on SSURGO data", { resrose <- rosettaPTF::run_rosetta(soildata[,varnames]) resrose$mukey <- resprop$mukey - rdf <- data.frame(mukey = as.numeric(terra::cats(res)[[1]][["category"]])) + rdf <- data.frame(mukey = as.numeric(terra::cats(res)[[1]][[1]])) rdf2 <- merge(rdf, resprop, by = "mukey", all.x = TRUE, sort = FALSE, incomparables = NA) rdf3 <- merge(rdf2, resrose, by = "mukey", all.x = TRUE, sort = FALSE, incomparables = NA) rdf3 <- rdf3[match(rdf3[["mukey"]], umukeys, incomparables = NA),][1:nrow(resprop),] - levels(res) <- data.frame(ID = 1:nrow(rdf3), rdf3) + terra::set.cats(res, 1, data.frame(ID = 1:nrow(rdf3), rdf3)) # @params x a SpatRaster with `levels()` set such that `cats(x)[[1]]` defines the mapping between raster values and one or more new attributes # @params columns character vector of column names to map from the categorical levels to raster values