diff --git a/.Rbuildignore b/.Rbuildignore
index c1adbe6..6ef99de 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -7,3 +7,4 @@
^CRAN-RELEASE$
^CRAN-SUBMISSION$
^man/figures
+^\.github$
diff --git a/.github/.gitignore b/.github/.gitignore
new file mode 100644
index 0000000..2d19fc7
--- /dev/null
+++ b/.github/.gitignore
@@ -0,0 +1 @@
+*.html
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
new file mode 100644
index 0000000..a3ac618
--- /dev/null
+++ b/.github/workflows/R-CMD-check.yaml
@@ -0,0 +1,49 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+ push:
+ branches: [main, master]
+ pull_request:
+ branches: [main, master]
+
+name: R-CMD-check
+
+jobs:
+ R-CMD-check:
+ runs-on: ${{ matrix.config.os }}
+
+ name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+ strategy:
+ fail-fast: false
+ matrix:
+ config:
+ - {os: macos-latest, r: 'release'}
+ - {os: windows-latest, r: 'release'}
+ - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'}
+ - {os: ubuntu-latest, r: 'release'}
+ - {os: ubuntu-latest, r: 'oldrel-1'}
+
+ env:
+ GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+ R_KEEP_PKG_SOURCE: yes
+
+ steps:
+ - uses: actions/checkout@v3
+
+ - uses: r-lib/actions/setup-pandoc@v2
+
+ - uses: r-lib/actions/setup-r@v2
+ with:
+ r-version: ${{ matrix.config.r }}
+ http-user-agent: ${{ matrix.config.http-user-agent }}
+ use-public-rspm: true
+
+ - uses: r-lib/actions/setup-r-dependencies@v2
+ with:
+ extra-packages: any::rcmdcheck
+ needs: check
+
+ - uses: r-lib/actions/check-r-package@v2
+ with:
+ upload-snapshots: true
diff --git a/DESCRIPTION b/DESCRIPTION
index acb5382..2d17419 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: envi
Type: Package
Title: Environmental Interpolation using Spatial Kernel Density Estimation
-Version: 0.1.15
-Date: 2022-08-30
+Version: 0.1.16.9000
+Date: 2022-12-16
Authors@R:
c(person(given = "Ian D.",
family = "Buller",
@@ -28,7 +28,7 @@ Description: Estimates an ecological niche using occurrence data, covariates, an
License: Apache License (>= 2.0)
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.2
Depends:
R (>= 3.5.0)
Imports:
@@ -43,15 +43,13 @@ Imports:
grDevices,
iterators,
pls,
- raster,
- rgeos,
ROCR,
- sp,
+ sf,
sparr,
spatstat.geom,
- stats
+ stats,
+ terra
Suggests:
- maptools,
R.rsp,
RStoolbox,
spatstat.data,
diff --git a/NAMESPACE b/NAMESPACE
index a1e529b..2303b28 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -34,23 +34,10 @@ importFrom(graphics,plot.new)
importFrom(graphics,title)
importFrom(iterators,icount)
importFrom(pls,cvsegments)
-importFrom(raster,crs)
-importFrom(raster,cut)
-importFrom(raster,extract)
-importFrom(raster,image)
-importFrom(raster,projectRaster)
-importFrom(raster,raster)
-importFrom(raster,rasterToPoints)
-importFrom(raster,reclassify)
-importFrom(raster,values)
-importFrom(rgeos,gBuffer)
-importFrom(sp,CRS)
-importFrom(sp,Polygon)
-importFrom(sp,Polygons)
-importFrom(sp,SpatialPolygons)
-importFrom(sp,bbox)
-importFrom(sp,coordinates)
-importFrom(sp,gridded)
+importFrom(sf,st_bbox)
+importFrom(sf,st_buffer)
+importFrom(sf,st_coordinates)
+importFrom(sf,st_polygon)
importFrom(sparr,risk)
importFrom(spatstat.geom,as.solist)
importFrom(spatstat.geom,im.apply)
@@ -65,3 +52,12 @@ importFrom(spatstat.geom,superimpose)
importFrom(stats,median)
importFrom(stats,na.omit)
importFrom(stats,sd)
+importFrom(terra,classify)
+importFrom(terra,crds)
+importFrom(terra,crs)
+importFrom(terra,extract)
+importFrom(terra,image)
+importFrom(terra,project)
+importFrom(terra,rast)
+importFrom(terra,res)
+importFrom(terra,values)
diff --git a/NEWS.md b/NEWS.md
index b202def..c10e582 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,14 @@
# envi (development version)
+# envi v0.1.16.9000
+* Migrated R-spatial dependencies
+* Replaced `raster`, `rgeos`, and `sp` packages in Imports with `terra` and `sf` because of imminent package retirements
+* Removed `maptools` from Suggests (replaced with new internal function `as.im.SpatRaster()`)
+* Thank you, [Roger Bivand](https://github.com/rsbivand), for the notice. Relates to [ndi Issue #3](https://github.com/lance-waller-lab/envi/issues/3)
+* Note: `raster` is a dependency of `RStoolbox` (at present) which is used in the vignette
+* Updated test, examples, vignette, and documentation throughout
+* Added GitHub R-CMD-check
+
# envi v0.1.15
* Uwe Ligges suggested (2022-08-26) that some parallel cluster is not cleanly closed in the tests
* Added `future::plan(future::multisession)` in tests to remove the files in temp directory
diff --git a/R/div_plot.R b/R/div_plot.R
index 80d0e99..971be00 100644
--- a/R/div_plot.R
+++ b/R/div_plot.R
@@ -1,8 +1,8 @@
-#' Prepare an 'im' or 'raster' object for plotting with diverging color palette
+#' Prepare an 'im' or 'SpatRaster' object for plotting with diverging color palette
#'
-#' Internal function to convert 'im' object or 'RasterLayer' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_obs}}, \code{\link{plot_predict}}, and \code{\link{plot_perturb}} functions.
+#' Internal function to convert 'im' object or 'SpatRaster' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_obs}}, \code{\link{plot_predict}}, and \code{\link{plot_perturb}} functions.
#'
-#' @param input An object of class 'im' or 'RasterLayer' from the \code{\link{lrren}} function.
+#' @param input An object of class 'im' or 'SpatRaster' from the \code{\link{lrren}} function.
#' @param plot_cols Character string of length three (3) specifying the colors for plotting: 1) presence, 2) neither, and 3) absence from the \code{\link{plot_obs}} function.
#' @param midpoint Numeric. The value to center the diverging color palette.
#' @param thresh_up Numeric. The upper value to concatenate the color key. The default (NULL) uses the maximum value from \code{input}.
@@ -20,7 +20,7 @@
#' }
#'
#' @importFrom grDevices colorRampPalette
-#' @importFrom raster raster
+#' @importFrom terra rast
#'
#' @keywords internal
@@ -33,15 +33,15 @@ div_plot <- function(input,
# Inputs
if (inherits(input, "im")) {
- out <- raster::raster(input)
+ out <- terra::rast(input)
} else { out <- input }
if (length(cols) != 3) {
stop("The 'cols' argument must be a vector of length 3")
}
- min_raw_value <- min(out[is.finite(out)], na.rm = TRUE) # minimum absolute value of raster
- max_raw_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of raster
+ min_raw_value <- min(out[is.finite(out)], na.rm = TRUE) # minimum absolute value of SpatRaster
+ max_raw_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of SpatRaster
# Restrict spurious log relative risk values
if (!is.null(thresh_low)) {
@@ -60,8 +60,8 @@ div_plot <- function(input,
# Identify ramp above and below midpoint
lowerhalf <- length(out[out < midpoint & !is.na(out)]) # values below 0
upperhalf <- length(out[out > midpoint & !is.na(out)]) # values above 0
- min_absolute_value <- min(out[is.finite(out)], na.rm = TRUE) # minimum absolute value of raster
- max_absolute_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of raster
+ min_absolute_value <- min(out[is.finite(out)], na.rm = TRUE) # minimum absolute value of SpatRaster
+ max_absolute_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of SpatRaster
# Color ramp parameters
## Colors
diff --git a/R/lrren.R b/R/lrren.R
index 14c9c36..a84cfcf 100644
--- a/R/lrren.R
+++ b/R/lrren.R
@@ -70,12 +70,11 @@
#' @importFrom grDevices chull
#' @importFrom iterators icount
#' @importFrom pls cvsegments
-#' @importFrom raster extract raster
-#' @importFrom rgeos gBuffer
-#' @importFrom sp bbox coordinates Polygon Polygons SpatialPolygons
+#' @importFrom sf st_bbox st_buffer st_coordinates st_polygon
#' @importFrom sparr risk
#' @importFrom spatstat.geom owin ppp
#' @importFrom stats na.omit
+#' @importFrom terra extract rast values
#' @export
#'
#' @examples
@@ -89,8 +88,8 @@
#' grad <- spatstat.data::bei.extra[[2]]
#' elev$v <- scale(elev)
#' grad$v <- scale(grad)
-#' elev_raster <- raster::raster(elev)
-#' grad_raster <- raster::raster(grad)
+#' elev_raster <- terra::rast(elev)
+#' grad_raster <- terra::rast(grad)
#'
#' # Presence data
#' presence <- spatstat.data::bei
@@ -115,8 +114,10 @@
#' obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
#'
#' # Prediction Data
-#' predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-#' predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+#' predict_xy <- terra::crds(elev_raster)
+#' predict_locs <- as.data.frame(predict_xy)
+#' predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+#' predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
#'
#' # Run lrren
#' test_lrren <- lrren(obs_locs = obs_locs,
@@ -138,7 +139,7 @@ lrren <- function(obs_locs,
n_core = 2,
poly_buffer = NULL,
obs_window = NULL,
- verbose = FALSE,
+ verbose = FALSE,
...) {
if (verbose == TRUE) { message("Estimating relative risk surfaces\n") }
@@ -148,17 +149,15 @@ lrren <- function(obs_locs,
# Compute spatial windows
## Calculate inner boundary polygon (extent of presence and absence locations in environmental space)
inner_chull <- concaveman::concaveman(as.matrix(obs_locs[ , 5:6]))
- inner_chull_pts <- sp::coordinates(inner_chull)
- inner_chull_pts <- rbind(inner_chull_pts, inner_chull_pts[1, ])
- inner_chull_poly <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(inner_chull_pts)), 1)))
-
+ inner_chull_poly <- sf::st_polygon(list(inner_chull))
+
if (is.null(poly_buffer)) {
- poly_buffer <- abs(min(diff(sp::bbox(inner_chull_poly)[1, ]), diff(sp::bbox(inner_chull_poly)[2, ])) / 100)
+ poly_buffer <- abs(min(diff(sf::st_bbox(inner_chull_poly)[c(1,3)]), diff(sf::st_bbox(inner_chull_poly)[c(2,4)])) / 100)
}
# add small buffer around polygon to include boundary points
- inner_chull_poly_buffer <- rgeos::gBuffer(inner_chull_poly, width = poly_buffer, byid = TRUE)
- inner_poly <- inner_chull_poly_buffer@polygons[[1]]@Polygons[[1]]@coords
+ inner_chull_poly_buffer <- sf::st_buffer(inner_chull_poly, dist = poly_buffer, byid = TRUE)
+ inner_poly <- sf::st_polygon(list(as.matrix(inner_chull_poly_buffer)))
if (is.null(predict_locs)) {
outer_chull_poly <- inner_chull_poly_buffer
@@ -166,17 +165,16 @@ lrren <- function(obs_locs,
} else {
## Calculate outer boundary polygon (full extent of geographical extent in environmental space)
if (nrow(predict_locs) > 5000000) { # convex hull
- outer_chull <- grDevices::chull(x = stats::na.omit(predict_locs)[ , 3], y = stats::na.omit(predict_locs)[ , 4])
- outer_chull_pts <- predict_locs[c(outer_chull, outer_chull[1]), 3:4]
+ predict_locs_woNAs <- stats::na.omit(predict_locs)
+ outer_chull <- grDevices::chull(x = predict_locs_woNAs[ , 3], y = predict_locs_woNAs[ , 4])
+ outer_chull_pts <- predict_locs_woNAs[c(outer_chull, outer_chull[1]), 3:4]
} else { # concave hull
- outer_chull <- concaveman::concaveman(as.matrix(stats::na.omit(predict_locs)[ , 3:4]))
- outer_chull_pts <- sp::coordinates(outer_chull)
+ outer_chull_pts <- concaveman::concaveman(as.matrix(stats::na.omit(predict_locs[ , 3:4])))
}
- outer_chull_pts <- rbind(outer_chull_pts, outer_chull_pts[1, ])
- outer_chull_poly <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(outer_chull_pts)), 1)))
- #add small buffer around polygon to include boundary points
- outer_chull_poly_buffer <- rgeos::gBuffer(outer_chull_poly, width = poly_buffer, byid = TRUE)
- outer_poly <- outer_chull_poly_buffer@polygons[[1]]@Polygons[[1]]@coords #extract coordinates of new polygon
+ outer_chull_poly <- sf::st_polygon(list(as.matrix(outer_chull_pts)))
+ # add small buffer around polygon to include boundary points
+ outer_chull_poly_buffer <- sf::st_buffer(outer_chull_poly, dist = poly_buffer, byid = TRUE)
+ outer_poly <- sf::st_polygon(list(as.matrix(outer_chull_poly_buffer)))
}
if (conserve == FALSE & is.null(predict_locs)) {
@@ -185,8 +183,8 @@ lrren <- function(obs_locs,
if (conserve == TRUE) { window_poly <- inner_poly } else { window_poly <- outer_poly }
if (is.null(obs_window)) {
- wind <- spatstat.geom::owin(poly = list(x = rev(window_poly[ , 1]),
- y = rev(window_poly[ , 2])))
+ wind <- spatstat.geom::owin(poly = list(x = rev(sf::st_coordinates(window_poly)[ , 1]),
+ y = rev(sf::st_coordinates(window_poly)[ , 2])))
} else { wind <- obs_window }
# Input Preparation
@@ -207,7 +205,7 @@ lrren <- function(obs_locs,
obs <- sparr::risk(f = ppp_presence,
g = ppp_absence,
tolerate = TRUE,
- verbose = verbose,
+ verbose = verbose,
...)
bandw <- obs$f$h0
@@ -228,17 +226,17 @@ lrren <- function(obs_locs,
# Project relative risk surface into geographic space
if (verbose == TRUE) { message("Predicting area of interest") }
- # Convert to semi-continuous raster
- rr_raster <- raster::raster(obs$rr)
+ # Convert to semi-continuous SpatRaster
+ rr_raster <- terra::rast(obs$rr)
- # Convert to categorical raster
- pval_raster <- raster::raster(obs$P)
+ # Convert to categorical SpatRaster
+ pval_raster <- terra::rast(obs$P)
# Prediction locations
extract_points <- cbind(predict_locs[ , 3], predict_locs[ , 4])
extract_predict <- data.frame("predict_locs" = predict_locs,
- "rr" = raster::extract(rr_raster, extract_points),
- "pval" = raster::extract(pval_raster, extract_points))
+ "rr" = terra::extract(rr_raster, extract_points)[ , 1],
+ "pval" = terra::extract(pval_raster, extract_points)[ , 1])
output <- list("obs" = obs,
"presence" = ppp_presence,
@@ -326,19 +324,19 @@ lrren <- function(obs_locs,
rand_lrr <- sparr::risk(f = ppp_presence_training,
g = ppp_absence_training,
tolerate = TRUE,
- verbose = FALSE,
+ verbose = FALSE,
...)
- ##### Convert to semi-continuous raster
- rr_raster <- raster::raster(rand_lrr$rr)
- rr_raster[is.na(rr_raster[])] <- 0 # if NA, assigned null value (log(rr) = 0)
+ ##### Convert to semi-continuous SpatRaster
+ rr_raster <- terra::rast(rand_lrr$rr)
+ rr_raster[is.na(terra::values(rr_raster))] <- 0 # if NA, assigned null value (log(rr) = 0)
##### Predict testing dataset
extract_testing <- testing[ , 5:6]
##### Output for each k-fold
###### Record category (semi-continuous) of testing data
- cv_predictions_rr <- raster::extract(rr_raster, extract_testing)
+ cv_predictions_rr <- terra::extract(rr_raster, extract_testing)[ , 2]
cv_labels <- testing[ , 4] # Record labels (marks) of testing data
par_results <- list("cv_predictions_rr" = cv_predictions_rr,
diff --git a/R/package.R b/R/package.R
index 328eff2..20f8b63 100644
--- a/R/package.R
+++ b/R/package.R
@@ -30,7 +30,7 @@
#' @aliases envi-package envi
#' @docType package
#'
-#' @section Dependencies: The 'envi' package relies heavily upon \code{\link{sparr}}, \code{\link{spatstat.geom}}, and \code{\link{raster}}. For a single species (presence/absence data), the spatial relative risk function uses the \code{\link[sparr]{risk}} function. Cross-validation is can be performed in parallel using the \code{\link{future}}, \code{\link{doFuture}}, \code{\link{doRNG}}, and \code{\link{foreach}} packages. Spatial perturbation is performed using the \code{\link[spatstat.geom]{rjitter}} function. Basic visualizations rely on the \code{\link[spatstat.geom]{plot.ppp}} and \code{\link[fields]{image.plot}} functions.
+#' @section Dependencies: The 'envi' package relies heavily upon \code{\link{sparr}}, \code{\link{spatstat.geom}}, \code{\link{sf}}, and \code{\link{terra}}. For a single species (presence/absence data), the spatial relative risk function uses the \code{\link[sparr]{risk}} function. Cross-validation is can be performed in parallel using the \code{\link{future}}, \code{\link{doFuture}}, \code{\link{doRNG}}, and \code{\link{foreach}} packages. Spatial perturbation is performed using the \code{\link[spatstat.geom]{rjitter}} function. Basic visualizations rely on the \code{\link[spatstat.geom]{plot.ppp}} and \code{\link[fields]{image.plot}} functions.
#'
#' @author Ian D. Buller\cr \emph{Environmental Health Sciences, Emory University, Atlanta, Georgia, USA.}\cr
#'
@@ -50,10 +50,9 @@ NULL
#' @importFrom grDevices chull colorRampPalette
#' @importFrom iterators icount
#' @importFrom pls cvsegments
-#' @importFrom raster crs cut extract image projectRaster raster rasterToPoints reclassify values
-#' @importFrom rgeos gBuffer
#' @importFrom ROCR performance prediction
-#' @importFrom sp bbox coordinates CRS gridded Polygon Polygons SpatialPolygons
+#' @importFrom sf st_bbox st_buffer st_coordinates st_polygon
#' @importFrom spatstat.geom as.solist im.apply marks owin pixellate plot.ppp ppp rjitter setmarks superimpose
#' @importFrom stats median na.omit sd
+#' @importFrom terra crds crs image project rast res classify values
NULL
diff --git a/R/perlrren.R b/R/perlrren.R
index 4ccfc67..8de29ab 100644
--- a/R/perlrren.R
+++ b/R/perlrren.R
@@ -58,9 +58,9 @@
#' @importFrom foreach %do% %dopar% foreach setDoPar
#' @importFrom future multisession plan
#' @importFrom iterators icount
-#' @importFrom raster crs
#' @importFrom spatstat.geom as.solist im.apply marks owin ppp rjitter superimpose
#' @importFrom stats median sd
+#' @importFrom terra crs
#' @export
#'
#' @examples
@@ -115,10 +115,6 @@ perlrren <- function(obs_ppp,
verbose = FALSE,
...) {
- if (!identical(raster::crs(obs_ppp), raster::crs(covariates))) {
- stop("The arguments 'obs_ppp' and 'covariates' must have the same coordinate reference system")
- }
-
if (is.null(radii)) {
radii <- rep(0, nlevels(spatstat.geom::marks(obs_ppp)[ , 5]))
message("The argument 'radii' is unspecified and the observation coordinates are not perturbed")
@@ -143,7 +139,7 @@ perlrren <- function(obs_ppp,
### Set function used in foreach
if (parallel == TRUE) {
oldplan <- doFuture::registerDoFuture()
- on.exit(with(oldplan, foreach::setDoPar(fun=fun, data=data, info=info)), add = TRUE)
+ on.exit(with(oldplan, foreach::setDoPar(fun = fun, data = data, info = info)), add = TRUE)
future::plan(future::multisession, workers = n_core)
`%fun%` <- doRNG::`%dorng%`
} else { `%fun%` <- foreach::`%do%` }
@@ -246,8 +242,8 @@ perlrren <- function(obs_ppp,
# Project relative risk surface into geographic space
if (verbose == TRUE) { message("\nPredicting area of interest") }
window_poly <- out_par[[3]][[1]]
- wind <- spatstat.geom::owin(poly = list(x = rev(window_poly[ , 1]),
- y = rev(window_poly[ , 2])))
+ wind <- spatstat.geom::owin(poly = list(x = rev(sf::st_coordinates(window_poly)[ , 1]),
+ y = rev(sf::st_coordinates(window_poly)[ , 2])))
xxxxx <- spatstat.geom::ppp(x = predict_locs[ , 3],
y = predict_locs[ , 4],
diff --git a/R/plot_cv.R b/R/plot_cv.R
index e841bef..eae6e01 100644
--- a/R/plot_cv.R
+++ b/R/plot_cv.R
@@ -24,8 +24,8 @@
#' grad <- spatstat.data::bei.extra[[2]]
#' elev$v <- scale(elev)
#' grad$v <- scale(grad)
-#' elev_raster <- raster::raster(elev)
-#' grad_raster <- raster::raster(grad)
+#' elev_raster <- terra::rast(elev)
+#' grad_raster <- terra::rast(grad)
#'
#' # Presence data
#' presence <- spatstat.data::bei
@@ -49,14 +49,8 @@
#' obs_locs$id <- seq(1, nrow(obs_locs), 1)
#' obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
#'
-#' # Prediction Data
-#' predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-#' predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
-#'
#' # Run lrren
#' test_lrren <- lrren(obs_locs = obs_locs,
-#' predict_locs = predict_locs,
-#' predict = TRUE,
#' cv = TRUE)
#'
#' # Run plot_cv
diff --git a/R/plot_obs.R b/R/plot_obs.R
index aefb5ca..e3b1742 100644
--- a/R/plot_obs.R
+++ b/R/plot_obs.R
@@ -14,8 +14,8 @@
#'
#' @importFrom fields image.plot
#' @importFrom graphics par
-#' @importFrom raster cut raster values
#' @importFrom spatstat.geom plot.ppp setmarks superimpose
+#' @importFrom terra rast values
#' @export
#'
#' @examples
@@ -29,8 +29,8 @@
#' grad <- spatstat.data::bei.extra[[2]]
#' elev$v <- scale(elev)
#' grad$v <- scale(grad)
-#' elev_raster <- raster::raster(elev)
-#' grad_raster <- raster::raster(grad)
+#' elev_raster <- terra::rast(elev)
+#' grad_raster <- terra::rast(grad)
#'
#' # Presence data
#' presence <- spatstat.data::bei
@@ -54,15 +54,9 @@
#' obs_locs$id <- seq(1, nrow(obs_locs), 1)
#' obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
#'
-#' # Prediction Data
-#' predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-#' predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
-#'
#' # Run lrren
#' test_lrren <- lrren(obs_locs = obs_locs,
-#' predict_locs = predict_locs,
-#' predict = TRUE,
-#' cv = TRUE)
+#' cv = FALSE)
#'
#' # Run plot_obs
#' plot_obs(input = test_lrren)
@@ -136,12 +130,12 @@ plot_obs <- function(input,
cex.axis = 0.67))
# Plot 3: Significant p-values
- pvalp <- raster::raster(input$out$obs$P) # create raster
- pvalp <- raster::cut(pvalp,
- breaks = c(-Inf, alpha / 2, 1 - alpha / 2, Inf),
- right = FALSE)
+ pvalp <- terra::rast(input$out$obs$P) # create raster
+ terra::values(pvalp) <- cut(terra::values(pvalp),
+ breaks = c(-Inf, alpha / 2, 1 - alpha / 2, Inf),
+ right = FALSE)
- if (all(raster::values(pvalp)[!is.na(raster::values(pvalp))] == 2)) {
+ if (all(terra::values(pvalp)[!is.na(terra::values(pvalp))] == 2)) {
pcols <- plot_cols[2]
brp <- c(1, 3)
atp <- 2
diff --git a/R/plot_perturb.R b/R/plot_perturb.R
index afcad85..b90fc08 100644
--- a/R/plot_perturb.R
+++ b/R/plot_perturb.R
@@ -19,9 +19,8 @@
#'
#' @importFrom fields image.plot
#' @importFrom graphics par
-#' @importFrom raster crs raster projectRaster
-#' @importFrom sp CRS
#' @importFrom spatstat.geom pixellate
+#' @importFrom terra crs project rast
#'
#' @export
#'
@@ -179,41 +178,45 @@ plot_perturb <- function(input,
if (predict == TRUE) {
- # Convert 'im' objects to spatially projected 'RasterLayer' objects
+ # Convert 'im' objects to spatially projected 'SpatRaster' objects
lrr_mean <- spatstat.geom::pixellate(input$predict,
weights = marks(input$predict)$lrr_mean)
- lrr_mean <- raster::raster(lrr_mean)
- raster::crs(lrr_mean) <- sp::CRS(SRS_string = cref0)
+ lrr_mean <- terra::rast(lrr_mean)
+ terra::crs(lrr_mean) <- cref0
if (!is.null(cref1)) {
- lrr_mean <- raster::projectRaster(lrr_mean,
- crs = sp::CRS(SRS_string = cref1))
+ lrr_mean <- terra::project(lrr_mean,
+ y = cref1,
+ method = "bilinear")
}
lrr_sd <- spatstat.geom::pixellate(input$predict,
weights = marks(input$predict)$lrr_sd)
- lrr_sd <- raster::raster(lrr_sd)
- raster::crs(lrr_sd) <- sp::CRS(SRS_string = cref0)
+ lrr_sd <- terra::rast(lrr_sd)
+ terra::crs(lrr_sd) <- cref0
if (!is.null(cref1)) {
- lrr_sd <- raster::projectRaster(lrr_sd,
- crs = sp::CRS(SRS_string = cref1))
+ lrr_sd <- terra::project(lrr_sd,
+ y = cref1,
+ method = "bilinear")
}
pval_mean <- spatstat.geom::pixellate(input$predict,
weights = marks(input$predict)$pval_mean)
- pval_mean <- raster::raster(pval_mean)
- raster::crs(pval_mean) <- sp::CRS(SRS_string = cref0)
+ pval_mean <- terra::rast(pval_mean)
+ terra::crs(pval_mean) <- cref0
if (!is.null(cref1)) {
- pval_mean <- raster::projectRaster(pval_mean,
- crs = sp::CRS(SRS_string = cref1))
+ pval_mean <- terra::project(pval_mean,
+ y = cref1,
+ method = "bilinear")
}
pval_prop <- spatstat.geom::pixellate(input$predict,
weights = marks(input$predict)$pval_prop)
- pval_prop <- raster::raster(pval_prop)
- raster::crs(pval_prop) <- sp::CRS(SRS_string = cref0)
+ pval_prop <- terra::rast(pval_prop)
+ terra::crs(pval_prop) <- cref0
if (!is.null(cref1)) {
- pval_prop <- raster::projectRaster(pval_prop,
- crs = sp::CRS(SRS_string = cref1))
+ pval_prop <- terra::project(pval_prop,
+ y = cref1,
+ method = "bilinear")
}
# Plot 5: mean log relative risk
diff --git a/R/plot_predict.R b/R/plot_predict.R
index e2ffa61..d9caa26 100644
--- a/R/plot_predict.R
+++ b/R/plot_predict.R
@@ -16,8 +16,7 @@
#'
#' @importFrom fields image.plot
#' @importFrom graphics par
-#' @importFrom raster crs cut image projectRaster raster reclassify values
-#' @importFrom sp coordinates CRS gridded
+#' @importFrom terra crs image project rast classify values
#' @export
#'
#' @examples
@@ -31,8 +30,8 @@
#' grad <- spatstat.data::bei.extra[[2]]
#' elev$v <- scale(elev)
#' grad$v <- scale(grad)
-#' elev_raster <- raster::raster(elev)
-#' grad_raster <- raster::raster(grad)
+#' elev_raster <- terra::rast(elev)
+#' grad_raster <- terra::rast(grad)
#'
#' # Presence data
#' presence <- spatstat.data::bei
@@ -57,14 +56,15 @@
#' obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
#'
#' # Prediction Data
-#' predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-#' predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+#' predict_xy <- terra::crds(elev_raster)
+#' predict_locs <- as.data.frame(predict_xy)
+#' predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+#' predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
#'
#' # Run lrren
#' test_lrren <- lrren(obs_locs = obs_locs,
#' predict_locs = predict_locs,
-#' predict = TRUE,
-#' cv = TRUE)
+#' predict = TRUE)
#'
#' # Run plot_predict
#' plot_predict(input = test_lrren, cref0 = "EPSG:5472")
@@ -90,59 +90,50 @@ plot_predict <- function(input,
op <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(op))
+
# Convert to geospatial rasters
- predict_risk <- data.frame("x" = input$out$predict$predict_locs.x,
- "y" = input$out$predict$predict_locs.y,
- "v" = input$out$predict$rr)
- naband <- predict_risk # save for next step
- sp::coordinates(predict_risk) <- ~ x + y # coordinates
- sp::gridded(predict_risk) <- TRUE # gridded
- predict_risk_raster <- raster::raster(predict_risk)
- raster::crs(predict_risk_raster) <- sp::CRS(SRS_string = cref0)
+ ## log relative risk
+ predict_risk <- data.frame("x" = input$out$predict[ , 1],
+ "y" = input$out$predict[ , 2],
+ "v" = input$out$predict$rr)
+ predict_risk_raster <- terra::rast(predict_risk)
+ terra::crs(predict_risk_raster) <- cref0
if (!is.null(cref1)) {
- predict_risk_raster <- raster::projectRaster(predict_risk_raster,
- crs = sp::CRS(SRS_string = cref1),
- method = "ngb",
- legacy = TRUE)
+ predict_risk_raster <- terra::project(predict_risk_raster,
+ y = cref1,
+ method = "bilinear")
+
}
-
- # Create separate layer for NAs (if any)
- naband$v <- ifelse(is.na(naband$v), 9999, naband$v)
- sp::coordinates(naband) <- ~ x + y # coordinates
- sp::gridded(naband) <- TRUE # gridded
- NA_risk_raster <- raster::raster(naband)
- raster::crs(NA_risk_raster) <- sp::CRS(SRS_string = cref0)
+
+ ## p-value
+ predict_tol <- data.frame("x" = input$out$predict[ , 1],
+ "y" = input$out$predict[ , 2],
+ "v" = input$out$predict$pval)
+ predict_tol_raster <- terra::rast(predict_tol)
+ terra::crs(predict_tol_raster) <- cref0
if (!is.null(cref1)) {
- NA_risk_raster <- raster::projectRaster(NA_risk_raster,
- crs = sp::CRS(SRS_string = cref1),
- method = "ngb",
- legacy = TRUE)
+ predict_tol_raster <- terra::project(predict_tol_raster,
+ y = cref1,
+ method = "bilinear")
}
- naband_reclass <- raster::reclassify(NA_risk_raster,
- c(-Inf, 9998, NA,
- 9998, Inf, 1))
- if (all(is.na(raster::values(naband_reclass)))) { naband_reclass <- NULL }
+ terra::values(predict_tol_raster) <- cut(terra::values(predict_tol_raster),
+ breaks = c(-Inf, alpha / 2, 1 - alpha / 2, Inf),
+ right = FALSE)
- # Convert to geospatial raster
- predict_tol <- data.frame("x" = input$out$predict$predict_locs.x,
- "y" = input$out$predict$predict_locs.y,
- "v" = input$out$predict$pval)
- sp::coordinates(predict_tol) <- ~ x + y # coordinates
- sp::gridded(predict_tol) <- TRUE # gridded
- predict_tol_raster <- raster::raster(predict_tol)
- raster::crs(predict_tol_raster) <- sp::CRS(SRS_string = cref0)
+ ## Separate layer for NAs (if any)
+ naband <- predict_risk
+ naband$v <- ifelse(is.na(naband$v), 9999, NA)
+ naband_raster <- terra::rast(naband)
+ terra::crs(naband_raster) <- cref0
if (!is.null(cref1)) {
- predict_tol_raster <- raster::projectRaster(predict_tol_raster,
- crs = sp::CRS(SRS_string = cref1),
- method = "ngb",
- legacy = TRUE)
+ naband_raster <- terra::project(naband_raster,
+ y = cref1,
+ method = "near")
}
-
- reclass_tol <- raster::cut(predict_tol_raster,
- breaks = c(-Inf, alpha / 2, 1 - alpha / 2, Inf),
- right = FALSE)
-
+
+ if (all(is.na(terra::values(naband_raster)))) { naband_raster <- NULL }
+
# Plot 1: log relative risk
rrp <- div_plot(input = predict_risk_raster,
cols = plot_cols[1:3],
@@ -164,12 +155,12 @@ plot_predict <- function(input,
las = 0,
labels = rrp$labels,
cex.axis = 0.67))
- if (!is.null(naband_reclass)) {
- raster::image(naband_reclass, col = plot_cols[4], add = TRUE)
+ if (!is.null(naband_raster)) {
+ terra::image(naband_raster, y = 1, col = plot_cols[4], add = TRUE)
}
# Plot 2: Significant p-values
- if (all(raster::values(reclass_tol)[!is.na(raster::values(reclass_tol))] == 2)) {
+ if (all(terra::values(predict_tol_raster)[!is.na(terra::values(predict_tol_raster))] == 2)) {
pcols <- plot_cols[2]
brp <- c(1, 3)
atp <- 2
@@ -181,7 +172,7 @@ plot_predict <- function(input,
labp <- c("presence", "insignificant", "absence")
}
- p2 <- fields::image.plot(reclass_tol,
+ p2 <- fields::image.plot(predict_tol_raster,
breaks = brp,
col = pcols,
axes = TRUE,
@@ -193,7 +184,7 @@ plot_predict <- function(input,
labels = labp,
las = 0,
cex.axis = 0.67))
- if (!is.null(naband_reclass)) {
- raster::image(naband_reclass, col = plot_cols[4], add = TRUE)
+ if (!is.null(naband_raster)) {
+ terra::image(naband_raster, y = 1, col = plot_cols[4], add = TRUE)
}
}
diff --git a/R/seq_plot.R b/R/seq_plot.R
index 0ef2be9..da739dc 100644
--- a/R/seq_plot.R
+++ b/R/seq_plot.R
@@ -1,8 +1,8 @@
-#' Prepare an 'im' or 'raster' object for plotting with sequential color palette
+#' Prepare an 'im' or 'SpatRaster' object for plotting with sequential color palette
#'
-#' Internal function to convert 'im' object or 'RasterLayer' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_perturb}} function.
+#' Internal function to convert 'im' object or 'SpatRaster' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_perturb}} function.
#'
-#' @param input An object of class 'im' or 'RasterLayer' from the \code{\link{perlrren}} function.
+#' @param input An object of class 'im' or 'SpatRaster' from the \code{\link{perlrren}} function.
#' @param plot_cols Character string of length three (3) specifying the colors for plotting: 1) presence, 2) neither, and 3) absence from the \code{\link{plot_obs}} function.
#' @param thresh_up Numeric. The upper value to concatenate the color key. The default (NULL) uses the maximum value from \code{input}.
#' @param digits Integer. The number of significant digits for the labels using the \code{\link[base]{round}} function (default is 1).
@@ -18,6 +18,7 @@
#' }
#'
#' @importFrom grDevices colorRampPalette
+#' @importFrom terra rast
#'
#' @keywords internal
#'
@@ -28,10 +29,10 @@ seq_plot <- function(input,
# Inputs
if (inherits(input, "im")) {
- out <- raster::raster(input)
+ out <- terra::rast(input)
} else { out <- input }
- max_raw_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of raster
+ max_raw_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of SpatRaster
# Restrict spurious standard deviation values
if (!is.null(thresh_up)) {
@@ -41,7 +42,7 @@ seq_plot <- function(input,
out[out >= thresh_up] <- thresh_up
}
- max_absolute_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of raster
+ max_absolute_value <- max(out[is.finite(out)], na.rm = TRUE) # maximum absolute value of SpatRaster
ncols <- length(out[!is.na(out)]) # number of values
## Colors
diff --git a/R/utils.R b/R/utils.R
index 3b10d06..ecf5807 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -32,3 +32,52 @@ progBar <- function(kk, N, per = 1) {
if (kk == N) cat("\r")
}
}
+
+# Temporary stopgap for converting 'SpatRaster' to 'im'
+## Modification of `as.im.RasterLayer` from `maptools` because `maptools` is retiring in 2023
+as.im.SpatRaster <- function (from, factor.col.name = NULL) {
+ if (!requireNamespace("spatstat.geom", quietly = TRUE)) {
+ stop(paste("package spatstat.geom required; please install it (or the full spatstat package) first"))
+ }
+ else {
+ spst_ver <- try(utils::packageVersion("spatstat"), silent = TRUE)
+ if (!inherits(spst_ver, "try-error") && spst_ver < 2 -
+ 0) {
+ stop(paste("You have an old version of spatstat installed which is",
+ " incompatible with spatstat.geom. Please update spatstat (or uninstall it).",
+ sep = ""))
+ }
+ }
+ if (!requireNamespace("sf", quietly = TRUE))
+ stop("package sf required for coercion")
+ if (!requireNamespace("terra", quietly = TRUE))
+ stop("package terra required for coercion")
+ if (!terra::hasValues(from))
+ stop("values required in SpatRaster object")
+ rs <- terra::res(from)
+ orig <- sf::st_bbox(from)[1:2] + 0.5 * rs
+ dm <- dim(from)[2:1]
+ xx <- unname(orig[1] + cumsum(c(0, rep(rs[1], dm[1] - 1))))
+ yy <- unname(orig[2] + cumsum(c(0, rep(rs[2], dm[2] - 1))))
+ val <- terra::values(from)
+ if (is.factor(from)) {
+ lev <- levels(from)[[1]]
+ if (!is.null(factor.col.name)) {
+ if (factor.col.name %in% colnames(lev)) {
+ factor.col <- which(colnames(lev) == factor.col.name)
+ }
+ else {
+ stop("'factor.col.name' is not a column name of the SpatRaster 'from'")
+ }
+ }
+ else {
+ factor.col <- length(lev)
+ }
+ val <- factor(val, levels = lev$ID, labels = lev[[factor.col]])
+ }
+ dim(val) <- dm
+ val <- spatstat.geom::transmat(val, from = list(x = "-i",
+ y = "j"), to = "spatstat")
+ im <- spatstat.geom::im(val, xcol = xx, yrow = yy)
+ return(im)
+}
diff --git a/README.md b/README.md
index e78d709..1c378d2 100644
--- a/README.md
+++ b/README.md
@@ -2,15 +2,18 @@ envi: Environmental Interpolation using Spatial Kernel Density Estimation ![]()
+[](https://github.com/lance-waller-lab/envi/actions/workflows/R-CMD-check.yaml)
[](https://cran.r-project.org/package=envi)
[](https://cran.r-project.org/package=envi)
-[](https://r-pkg.org/pkg/envi)
+[](https://r-pkg.org/pkg/envi)
+[](https://www.r-pkg.org:443/pkg/envi)
[](https://opensource.org/licenses/Apache-2.0)
[](https://doi.org/10.5281/zenodo.5347826)
+
-**Date repository last updated**: August 30, 2022
+**Date repository last updated**: December 16, 2022
@@ -106,10 +109,10 @@ set.seed(1234) # for reproducibility
# ------------------ #
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
# -------------- #
# Prepare inputs #
@@ -122,14 +125,14 @@ elev <- spatstat.data::bei.extra[[1]]
grad <- spatstat.data::bei.extra[[2]]
elev$v <- scale(elev)
grad$v <- scale(grad)
-elev_raster <- raster::raster(elev)
-grad_raster <- raster::raster(grad)
+elev_raster <- terra::rast(elev)
+grad_raster <- terra::rast(grad)
# Presence data
presence <- spatstat.data::bei
spatstat.geom::marks(presence) <- data.frame("presence" = rep(1, presence$n),
- "lon" = presence$x,
- "lat" = presence$y)
+ "lon" = presence$x,
+ "lat" = presence$y)
spatstat.geom::marks(presence)$elev <- elev[presence]
spatstat.geom::marks(presence)$grad <- grad[presence]
@@ -148,8 +151,10 @@ obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
-predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(elev_raster)
+predict_locs <- as.data.frame(predict_xy)
+predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# ----------- #
# Run lrren() #
@@ -180,7 +185,6 @@ envi::plot_predict(test1,
# ------------- #
envi::plot_cv(test1)
-
```
@@ -224,10 +228,10 @@ set.seed(1234) # for reproducibility
# ------------------ #
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
# -------------- #
# Prepare inputs #
@@ -269,7 +273,7 @@ spatstat.geom::marks(obs_locs) <- spatstat.geom::marks(obs_locs)[ , c(4, 2, 3, 1
test3 <- envi::perlrren(obs_ppp = obs_locs,
covariates = ims,
- radii = c(10,100,500),
+ radii = c(10, 100, 500),
verbose = FALSE, # may not be availabe if parallel = TRUE
parallel = TRUE,
n_sim = 100)
@@ -282,7 +286,6 @@ envi::plot_perturb(test3,
cref0 = "EPSG:5472",
cref1 = "EPSG:4326",
cov_labs = c("elev", "grad"))
-
```
diff --git a/cran-comments.md b/cran-comments.md
index a900b1b..ded53b7 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,9 +1,10 @@
-## This is the fifteenth resubmission
+## This is the sixteenth resubmission
* Updates since previous submission:
- * Uwe Ligges suggested (2022-08-26) that some parallel cluster is not cleanly closed in the tests
- * Added `future::plan(future::multisession)` in tests to remove the files in temp directory
- * Added links to `sparr` package in 'package.R'
+ * Migrated R-spatial dependencies
+ * Replaced `raster`, `rgeos`, and `sp` packages in Imports with `terra` and `sf` because of imminent package retirements
+ * Removed `maptools` from Suggests (replaced with new internal function `as.im.SpatRaster()`)
+ * Updated test, examples, vignette, and documentation throughout
* Documentation for "envi-package.Rd", "pval_correct.Rd", DESCRIPTION, and vignette references the following DOIs, which throws a NOTE but are valid URLs:
*
@@ -13,7 +14,7 @@
*
## Test environments
-* local OS X install, R 4.2.1
+* local Windows install, R 4.2.1
* win-builder, (devel, release, oldrelease)
* Rhub
* Fedora Linux, R-devel, clang, gfortran
diff --git a/inst/CITATION b/inst/CITATION
index 879771f..2c92361 100755
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -5,7 +5,7 @@ citEntry(entry = "manual",
author = personList(as.person("Ian D. Buller")),
publisher = "The Comprehensive R Archive Network",
year = "2022",
- number = "0.1.15",
+ number = "0.1.16.9000",
doi = "10.5281/zenodo.5347826",
url = "https://cran.r-project.org/package=envi",
@@ -13,7 +13,7 @@ citEntry(entry = "manual",
paste("Ian D. Buller (2022).",
"envi: Environmental Interpolation using Spatial Kernel Density Estimation.",
"The Comprehensive R Archive Network.",
- "v0.1.15.",
+ "v0.1.16.9000.",
"DOI:10.5281/zenodo.5347826",
"Accessed by: https://cran.r-project.org/package=envi")
)
diff --git a/man/div_plot.Rd b/man/div_plot.Rd
index da8c00b..877c2c8 100644
--- a/man/div_plot.Rd
+++ b/man/div_plot.Rd
@@ -2,7 +2,7 @@
% Please edit documentation in R/div_plot.R
\name{div_plot}
\alias{div_plot}
-\title{Prepare an 'im' or 'raster' object for plotting with diverging color palette}
+\title{Prepare an 'im' or 'SpatRaster' object for plotting with diverging color palette}
\usage{
div_plot(
input,
@@ -14,7 +14,7 @@ div_plot(
)
}
\arguments{
-\item{input}{An object of class 'im' or 'RasterLayer' from the \code{\link{lrren}} function.}
+\item{input}{An object of class 'im' or 'SpatRaster' from the \code{\link{lrren}} function.}
\item{midpoint}{Numeric. The value to center the diverging color palette.}
@@ -38,6 +38,6 @@ An object of class 'list'. This is a named list with the following components:
}
}
\description{
-Internal function to convert 'im' object or 'RasterLayer' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_obs}}, \code{\link{plot_predict}}, and \code{\link{plot_perturb}} functions.
+Internal function to convert 'im' object or 'SpatRaster' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_obs}}, \code{\link{plot_predict}}, and \code{\link{plot_perturb}} functions.
}
\keyword{internal}
diff --git a/man/envi-package.Rd b/man/envi-package.Rd
index a2cb387..06b8c13 100644
--- a/man/envi-package.Rd
+++ b/man/envi-package.Rd
@@ -34,7 +34,7 @@ Key content of the 'envi' package include:\cr
\code{\link{plot_perturb}} Visualizes the \code{\link{perlrren}} output, specifically four summary statistics of the iterations, including mean log relative risk, standard deviation of the log relative risk, mean p-value, and proportion of iterations the p-value was significant based on an alpha-level threshold. It also can predict the spatial distribution of the summary statistics.
}
\section{Dependencies}{
- The 'envi' package relies heavily upon \code{\link{sparr}}, \code{\link{spatstat.geom}}, and \code{\link{raster}}. For a single species (presence/absence data), the spatial relative risk function uses the \code{\link[sparr]{risk}} function. Cross-validation is can be performed in parallel using the \code{\link{future}}, \code{\link{doFuture}}, \code{\link{doRNG}}, and \code{\link{foreach}} packages. Spatial perturbation is performed using the \code{\link[spatstat.geom]{rjitter}} function. Basic visualizations rely on the \code{\link[spatstat.geom]{plot.ppp}} and \code{\link[fields]{image.plot}} functions.
+ The 'envi' package relies heavily upon \code{\link{sparr}}, \code{\link{spatstat.geom}}, \code{\link{sf}}, and \code{\link{terra}}. For a single species (presence/absence data), the spatial relative risk function uses the \code{\link[sparr]{risk}} function. Cross-validation is can be performed in parallel using the \code{\link{future}}, \code{\link{doFuture}}, \code{\link{doRNG}}, and \code{\link{foreach}} packages. Spatial perturbation is performed using the \code{\link[spatstat.geom]{rjitter}} function. Basic visualizations rely on the \code{\link[spatstat.geom]{plot.ppp}} and \code{\link[fields]{image.plot}} functions.
}
\author{
diff --git a/man/lrren.Rd b/man/lrren.Rd
index 678e7f6..ec7f942 100644
--- a/man/lrren.Rd
+++ b/man/lrren.Rd
@@ -113,8 +113,8 @@ if (interactive()) {
grad <- spatstat.data::bei.extra[[2]]
elev$v <- scale(elev)
grad$v <- scale(grad)
- elev_raster <- raster::raster(elev)
- grad_raster <- raster::raster(grad)
+ elev_raster <- terra::rast(elev)
+ grad_raster <- terra::rast(grad)
# Presence data
presence <- spatstat.data::bei
@@ -139,8 +139,10 @@ if (interactive()) {
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
- predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
- predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+ predict_xy <- terra::crds(elev_raster)
+ predict_locs <- as.data.frame(predict_xy)
+ predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+ predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Run lrren
test_lrren <- lrren(obs_locs = obs_locs,
diff --git a/man/plot_cv.Rd b/man/plot_cv.Rd
index a3ac764..6e9ed92 100644
--- a/man/plot_cv.Rd
+++ b/man/plot_cv.Rd
@@ -28,8 +28,8 @@ if (interactive()) {
grad <- spatstat.data::bei.extra[[2]]
elev$v <- scale(elev)
grad$v <- scale(grad)
- elev_raster <- raster::raster(elev)
- grad_raster <- raster::raster(grad)
+ elev_raster <- terra::rast(elev)
+ grad_raster <- terra::rast(grad)
# Presence data
presence <- spatstat.data::bei
@@ -53,14 +53,8 @@ if (interactive()) {
obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
-# Prediction Data
- predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
- predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
-
# Run lrren
test_lrren <- lrren(obs_locs = obs_locs,
- predict_locs = predict_locs,
- predict = TRUE,
cv = TRUE)
# Run plot_cv
diff --git a/man/plot_obs.Rd b/man/plot_obs.Rd
index ecbbf02..480c2ab 100644
--- a/man/plot_obs.Rd
+++ b/man/plot_obs.Rd
@@ -46,8 +46,8 @@ if (interactive()) {
grad <- spatstat.data::bei.extra[[2]]
elev$v <- scale(elev)
grad$v <- scale(grad)
- elev_raster <- raster::raster(elev)
- grad_raster <- raster::raster(grad)
+ elev_raster <- terra::rast(elev)
+ grad_raster <- terra::rast(grad)
# Presence data
presence <- spatstat.data::bei
@@ -71,15 +71,9 @@ if (interactive()) {
obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
-# Prediction Data
- predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
- predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
-
# Run lrren
test_lrren <- lrren(obs_locs = obs_locs,
- predict_locs = predict_locs,
- predict = TRUE,
- cv = TRUE)
+ cv = FALSE)
# Run plot_obs
plot_obs(input = test_lrren)
diff --git a/man/plot_predict.Rd b/man/plot_predict.Rd
index 8f3cad1..4b6cdb1 100644
--- a/man/plot_predict.Rd
+++ b/man/plot_predict.Rd
@@ -52,8 +52,8 @@ if (interactive()) {
grad <- spatstat.data::bei.extra[[2]]
elev$v <- scale(elev)
grad$v <- scale(grad)
- elev_raster <- raster::raster(elev)
- grad_raster <- raster::raster(grad)
+ elev_raster <- terra::rast(elev)
+ grad_raster <- terra::rast(grad)
# Presence data
presence <- spatstat.data::bei
@@ -78,14 +78,15 @@ if (interactive()) {
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
- predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
- predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+ predict_xy <- terra::crds(elev_raster)
+ predict_locs <- as.data.frame(predict_xy)
+ predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+ predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Run lrren
test_lrren <- lrren(obs_locs = obs_locs,
predict_locs = predict_locs,
- predict = TRUE,
- cv = TRUE)
+ predict = TRUE)
# Run plot_predict
plot_predict(input = test_lrren, cref0 = "EPSG:5472")
diff --git a/man/seq_plot.Rd b/man/seq_plot.Rd
index 8a44b64..2e3da9e 100644
--- a/man/seq_plot.Rd
+++ b/man/seq_plot.Rd
@@ -2,12 +2,12 @@
% Please edit documentation in R/seq_plot.R
\name{seq_plot}
\alias{seq_plot}
-\title{Prepare an 'im' or 'raster' object for plotting with sequential color palette}
+\title{Prepare an 'im' or 'SpatRaster' object for plotting with sequential color palette}
\usage{
seq_plot(input, cols, thresh_up = NULL, digits = 1)
}
\arguments{
-\item{input}{An object of class 'im' or 'RasterLayer' from the \code{\link{perlrren}} function.}
+\item{input}{An object of class 'im' or 'SpatRaster' from the \code{\link{perlrren}} function.}
\item{thresh_up}{Numeric. The upper value to concatenate the color key. The default (NULL) uses the maximum value from \code{input}.}
@@ -27,6 +27,6 @@ An object of class 'list'. This is a named list with the following components:
}
}
\description{
-Internal function to convert 'im' object or 'RasterLayer' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_perturb}} function.
+Internal function to convert 'im' object or 'SpatRaster' object to values readable by \code{\link[fields]{image.plot}} function within the \code{\link{plot_perturb}} function.
}
\keyword{internal}
diff --git a/tests/testthat/test-lrren.R b/tests/testthat/test-lrren.R
index c68a664..2872369 100644
--- a/tests/testthat/test-lrren.R
+++ b/tests/testthat/test-lrren.R
@@ -7,10 +7,11 @@ context("lrren")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
+library(sf)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
set.seed(1234)
# -------------- #
@@ -23,8 +24,8 @@ elev <- spatstat.data::bei.extra$elev
grad <- spatstat.data::bei.extra$grad
elev$v <- scale(elev)
grad$v <- scale(grad)
-elev_raster <- raster::raster(elev)
-grad_raster <- raster::raster(grad)
+elev_raster <- terra::rast(elev)
+grad_raster <- terra::rast(grad)
## Presence Locations
presence <- spatstat.data::bei
@@ -50,19 +51,18 @@ obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
-predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(elev_raster)
+predict_locs <- as.data.frame(predict_xy)
+colnames(predict_locs) <- c("lon", "lat")
+predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Test custom window
custom_chull <- grDevices::chull(x = obs_locs[ , 5], y = obs_locs[ , 6])
custom_chull_pts <- obs_locs[c(custom_chull, custom_chull[1]), 5:6]
-custom_chull_pts <- rbind(custom_chull_pts, custom_chull_pts[1, ])
-custom_chull_poly <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(custom_chull_pts)), 1)))
-#add small buffer around polygon to include boundary points
-custom_poly <- custom_chull_poly@polygons[[1]]@Polygons[[1]]@coords #extract coordinates of new polygon
-custom_owin <- spatstat.geom::owin(poly = list(x = rev(custom_poly[ , 1]),
- y = rev(custom_poly[ , 2])))
-
+custom_poly <- sf::st_polygon(list(as.matrix(custom_chull_pts)))
+custom_owin <- spatstat.geom::owin(poly = list(x = rev(sf::st_coordinates(custom_poly)[ , 1]),
+ y = rev(sf::st_coordinates(custom_poly)[ , 2])))
test_that("lrren throws error with invalid arguments", {
diff --git a/tests/testthat/test-perlrren.R b/tests/testthat/test-perlrren.R
index 755b4de..a18b6ea 100644
--- a/tests/testthat/test-perlrren.R
+++ b/tests/testthat/test-perlrren.R
@@ -7,7 +7,6 @@ context("perlrren")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
diff --git a/tests/testthat/test-plot_cv.R b/tests/testthat/test-plot_cv.R
index 7a3b45e..9b4fe80 100644
--- a/tests/testthat/test-plot_cv.R
+++ b/tests/testthat/test-plot_cv.R
@@ -7,10 +7,10 @@ context("plot_cv")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
set.seed(1234)
# -------------- #
@@ -23,8 +23,8 @@ elev <- spatstat.data::bei.extra$elev
grad <- spatstat.data::bei.extra$grad
elev$v <- scale(elev)
grad$v <- scale(grad)
-elev_raster <- raster::raster(elev)
-grad_raster <- raster::raster(grad)
+elev_raster <- terra::rast(elev)
+grad_raster <- terra::rast(grad)
## Presence Locations
presence <- spatstat.data::bei
@@ -50,8 +50,11 @@ obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
-predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(elev_raster)
+predict_locs <- as.data.frame(predict_xy)
+colnames(predict_locs) <- c("lon", "lat")
+predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Run lrren
test_lrren <- envi::lrren(obs_locs = obs_locs,
diff --git a/tests/testthat/test-plot_obs.R b/tests/testthat/test-plot_obs.R
index 63bb430..fac1293 100644
--- a/tests/testthat/test-plot_obs.R
+++ b/tests/testthat/test-plot_obs.R
@@ -7,10 +7,10 @@ context("plot_obs")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
set.seed(1234)
# -------------- #
@@ -23,8 +23,8 @@ elev <- spatstat.data::bei.extra$elev
grad <- spatstat.data::bei.extra$grad
elev$v <- scale(elev)
grad$v <- scale(grad)
-elev_raster <- raster::raster(elev)
-grad_raster <- raster::raster(grad)
+elev_raster <- terra::rast(elev)
+grad_raster <- terra::rast(grad)
## Presence Locations
presence <- spatstat.data::bei
@@ -50,8 +50,11 @@ obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
-predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(elev_raster)
+predict_locs <- as.data.frame(predict_xy)
+colnames(predict_locs) <- c("lon", "lat")
+predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Run lrren
test_lrren <- envi::lrren(obs_locs = obs_locs,
diff --git a/tests/testthat/test-plot_perturb.R b/tests/testthat/test-plot_perturb.R
index c2edf0c..b0ef9fe 100644
--- a/tests/testthat/test-plot_perturb.R
+++ b/tests/testthat/test-plot_perturb.R
@@ -7,7 +7,6 @@ context("plot_perturb")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
@@ -92,7 +91,7 @@ test_that("plot_obs works", {
)
# With spatial transformation
- expect_warning(
+ expect_silent(
plot_perturb(input = test_perlrren,
cref0 = "EPSG:5472",
cref1 = "EPSG:4326")
diff --git a/tests/testthat/test-plot_predict.R b/tests/testthat/test-plot_predict.R
index fa3b6f1..2f902c0 100644
--- a/tests/testthat/test-plot_predict.R
+++ b/tests/testthat/test-plot_predict.R
@@ -7,10 +7,10 @@ context("plot_predict")
# Generate testing data
## Environmental Covariates
library(envi)
-library(raster)
library(spatstat.data)
library(spatstat.geom)
library(spatstat.random)
+library(terra)
set.seed(1234)
# -------------- #
@@ -23,8 +23,8 @@ elev <- spatstat.data::bei.extra$elev
grad <- spatstat.data::bei.extra$grad
elev$v <- scale(elev)
grad$v <- scale(grad)
-elev_raster <- raster::raster(elev)
-grad_raster <- raster::raster(grad)
+elev_raster <- terra::rast(elev)
+grad_raster <- terra::rast(grad)
## Presence Locations
presence <- spatstat.data::bei
@@ -50,8 +50,11 @@ obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
# Prediction Data
-predict_locs <- data.frame(raster::rasterToPoints(elev_raster))
-predict_locs$layer2 <- raster::extract(grad_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(elev_raster)
+predict_locs <- as.data.frame(predict_xy)
+colnames(predict_locs) <- c("lon", "lat")
+predict_locs$elev <- terra::extract(elev_raster, predict_xy)[ , 1]
+predict_locs$grad <- terra::extract(grad_raster, predict_xy)[ , 1]
# Run lrren
test_lrren <- envi::lrren(obs_locs = obs_locs,
@@ -125,7 +128,7 @@ test_that("plot_predict works", {
skip_on_cran()
# Full run
- expect_warning(
+ expect_silent(
plot_predict(input = test_lrren,
plot_cols = c("#8b3a3a", "#cccccc", "#0000cd", "#ffff00"),
alpha = 0.05,
@@ -143,7 +146,7 @@ test_that("plot_predict works", {
)
# With spatial transformation
- expect_warning(
+ expect_silent(
plot_predict(input = test_lrren,
plot_cols = c("#8b3a3a", "#cccccc", "#0000cd", "#ffff00"),
alpha = 0.05,
diff --git a/vignettes/vignette.Rmd b/vignettes/vignette.Rmd
index 151a7ef..ac812f0 100644
--- a/vignettes/vignette.Rmd
+++ b/vignettes/vignette.Rmd
@@ -17,7 +17,7 @@ knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = FAL
Start with the necessary packages and seed for the vignette.
```{r packages}
-loadedPackages <- c("envi", "maptools", "raster", "RStoolbox", "spatstat.data", "spatstat.geom", "spatstat.random")
+loadedPackages <- c("envi", "RStoolbox", "sf", "spatstat.data", "spatstat.geom", "spatstat.random", "terra")
invisible(lapply(loadedPackages, library, character.only = TRUE))
set.seed(1234) # for reproducibility
```
@@ -43,11 +43,11 @@ slopeangle$v <- scale(slopeangle)
waterdist$v <- scale(waterdist)
```
-Convert the covariate data to class `RasterLayer`.
+Convert the covariate data to class `SpatRaster`.
```{r rasterize}
-slopeangle_raster <- raster::raster(slopeangle)
-waterdist_raster <- raster::raster(waterdist)
+slopeangle_raster <- terra::rast(slopeangle)
+waterdist_raster <- terra::rast(waterdist)
```
Add appropriate marks to the `gorillas` data from `spatstat.data` package. These points are considered our "presence" locations.
@@ -103,8 +103,10 @@ obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
Extract the prediction locations within the study area from one of the covariates.
```{r predict_locs}
-predict_locs <- data.frame(raster::rasterToPoints(slopeangle_raster))
-predict_locs$layer2 <- raster::extract(waterdist_raster, predict_locs[, 1:2])
+predict_xy <- terra::crds(slopeangle_raster)
+predict_locs <- as.data.frame(predict_xy)
+predict_locs$slopeangle <- terra::extract(slopeangle_raster, predict_xy)[ , 1]
+predict_locs$waterdist <- terra::extract(waterdist_raster, predict_xy)[ , 1]
```
```{r cleanup2, echo = FALSE}
@@ -191,33 +193,24 @@ slopeangle$v <- scale(slopeangle$v)
waterdist$v <- scale(waterdist$v)
# Create rasters
-aspect <- raster::raster(aspect)
-elevation <- raster::raster(elevation)
-heat <- raster::raster(heat)
-slopeangle <- raster::raster(slopeangle)
-slopetype <- raster::raster(slopetype)
-vegetation <- raster::raster(vegetation)
-waterdist <- raster::raster(waterdist)
+aspect <- terra::rast(aspect)
+elevation <- terra::rast(elevation)
+heat <- terra::rast(heat)
+slopeangle <- terra::rast(slopeangle)
+slopetype <- terra::rast(slopetype)
+vegetation <- terra::rast(vegetation)
+waterdist <- terra::rast(waterdist)
# Reorder aspect to center by "N" instead of "S"
-raster::values(aspect) <- factor(values(aspect),
- levels = c("5", "6", "7", "8", "1", "2", "3", "4"))
+terra::values(aspect) <- factor(terra::values(aspect),
+ levels = c("5", "6", "7", "8", "1", "2", "3", "4"))
# Reorder slope types to order flatter types next to each other
-raster::values(slopetype) <- factor(values(slopetype),
- levels = c("1", "6", "3", "2", "4", "5"))
+terra::values(slopetype) <- factor(terra::values(slopetype),
+ levels = c("1", "6", "3", "2", "4", "5"))
-# Save names for RasterLayers
-names(aspect) <- "aspect"
-names(elevation) <- "elevation"
-names(heat) <- "heat"
-names(slopeangle) <- "slopeangle"
-names(slopetype) <- "slopetype"
-names(vegetation) <- "vegetation"
-names(waterdist) <- "waterdist"
-
-# RasterStack
-park <- raster::stack(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist)
+# Stack of SpatRasters
+park <- c(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist)
# Principal Component Analysis
pca_park <- RStoolbox::rasterPCA(park)
@@ -228,27 +221,27 @@ summary(pca_park$model) # PCA components
rm(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist, park) # conserve memory
```
-The first two components of the PCA explain over 57% of the variation across the seven covariates.
+The first two components of the PCA explain almost 64% of the variation across the seven covariates.
```{r loadings}
pca_park$model$loadings # PCA loadings
```
-The loadings of the first component is almost entirely the centered 'aspect' variable. The loadings of the second component is a combination of 'vegetation type,' centered 'slope type,' centered 'elevation,' and centered 'distance from water.'
+The loadings of the first component is almost entirely the centered 'aspect' variable. The loadings of the second component is a combination of centered 'elevation,' 'vegetation type,' centered 'distance from water,' and centered 'slope type.'
```{r bands}
# Extract Bands from PCA
pca_bands <- pca_park$map
pc1 <- pca_bands[[1]] # PC1
pc2 <- pca_bands[[2]] # PC2
-pc1 <- spatstat.geom::as.im(pc1) # convert to class 'im'
-pc2 <- spatstat.geom::as.im(pc2) # convert to class 'im'
+pc1 <- envi:::as.im.SpatRaster(terra::rast(pc1)) # convert to class 'im'
+pc2 <- envi:::as.im.SpatRaster(terra::rast(pc2)) # convert to class 'im'
spatstat.geom::plot.im(pc1,
main = 'Principal Component 1\nprimarily aspect (centered at "North")',
ann = TRUE,
axes = TRUE)
spatstat.geom::plot.im(pc2,
- main = 'Principal Component 2\ncombination of vegetation type, slope type,\nelevation, and distance from water',
+ main = 'Principal Component 2\ncombination of elevation, vegetation type,\ndistance from water, and slope type',
ann = TRUE,
axes = TRUE)
```
@@ -257,7 +250,7 @@ spatstat.geom::plot.im(pc2,
rm(pca_park) # conserve memory
```
-For the first component (centered aspect), we can observe northern (i.e., North, North East, North West) aspects are primarily located in the northern section of the park. For the second component (combination), we can see a large cluster of middle loading values (between -2 and 0) in the center of the park surrounded by high loading values (> 0) and a few smaller clusters of low loading values (< -2) within the large middle loading value cluster.
+For the first component (centered aspect), we can observe northern (i.e., North, North East, North West) aspects are primarily located in the northern section of the park. For the second component (combination), we can see a large cluster of middle loading values (between 0 and 2) in the center of the park surrounded by lower loading values (< 0) and a few smaller clusters of high loading values (> 2) and low loading values (< -2) within the large middle loading value cluster.
We prepare inputs for a new `lrren` run similar to the first example above. For consistency, we also reset the random number generator to select similar pseudo-absence locations as the example above.
@@ -283,8 +276,10 @@ obs_locs <- spatstat.geom::marks(obs_locs)
obs_locs$id <- seq(1, nrow(obs_locs), 1)
obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
-predict_locs <- data.frame(raster::rasterToPoints(pca_bands[[1]]))
-predict_locs$layer2 <- raster::extract(pca_bands[[2]], predict_locs[, 1:2])
+predict_xy <- terra::crds(terra::rast(pca_bands[[1]]))
+predict_locs <- as.data.frame(predict_xy)
+predict_locs$pc1 <- terra::extract(terra::rast(pca_bands[[1]]), predict_xy)[ , 1]
+predict_locs$pc2 <- terra::extract(terra::rast(pca_bands[[2]]), predict_xy)[ , 1]
```
```{r cleanup6, echo = FALSE}
@@ -317,7 +312,9 @@ plot_obs(pca_lrren,
We display the estimated ecological niche predicted to the study area within geographic space using the `plot_predict` function. We use the default two-tailed alpha level (`alpha = 0.05`), and we limit the color key for the log relative risk prediction to (-1, 1).
```{r pca_plot_predict}
-plot_predict(pca_lrren, cref0 = "EPSG:32632", cref1 = "EPSG:4326",
+plot_predict(pca_lrren,
+ cref0 = "EPSG:32632",
+ cref1 = "EPSG:4326",
lower_lrr = -1,
upper_lrr = 1)
```
@@ -332,9 +329,9 @@ plot_cv(pca_lrren)
rm(pca_lrren, obs_locs, predict_locs) # conserve memory
```
-Based on only the first two components of the Principal Component Analysis of seven covariates, we detected one ecological niche of gorillas compared to pseudo-absence points drawn randomly from within the park. Presence and absences appear separated more by Principal Component 2 (loading values between -4 and -1) than Principal Component 1. This is reflected in geographic space where the gorilla niche is located in an area similar to the large cluster of middle loading values of Principal Component 2, which is located in the central, northwestern section of the park. The log relative risk estimate accurately predicted about 73% of the gorilla nesting sites, which is markedly improved from our first example above (60%).
+Based on only the first two components of the Principal Component Analysis of seven covariates, we detected one ecological niche of gorillas compared to pseudo-absence points drawn randomly from within the park. Presence and absences appear separated more by Principal Component 2 (loading values between 1 and 4) than Principal Component 1. This is reflected in geographic space where the gorilla niche is located in an area similar to the large cluster of middle loading values of Principal Component 2, which is located in the central, northwestern section of the park. The log relative risk estimate accurately predicted about 75% of the gorilla nesting sites, which is markedly improved from our first example above (60%).
-The two components captured about 57% of the variation across the seven covariates, which may not be important in separating presence locations from absence locations within covariate space. Further studies can assess different combinations of principal components (e.g., Principal Component 2 and Principal Component 3) or include another yet unavailable covariate that may be ecologically important for gorillas.
+The two components captured about 64% of the variation across the seven covariates, which may not be important in separating presence locations from absence locations within covariate space. Further studies can assess different combinations of principal components (e.g., Principal Component 2 and Principal Component 3) or include another yet unavailable covariate that may be ecologically important for gorillas.
The original study by [Funwi-Gabga and Mateu (2012)](https://doi.org/10.1007/s00477-011-0541-1) used kernel density estimation techniques in geographic space (see: Figure 5a in reference). Our method applies a kernel density estimation technique in covariate space and shows results similar to the inhomogeneous spatial point process model results of the original study (see: Figure 10a in reference), which incorporates the same covariate information to predict the spatial distribution of gorilla nesting sites. Our approach predicted more nesting sites in the western section of the park than in the original study.
diff --git a/vignettes/vignette.html b/vignettes/vignette.html
index 2c97b30..ac05201 100644
--- a/vignettes/vignette.html
+++ b/vignettes/vignette.html
@@ -12,7 +12,7 @@
-
+
envi: Environmental Interpolation using Spatial Kernel Density Estimation
@@ -31,23 +31,23 @@
+code{white-space: pre-wrap;}
+span.smallcaps{font-variant: small-caps;}
+span.underline{text-decoration: underline;}
+div.column{display: inline-block; vertical-align: top; width: 50%;}
+div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
+ul.task-list{list-style: none;}
+