From 5d48e4d4c83853aa482eef9d7964981fd03bd638 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Darius=20G=C3=B6rgen?= Date: Tue, 20 Apr 2021 15:25:27 +0200 Subject: [PATCH] coerce terra SpatRaster to RasterStack --- DESCRIPTION | 4 +- R/aoa.R | 12 +++++- R/calibrate_aoa.R | 15 ++++++++ man/aoa.Rd | 2 +- vignettes/terra-showcase.Rmd | 71 ++++++++++++++++++++++++++++++++++++ 5 files changed, 100 insertions(+), 4 deletions(-) create mode 100644 vignettes/terra-showcase.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index a9bd236d..63a345d7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,6 +18,6 @@ Depends: R (>= 3.5.0) Imports: caret, stats, utils, ggplot2, graphics, reshape, FNN, plyr, zoo, methods, grDevices, data.table, lattice Suggests: doParallel, randomForest, lubridate, raster, sp, knitr, mapview, rmarkdown, sf, scales, parallel, latticeExtra, virtualspecies, gridExtra, - viridis, rgeos, stars, scam -RoxygenNote: 7.1.0 + viridis, rgeos, stars, scam, terra +RoxygenNote: 7.1.1 VignetteBuilder: knitr diff --git a/R/aoa.R b/R/aoa.R index 9845cdbd..c6145206 100644 --- a/R/aoa.R +++ b/R/aoa.R @@ -8,7 +8,7 @@ #' variable importance of the machine learning algorithm used for model training. #' The AOA is derived by applying a threshold on the DI which is the (outlier-removed) #' maximum DI of the cross-validated training data. -#' @param newdata A RasterStack, RasterBrick or data.frame containing the data +#' @param newdata A RasterStack, RasterBrick, SpatRaster or data.frame containing the data #' the model was meant to make predictions for. #' @param model A train object created with caret used to extract weights from (based on variable importance) as well as cross-validation folds #' @param cl A cluster object e.g. created with doParallel. Should only be used if newdata is large. @@ -119,12 +119,20 @@ aoa <- function(newdata, } } as_stars <- FALSE + as_terra <- FALSE if (inherits(newdata, "stars")) { if (!requireNamespace("stars", quietly = TRUE)) stop("package stars required: install that first") newdata = methods::as(newdata, "Raster") as_stars <- TRUE } + + if (inherits(newdata, "SpatRaster")) { + if (!requireNamespace("terra", quietly = TRUE)) + stop("package terra required: install that first") + newdata = methods::as(newdata, "Raster") + as_terra <- TRUE + } #### Prepare output as either as RasterLayer or vector: out <- NA if (inherits(newdata, "Raster")) @@ -299,6 +307,8 @@ aoa <- function(newdata, out <- raster::stack(out,AOA) if (as_stars) out <- split(stars::st_as_stars(out), "band") + if(as_terra) + out = methods::as(out, "SpatRaster") }else{ out <- DI_out AOA <- rep(1,length(out)) diff --git a/R/calibrate_aoa.R b/R/calibrate_aoa.R index 9ea6ded5..3287bf61 100644 --- a/R/calibrate_aoa.R +++ b/R/calibrate_aoa.R @@ -63,6 +63,7 @@ calibrate_aoa <- function(AOA,model, window.size=5, calib="scam",multiCV=FALSE, length.out = 10, maskAOA=TRUE, showPlot=TRUE,k=6,m=2){ as_stars <- FALSE + as_terra <- FALSE if (inherits(AOA, "stars")) { if (!requireNamespace("stars", quietly = TRUE)) stop("package stars required: install that first") @@ -72,6 +73,15 @@ calibrate_aoa <- function(AOA,model, window.size=5, calib="scam",multiCV=FALSE, as_stars <- TRUE } + if (inherits(AOA, "SpatRaster")) { + if (!requireNamespace("terra", quietly = TRUE)) + stop("package terra required: install that first") + attr <- attributes(AOA)[c("aoa_stats","TrainDI")] + AOA <- methods::as(AOA, "Raster") + attributes(AOA)<- c(attributes(AOA),attr) + as_terra <- TRUE + } + if(multiCV){ preds_all <- data.frame() train_predictors <- model$trainingData[,-which(names(model$trainingData)==".outcome")] @@ -274,6 +284,11 @@ calibrate_aoa <- function(AOA,model, window.size=5, calib="scam",multiCV=FALSE, AOA <- split(stars::st_as_stars(AOA), "band") attributes(AOA)<- c(attributes(AOA),attr) } + + if(as_terra){ + AOA <- methods::as(AOA, "SpatRaster") + attributes(AOA)<- c(attributes(AOA),attr) + } names(AOA)[names(AOA)=="expectedError"] <- paste0("expected_",model$metric) #return(AOA) diff --git a/man/aoa.Rd b/man/aoa.Rd index 0a39edb7..52f71fc1 100644 --- a/man/aoa.Rd +++ b/man/aoa.Rd @@ -16,7 +16,7 @@ aoa( ) } \arguments{ -\item{newdata}{A RasterStack, RasterBrick or data.frame containing the data +\item{newdata}{A RasterStack, RasterBrick, SpatRaster or data.frame containing the data the model was meant to make predictions for.} \item{model}{A train object created with caret used to extract weights from (based on variable importance) as well as cross-validation folds} diff --git a/vignettes/terra-showcase.Rmd b/vignettes/terra-showcase.Rmd new file mode 100644 index 00000000..f26a901c --- /dev/null +++ b/vignettes/terra-showcase.Rmd @@ -0,0 +1,71 @@ +--- +title: "Showcase with terra SpatRaster" +date: "`r Sys.Date()`" +output: + rmarkdown::html_document: + toc: true + theme: united +vignette: > + %\VignetteIndexEntry{Showcase with terra SpatRaster} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, echo=FALSE} +knitr::opts_chunk$set(fig.width = 8.83) +``` + +--- + +# Introduction + +This vignette is meant to showcase that the CAST package now is able to handle +terra SpatRaster objects. Changes to the code were implemented only in the `aoa()` and +`calibrate_aoa()` function and are restricted to coercing a SpatRaster input +to a RasterStack as well as coercing the output back to SpatRaster. +Below is a minimal example based on the AOA-tutorial vignette. + +```{r message = FALSE, warning=FALSE} +library(CAST) +library(terra) +library(raster) +library(lubridate) +library(methods) +``` + +```{r} +predictors_sp <- stack(system.file("extdata","predictors_2012-03-25.grd",package="CAST")) +data <- get(load(system.file("extdata","Cookfarm.RData",package="CAST"))) +trainDat <- data[data$altitude==-0.3& + year(data$Date)==2012& + week(data$Date)%in%c(10:12),] + +library(caret) +predictors <- c("DEM","TWI","Precip_cum","cday", + "MaxT_wrcc","Precip_wrcc","BLD", + "Northing","Easting","NDRE.M") +set.seed(10) +model <- train(trainDat[,predictors],trainDat$VW, + method="rf",tuneGrid=data.frame("mtry"=2), + importance=TRUE,ntree=50, + trControl=trainControl(method="cv",number=3)) +``` + +From here we will coerce the Raster object to SpatRaster objects to show +that the aoa and calibrate_aoa functions work on SpatRasters as well. The +application might not make sense on a thematic level. The example was drafted +just to show that the functions accept and return SpatRaster objects. + +```{r} +predictors_terra = as(predictors_sp, "SpatRaster") +AOA <- aoa(predictors_terra, model) +class(AOA) +attributes(AOA)$aoa_stats +``` + +```{r} +AOA_calib <- calibrate_aoa(AOA,model,window.size = 5,length.out = 5, multiCV=TRUE,showPlot=FALSE) +class(AOA_calib$AOA) +AOA_calib$plot +``` +