Merge pull request #85 from carlesmila/master

NNDM modeldomain as SpatRaster (+ updates in docs)

R/nndm.R

@@ -4,7 +4,7 @@
 #' indices to perform a NNDM LOO CV for map validation.
 #' @author Carles Milà
 #' @param tpoints sf or sfc point object. Contains the training points samples.
-#' @param modeldomain sf polygon object defining the prediction area (see Details).
+#' @param modeldomain sf polygon object or SpatRaster defining the prediction area (see Details).
 #' @param ppoints sf or sfc point object. Contains the target prediction points.
 #' Optional. Alternative to modeldomain (see Details).
 #' @param samplesize numeric. How many points in the modeldomain should be sampled as prediction points?
@@ -34,12 +34,13 @@
 #' When \emph{phi} is set to "max", nearest neighbor distance matching is performed for the entire prediction area. Euclidean distances are used for projected
 #' and non-defined CRS, great circle distances are used for geographic CRS (units in meters).
 #'
-#' The \emph{modeldomain} is a sf polygon that defines the prediction area. The function takes a regular point sample (amount defined by \emph{samplesize)} from the spatial extent.
+#' The \emph{modeldomain} is either a sf polygon that defines the prediction area, or alternatively a SpatRaster out of which a polygon,
+#' transformed into the CRS of the training points, is defined as the outline of all non-NA cells.
+#' Then, the function takes a regular point sample (amount defined by \emph{samplesize)} from the spatial extent.
 #' As an alternative use \emph{ppoints} instead of \emph{modeldomain}, if you have already defined the prediction locations (e.g. raster pixel centroids).
 #' When using either \emph{modeldomain} or \emph{ppoints}, we advise to plot the study area polygon and the training/prediction points as a previous step to ensure they are aligned.
 #'
-#' @note NNDM is a variation of LOOCV and therefore may take a long time for large training data sets.
-#' A k-fold variant will be implemented shortly.
+#' @note NNDM is a variation of LOOCV and therefore may take a long time for large training data sets. See \code{\link{knndm}} for a more efficient k-fold variant of the method.
 #' @seealso \code{\link{geodist}}, \code{\link{knndm}}
 #' @references
 #' \itemize{
@@ -97,8 +98,8 @@
 #' plot(nndm_pred)
 #'
 #' ########################################################################
-#' # Example 3: Real- world example; using a modeldomain instead of previously
-#' # sampled prediction locations
+#' # Example 3: Real- world example; using a SpatRast modeldomain instead
+#' # of previously sampled prediction locations
 #' ########################################################################
 #' \dontrun{
 #' library(sf)
@@ -112,15 +113,11 @@
 #' pts <- st_as_sf(pts,coords=c("Easting","Northing"))
 #' st_crs(pts) <- 26911
 #' studyArea <- rast(system.file("extdata","predictors_2012-03-25.tif",package="CAST"))
-#' studyArea[!is.na(studyArea)] <- 1
-#' studyArea <- as.polygons(studyArea, values = FALSE, na.all = TRUE) |>
-#'     st_as_sf() |>
-#'     st_union()
 #' pts <- st_transform(pts, crs = st_crs(studyArea))
-#' plot(studyArea)
-#' plot(st_geometry(pts), add = TRUE, col = "red")
+#' terra::plot(studyArea[["DEM"]])
+#' terra::plot(vect(pts), add = T)
 #'
-#' nndm_folds <- nndm(pts, modeldomain= studyArea)
+#' nndm_folds <- nndm(pts, modeldomain = studyArea)
 #' plot(nndm_folds)
 #'
 #' #use for cross-validation:
@@ -144,10 +141,22 @@ nndm <- function(tpoints, modeldomain = NULL, ppoints = NULL, samplesize = 1000,
   # create sample points from modeldomain
   if(is.null(ppoints)&!is.null(modeldomain)){
 
-    # Check modeldomain is indeed a polygon
-    if(!any(c("sfc", "sf") %in% class(modeldomain))){
-      stop("modeldomain must be a sf/sfc object.")
-    }else if(!any(class(sf::st_geometry(modeldomain)) %in% c("sfc_POLYGON", "sfc_MULTIPOLYGON"))){
+    # Check modeldomain is indeed a sf/SpatRaster
+    if(!any(c("sfc", "sf", "SpatRaster") %in% class(modeldomain))){
+      stop("modeldomain must be a sf/sfc object or a 'SpatRaster' object.")
+    }
+
+    # If modeldomain is a SpatRaster, transform into polygon
+    if(any(class(modeldomain) == "SpatRaster")){
+      modeldomain[!is.na(modeldomain)] <- 1
+      modeldomain <- terra::as.polygons(modeldomain, values = FALSE, na.all = TRUE) |>
+        sf::st_as_sf() |>
+        sf::st_union()
+      modeldomain <- sf::st_transform(modeldomain, crs = sf::st_crs(tpoints))
+    }
+
+    # Check modeldomain is indeed a polygon sf
+    if(!any(class(sf::st_geometry(modeldomain)) %in% c("sfc_POLYGON", "sfc_MULTIPOLYGON"))){
       stop("modeldomain must be a sf/sfc polygon object.")
     }
 

tests/testthat/test-nndm.R

@@ -37,7 +37,7 @@ test_that("NNDM detects wrong data and geometry types", {
 
   # model domain
   expect_error(suppressWarnings(nndm(tpoints_sfc, modeldomain = 1)),
-               "modeldomain must be a sf/sfc object.")
+               "modeldomain must be a sf/sfc object or a 'SpatRaster' object.")
   expect_error(nndm(tpoints_sfc, modeldomain = ppoints_sfc),
                "modeldomain must be a sf/sfc polygon object.")
 })
@@ -106,5 +106,24 @@ test_that("NNDM yields the expected results for all CRS", {
   # Geographic
   tpoints_sf_4326 <- sf::st_set_crs(tpoints_sfc, 4326)
   ppoints_sf_4326 <- sf::st_set_crs(ppoints_sfc, 4326)
-  expect_equal(as.numeric(nndm(tpoints_sf_4326, ppoints = ppoints_sf_4326, phi = 1000000)$Gjstar[20], 4), 355614.94)
+  expect_equal(as.numeric(nndm(tpoints_sf_4326, ppoints = ppoints_sf_4326, phi = 1000000)$Gjstar[20]), 355614.94)
+})
+
+test_that("NNDM yields the expected results with SpatRast modeldomain", {
+
+  set.seed(1234)
+
+  # prepare sample data
+  dat <- readRDS(system.file("extdata","Cookfarm.RDS",package="CAST"))
+  dat <- terra::aggregate(dat[,c("DEM","TWI", "NDRE.M", "Easting", "Northing","VW")],
+                          by=list(as.character(dat$SOURCEID)),mean)
+  pts <- dat[,-1]
+  pts <- sf::st_as_sf(pts,coords=c("Easting","Northing"))
+  sf::st_crs(pts) <- 26911
+  studyArea <- terra::rast(system.file("extdata","predictors_2012-03-25.tif",package="CAST"))
+  pts <- sf::st_transform(pts, crs = sf::st_crs(studyArea))
+
+  nndm_folds <- nndm(pts, modeldomain = studyArea, phi = 150)
+  expect_equal(as.numeric(nndm(pts, modeldomain = studyArea, phi = 150)$Gjstar[5]), 63.828663)
+
 })