Merge pull request #20 from adamlilith/solstice_2022_2023

Solstice 2022 2023

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: enmSdmX
 Type: Package
 Title: Species Distribution Modeling and Ecological Niche Modeling
-Version: 1.0.0
-Date: 2022-12-05
+Version: 1.0.1
+Date: 2023-01-25
 Authors@R: 
 	c(
 		person(

NEWS.md

@@ -1,3 +1,10 @@
+enmSdmX 1.0.1
+===========
+o Fixed bug experienced by some users using predictEnmSdm() and predictMaxNet() (Thank you, Nikki!)
+o Fixed bug in trainByCrossValid() using improper call to evalContBoyce()
+o Fixed extract() bug in some examples
+o summarizeByCrossValid() now summarizes natural spline (NS) models
+
 enmSdmX 1.0.0
 ===========
 o First release on CRAN

R/geoFold.r

@@ -22,7 +22,7 @@
 #'
 #' @return A vector of integers the same length as the number of points in \code{x}. Each integer indicates which fold a point in \code{x} belongs to.
 #'
-#' @example man/examples/pointGeoFold_examples.r
+#' @example man/examples/geoFold_examples.r
 #' 
 #' @export
 geoFold <- function(

R/predictEnmSdm.r

@@ -138,7 +138,8 @@ predictEnmSdm <- function(
 				notNas <- which(stats::complete.cases(nd))
 				nd <- nd[notNas, , drop=FALSE]
 				preds <- gbm::predict.gbm(model, newdata, n.trees=model$gbm.call$n.trees, type='response', ...)
-				out <- newdata[[1L]] * NA
+				out <- newdata[[1L]]
+				out[] <- NA
 				out <- setValueByCell(out, preds, cell=notNas, format='raster')
 			} else {
 				out <- gbm::predict.gbm(model, newdata, n.trees=model$gbm.call$n.trees, type='response', ...)
@@ -179,7 +180,8 @@ predictEnmSdm <- function(
 			preds <- preds[ , '1']
 
 			if (inherits(newdata, 'SpatRaster')) {
-				out <- newdata[[1L]] * NA
+				out <- newdata[[1L]]
+				out[] <- NA
 				out <- setValueByCell(out, preds, cell=notNas, format='raster')
 			} else {
 				out <- rep(NA, nrow(newdata))

R/predictMaxNet.r

@@ -42,7 +42,8 @@ predictMaxNet <- function(
 		preds <- do.call(predictMaxNet, args = list(model = model, newdata = nd, clamp = clamp, type = type, ...))
 
 		if (inherits(newdata, 'SpatRaster')) {
-			out <- newdata[[1L]] * NA
+			out <- newdata[[1L]]
+			out[] <- NA
 			out <- setValueByCell(out, preds, cell=notNas, format='raster')
 		} else {
 			out <- rep(NA, nrow(newdata))

R/summaryByCrossValid.r

@@ -170,12 +170,24 @@ summaryByCrossValid <- function(
 
 			thisTuning <- tuning[[k]]
 
-			thisTerm <- thisTuning$model[1]
-			thisTerm <- strsplit(thisTerm, ' ')[[1]]
-			thisTerm <- thisTerm[3:length(thisTerm)]
-			if (any(thisTerm == '1')) thisTerm <- thisTerm[-which(thisTerm == '1')]
-			if (any(thisTerm == '+')) thisTerm <- thisTerm[-which(thisTerm == '+')]
-			if (any(thisTerm == '-')) thisTerm <- thisTerm[-which(thisTerm == '-')]
+			thisTerm <- thisTuning$model[1L]
+
+			# natural splines model
+			if (grepl(thisTerm, pattern='splines')) {
+
+				thisTerm <- strsplit(thisTerm, ' \\+\\ ')[[1L]]
+				if (any(thisTerm == '1')) thisTerm <- thisTerm[-which(thisTerm == '1')]
+
+			# "normal" GLM model
+			} else  {
+
+				thisTerm <- strsplit(thisTerm, ' ')[[1]]
+				thisTerm <- thisTerm[3:length(thisTerm)]
+				if (any(thisTerm == '1')) thisTerm <- thisTerm[-which(thisTerm == '1')]
+				if (any(thisTerm == '+')) thisTerm <- thisTerm[-which(thisTerm == '+')]
+				if (any(thisTerm == '-')) thisTerm <- thisTerm[-which(thisTerm == '-')]
+
+			}
 
 			for (countTerm in seq_along(thisTerm)) {
 				whichTerm <- which(out$term == thisTerm[countTerm])

R/trainByCrossValid.r

@@ -264,11 +264,11 @@ trainByCrossValid <- function(
 					kTuning$logLossDelta[countModel] <- metricTrain - metricTest
 
 					# CBI
-					# metricTrain <- evalContBoyce(pres=predToTrainPres, bg=predToTrainContrast, presWeight=trainPresWeights, contrastWeight=trainContrastWeights, na.rm=na.rm, ...)
-					# metricTest <- evalContBoyce(pres=predToTestPres, bg=predToTestContrast, presWeight=testPresWeights, contrastWeight=testContrastWeights, na.rm=na.rm, ...)
+					# metricTrain <- evalContBoyce(pres=predToTrainPres, contrast=predToTrainContrast, presWeight=trainPresWeights, contrastWeight=trainContrastWeights, na.rm=na.rm, ...)
+					# metricTest <- evalContBoyce(pres=predToTestPres, contrast=predToTestContrast, presWeight=testPresWeights, contrastWeight=testContrastWeights, na.rm=na.rm, ...)
 
-					metricTrain <- evalContBoyce(pres=predToTrainPres, bg=predToTrainContrast, presWeight=trainPresWeights, contrastWeight=trainContrastWeights, na.rm=na.rm)
-					metricTest <- evalContBoyce(pres=predToTestPres, bg=predToTestContrast, presWeight=testPresWeights, contrastWeight=testContrastWeights, na.rm=na.rm)
+					metricTrain <- evalContBoyce(pres=predToTrainPres, contrast=predToTrainContrast, presWeight=trainPresWeights, contrastWeight=trainContrastWeights, na.rm=na.rm)
+					metricTest <- evalContBoyce(pres=predToTestPres, contrast=predToTestContrast, presWeight=testPresWeights, contrastWeight=testContrastWeights, na.rm=na.rm)
 
 					if (countModel == 1) kTuning$cbiDelta <- kTuning$cbiTest <- kTuning$cbiTrain <- NA
 					kTuning$cbiTrain[countModel] <- metricTrain

man/examples/bioticVelocity_examples.r

@@ -36,15 +36,20 @@ madClim2050 <- project(madClim2050, madAlbers)
 madClim2070 <- project(madClim2070, madAlbers)
 madClim2090 <- project(madClim2090, madAlbers)
 
+# Coordinate reference systems:
+wgs84 <- getCRS('WGS84') # WGS84
+madAlbers <- getCRS('madAlbers') # Madagascar Albers
+
 # lemur occurrence data
 data(lemurs)
-wgs84 <- getCRS('WGS84')
 occs <- lemurs[lemurs$species == 'Eulemur fulvus', ]
-occs <- st_as_sf(occs, coords = c('longitude', 'latitude'), crs = wgs84)
-occs <- st_transform(occs, getCRS('madAlbers'))
+occs <- vect(occs, geom=c('longitude', 'latitude'), crs=wgs84)
+occs <- project(occs, madAlbers)
 
+# eliminate cell duplicates
 occs <- elimCellDuplicates(occs, madClim)
 
+# extract environment at occurrences
 occEnv <- extract(madClim, occs, ID = FALSE)
 occEnv <- occEnv[complete.cases(occEnv), ]
 

man/examples/nearestEnvPoints_example.r

@@ -7,35 +7,47 @@
 library(sf)
 library(terra)
 
+# coordinate reference system
+wgs84 <- getCRS('WGS84')
+
+# lemur point data
 data(lemurs)
 precise <- lemurs[lemurs$species == 'Eulemur rubriventer', ]
 ll <- c('longitude', 'latitude')
-wgs84 <- getCRS('WGS84')
 precise <- sf::st_as_sf(precise[ , ll], coords=ll, crs=wgs84)
 
+# *fake* lemur administrative unit-level data
 faritras <- c('Vakinankaratra', 'Haute matsiatra', 'Ihorombe',
 'Vatovavy Fitovinany', 'Alaotra-Mangoro', 'Analanjirofo', 'Atsinanana',
 'Analamanga', 'Itasy')
+
 data(mad1)
 imprecise <- mad1[mad1$NAME_2 %in% faritras, ]
 
+# climate predictors
 rastFile <- system.file('extdata/madClim.tif', package='enmSdmX')
 rasts <- rast(rastFile)
 
 ### Plot environment of points and environments of each polygon closest to
 ### centroid of environments of points. In this example, we use the first two
 ### principal component axes to characterize the niche.
 
+# apply Nearest Environmental Point method
 envPtsPolys <- nearestEnvPoints(rasts, pts = precise, polys = imprecise,
 	pca = TRUE,	numPcs = 2)
 envPolys <- nearestEnvPoints(rasts, pts = precise, polys = imprecise, numPcs = 2,
 	out = 'polys')
 envPts <- nearestEnvPoints(rasts, pts = precise, polys = imprecise, numPcs = 2,
 	out = 'pts')
+
 allPolyEnvs <- extract(rasts, imprecise)
+
+# plot occurrences in environmental space
 plot(envPtsPolys$PC1, envPtsPolys$PC2, pch=16, col='black',
 	xlab='PC1', ylab='PC2')
+
 points(envPolys$PC1, envPolys$PC2, pch=21, bg='orange')
+
 legend(
 	'bottomleft',
 	inset = 0.01,
@@ -47,13 +59,15 @@ legend(
 
 ### compare identified environments to all environments across all polygons
 ###########################################################################
+
 env <- as.data.frame(rasts)
 pca <- stats::prcomp(env, center=TRUE, scale.=TRUE)
 
 allPolyEnvs <- extract(rasts, imprecise, ID = FALSE)
 allPolyEnvsPcs <- predict(pca, allPolyEnvs)
 allPolyEnvs <- cbind(allPolyEnvs, allPolyEnvsPcs)
 
+# plot in environmental space
 plot(allPolyEnvs$PC1, allPolyEnvs$PC2, pch=16, col='orange',
 	xlab='PC1', ylab='PC2')
 points(envPts$PC1, envPts$PC2, pch=16)
@@ -70,6 +84,7 @@ legend(
 ### display niches (minimum convex hulls) estimated
 ### using just precise or precise + imprecise records
 #####################################################
+
 pcs <- c('PC1', 'PC2')
 preciseIndices <- chull(envPts[ , pcs])
 preciseImpreciseIndices <- chull(envPtsPolys[ , pcs])
@@ -81,7 +96,7 @@ preciseImpreciseIndices <- c(preciseImpreciseIndices,
 preciseOnlyNiche <- envPts[preciseIndices, pcs]
 preciseImpreciseNiche <- envPtsPolys[preciseImpreciseIndices, pcs]
 
-# plot
+# plot in environmental space
 plot(allPolyEnvs$PC1, allPolyEnvs$PC2, pch=16, col='orange',
 	xlab='PC1', ylab='PC2')
 points(envPts$PC1, envPts$PC2, pch=16)

man/examples/nearestGeogPoints_example.r

@@ -5,34 +5,41 @@ library(terra)
 ### example using SpatVector inputs (terra package)
 ###################################################
 
-# Tananarive (Paris) / Laborde Grid - EPSG:29701
+# Get coordinate reference systems:
+# * WGS84
+# * Tananarive (Paris) / Laborde Grid - EPSG:29701
 wgs84 <- getCRS('WGS84')
 madProj <- getCRS('Madagascar Albers')
 
+# outline of Madagascar faritras
 data(mad1)
 mad1 <- vect(mad1)
 mad1 <- project(mad1, madProj)
 
+# lemur point data
 data(lemurs)
 redBelly <- lemurs[lemurs$species == 'Eulemur rubriventer', ]
 ll <- c('longitude', 'latitude')
-redBelly <- vect(redBelly[ , ll], geom=ll, crs=wgs84)
+redBelly <- vect(redBelly, geom=ll, crs=wgs84)
 redBelly <- project(redBelly, madProj)
 
+# *fake* lemur farita-level data
 faritras <- c('Vakinankaratra', 'Haute matsiatra', 'Ihorombe',
 'Vatovavy Fitovinany', 'Alaotra-Mangoro', 'Analanjirofo', 'Atsinanana',
 'Analamanga', 'Itasy')
 polys <- mad1[mad1$NAME_2 %in% faritras, ]
 
+# apply Nearest Geographic Point method
 mcpPolys <- nearestGeogPoints(polys = polys)
 mcpPts <- nearestGeogPoints(pts = redBelly, polys = NULL)
 mcpPolysPoints <- nearestGeogPoints(pts = redBelly, polys = polys)
 
-# extent of occurrence in m2
+# compare extent of occurrence (EOO) in m2
 expanse(mcpPolys)
 expanse(mcpPts)
 expanse(mcpPolysPoints)
 
+# plot minimum convex polygons
 plot(mad1, border='gray')
 plot(polys, col='gray80', add=TRUE)
 plot(mcpPolysPoints, col=scales::alpha('green', 0.4), add=TRUE)
@@ -52,34 +59,42 @@ border=c(NA, 'black', 'black', 'black', 'black'))
 ### example using sf* inputs (sf package)
 #########################################
 
-# Tananarive (Paris) / Laborde Grid - EPSG:29701
+# Get coordinate reference systems:
+# * WGS84
+# * Tananarive (Paris) / Laborde Grid - EPSG:29701
 madProj <- sf::st_crs(getCRS('Madagascar Albers'))
 wgs84 <- getCRS('WGS84')
 
+# outline of Madagascar faritras
 data(mad1)
 mad1 <- sf::st_transform(mad1, madProj)
 
+# lemur point occurrence data
 data(lemurs)
 redBelly <- lemurs[lemurs$species == 'Eulemur rubriventer', ]
 ll <- c('longitude', 'latitude')
 redBelly <- sf::st_as_sf(redBelly[ , ll], crs=wgs84, coords=ll)
 redBelly <- sf::st_transform(redBelly, madProj)
 
+# *fake* farita-level occurrences
 faritras <- c('Vakinankaratra', 'Haute matsiatra', 'Ihorombe',
 'Vatovavy Fitovinany', 'Alaotra-Mangoro', 'Analanjirofo', 'Atsinanana',
 'Analamanga', 'Itasy')
 polys <- mad1[mad1$NAME_2 %in% faritras, ]
 
+
+# apply Nearest Geographic Point method
 mcpPolys <- nearestGeogPoints(polys = polys, terra = FALSE)
 mcpPts <- nearestGeogPoints(pts = redBelly, polys = NULL, terra = FALSE)
 mcpPolysPoints <- nearestGeogPoints(pts = redBelly, polys = polys,
 terra = FALSE)
 
-# extent of occurrence in m2
+# extent of occurrence (EOO) in m2
 sf::st_area(mcpPolys)
 sf::st_area(mcpPts)
 sf::st_area(mcpPolysPoints)
 
+# plot minimum convex polygons
 plot(sf::st_geometry(mad1))
 plot(sf::st_geometry(polys), col='gray80', add=TRUE)
 plot(sf::st_geometry(mcpPolysPoints), col=scales::alpha('green', 0.4),

man/examples/trainByCrossValid_examples.r

@@ -21,16 +21,17 @@ set.seed(123)
 rastFile <- system.file('extdata/madClim.tif', package='enmSdmX')
 madClim <- rast(rastFile)
 
-crs <- sf::st_crs(madClim)
+# coordinate reference system
+wgs84 <- getCRS('WGS84')
 
 # lemur occurrence data
 data(lemurs)
 occs <- lemurs[lemurs$species == 'Eulemur fulvus', ]
-occs <- sf::st_as_sf(occs, coords=c('longitude', 'latitude'), crs=crs)
+occs <- vect(occs, geom=c('longitude', 'latitude'), crs=wgs84)
 
 occs <- elimCellDuplicates(occs, madClim)
 
-occEnv <- extract(madClim, occs, ID=FALSE)
+occEnv <- extract(madClim, occs, ID = FALSE)
 occEnv <- occEnv[complete.cases(occEnv), ]
 
 # create background sites (using just 1000 to speed things up!)
@@ -57,7 +58,7 @@ folds <- dismo::kfold(env, 3) # just 3 folds (for speed)
 ### calibrate models
 ####################
 
-cores <- 1 # increase this to go faster, if your computer can handle it
+cores <- 1 # increase this to go faster, if your computer handles it
 
 ## MaxEnt
 mxx <- trainByCrossValid(
@@ -149,7 +150,7 @@ brtx <- trainByCrossValid(
 	cores = cores
 )
 
-# summarize BRT parameters in best models
+# summarize BRT parameters across best models
 summaryByCrossValid(brtx)
 
 ## random forests