From b6466fba502122d40ed274fa14a159ff6b726386 Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 09:30:03 -0800 Subject: [PATCH 1/6] run_rosetta: close connections on error - closes #11 --- R/run_rosetta.R | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/R/run_rosetta.R b/R/run_rosetta.R index 885204f..7ba291b 100644 --- a/R/run_rosetta.R +++ b/R/run_rosetta.R @@ -180,7 +180,15 @@ run_rosetta.SpatRaster <- function(soildata, file = paste0(tempfile(),".grd"), nrows = nrow(soildata), 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,6 +199,12 @@ 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)) @@ -226,11 +240,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 } From c1807711ff0999cd51c2da0b98515a7f66b7a48c Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 09:36:00 -0800 Subject: [PATCH 2/6] run_rosetta: better heuristics for `nrows` - for #12 - use `core_thresh` internally for within-block processing --- R/run_rosetta.R | 22 ++++++++++++++-------- man/run_rosetta.Rd | 13 +++++++++---- 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/R/run_rosetta.R b/R/run_rosetta.R index 7ba291b..82b9b49 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, + core_thresh = 20000L, file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + 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, + core_thresh = 20000L, file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + 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,8 +181,9 @@ run_rosetta.SpatRaster <- function(soildata, vars = NULL, rosetta_version = 3, cores = 1, + core_thresh = 20000L, file = paste0(tempfile(),".grd"), - nrows = nrow(soildata), + nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { if (!terra::inMemory(soildata)) { @@ -208,12 +213,12 @@ run_rosetta.SpatRaster <- function(soildata, 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) @@ -221,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))) @@ -230,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, diff --git a/man/run_rosetta.Rd b/man/run_rosetta.Rd index a053df0..8491c69 100644 --- a/man/run_rosetta.Rd +++ b/man/run_rosetta.Rd @@ -21,8 +21,9 @@ vars = NULL, rosetta_version = 3, cores = 1, + core_thresh = 20000L, file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -31,8 +32,9 @@ vars = NULL, rosetta_version = 3, cores = 1, + core_thresh = 20000L, file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -41,8 +43,9 @@ vars = NULL, rosetta_version = 3, cores = 1, + core_thresh = 20000L, file = paste0(tempfile(), ".grd"), - nrows = nrow(soildata), + 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} } From 130d53c525068653c84bf7b4b16bea732ae6476b Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 09:54:49 -0800 Subject: [PATCH 3/6] run_rosetta: fix test related to terra categories --- tests/testthat/test-rosetta.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 7aeb83b0816320d89e96ef6c9c82d035c333a6b5 Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 10:22:52 -0800 Subject: [PATCH 4/6] run_rosetta: use geotiff tempfiles --- R/run_rosetta.R | 6 +++--- man/run_rosetta.Rd | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/run_rosetta.R b/R/run_rosetta.R index 82b9b49..6a87a2a 100644 --- a/R/run_rosetta.R +++ b/R/run_rosetta.R @@ -124,7 +124,7 @@ run_rosetta.RasterStack <- function(soildata, rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(),".grd"), + 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 @@ -156,7 +156,7 @@ run_rosetta.RasterBrick <- function(soildata, rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(),".grd"), + file = paste0(tempfile(), ".tif"), nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { run_rosetta(terra::rast(soildata), @@ -182,7 +182,7 @@ run_rosetta.SpatRaster <- function(soildata, rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(),".grd"), + file = paste0(tempfile(), ".tif"), nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh), overwrite = TRUE) { diff --git a/man/run_rosetta.Rd b/man/run_rosetta.Rd index 8491c69..e5a7c98 100644 --- a/man/run_rosetta.Rd +++ b/man/run_rosetta.Rd @@ -22,7 +22,7 @@ rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(), ".grd"), + file = paste0(tempfile(), ".tif"), nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -33,7 +33,7 @@ rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(), ".grd"), + file = paste0(tempfile(), ".tif"), nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) @@ -44,7 +44,7 @@ rosetta_version = 3, cores = 1, core_thresh = 20000L, - file = paste0(tempfile(), ".grd"), + file = paste0(tempfile(), ".tif"), nrows = nrow(soildata)/(terra::ncell(soildata)/core_thresh), overwrite = TRUE ) From a5c2ba13d1efa6b363ec1a8069c8a8cd8cc2772a Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 10:24:05 -0800 Subject: [PATCH 5/6] version bump --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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. From bbdf5d7c7f6eca947e6bc1db63c885ed3316a908 Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Wed, 27 Dec 2023 11:15:05 -0800 Subject: [PATCH 6/6] silence --- R/find_python.R | 4 ++-- R/install_rosetta.R | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) 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()