Move raster to suggests

17 files changed, 122 insertions, 76 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 86f9984..02a58d4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,10 +1,11 @@
 Package: AHMbook
 Type: Package
 Title: Functions and Data for the Book 'Applied Hierarchical Modeling in Ecology' Vols 1 and 2
-Version: 0.2.10
-Date: 2023-11-04
+Version: 0.2.11
+Date: 2024-01-14
 Depends: R (>= 2.10)
-Imports: grDevices, graphics, methods, stats, utils, raster, sp, fields, coda, unmarked (>= 0.12.2), mvtnorm, spdep
+Imports: grDevices, graphics, methods, stats, utils, sp, fields, coda, unmarked (>= 0.12.2), mvtnorm, spdep
+Suggests: raster
 Authors@R: c(
   person("Marc", "Kéry", role="aut"),
   person("Andy", "Royle", role="aut"),
diff --git a/NAMESPACE b/NAMESPACE
index 66273fb..5330d0e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,9 +16,6 @@ importFrom("utils", "head", "stack", "packageVersion", "data")
 
 importFrom("unmarked", "unmarkedFrameOccu", "unmarkedFramePCount", "unmarkedFrameDS",
           "occuRN", "bup", "ranef", "coef", "fitted", "residuals", "vcov", "confint")
-importFrom("raster", "raster", "extent", "extent<-", "rasterFromXYZ",
-          "getValues", "image", "addLayer",  "dropLayer", "ncell",
-          "values",  "values<-", "xyFromCell", "aggregate", "subset")
 # importFrom("plotrix", "draw.circle") # Now in file draw_circle.R
 # importFrom("RandomFields", "RFoptions", "RFsimulate", "RMexp")
 importFrom("fields",  "circulantEmbeddingSetup", "circulantEmbedding")
diff --git a/R/checkNamespace.R b/R/checkNamespace.R
new file mode 100644
index 0000000..386f046
--- /dev/null
+++ b/R/checkNamespace.R
@@ -0,0 +1,6 @@
+checkNamespace <- function(pkg){
+  if(!requireNamespace(pkg, quietly = TRUE)){
+    stop(paste("Install", pkg, "package to use this function", call.=FALSE))
+  }
+  TRUE
+}
diff --git a/R/simCJS_AHM2_3-2-2.R b/R/simCJS_AHM2_3-2-2.R
index e379b74..a22c417 100644
--- a/R/simCJS_AHM2_3-2-2.R
+++ b/R/simCJS_AHM2_3-2-2.R
@@ -22,6 +22,9 @@ simCJS <- function(
   #  ... and catch errors.
   #  Modified code for generating z and ch (gives different values with set.seed)
   #  Restore ask and par on exit.
+  
+  # Load raster package-------------------------------------
+  checkNamespace("raster")
 
   # Check and fixes for input data -------------------------
   n.occ <- round(n.occ[1])
@@ -86,12 +89,12 @@ simCJS <- function(
       par(mfrow = c(2, 2))
       # Plot the true alive/dead pattern (z)
       mapPalette <- colorRampPalette(c("white", "black"))
-      image(x = 1:n.occ, y = 1:n.ind, z = t(z), col = mapPalette(10), axes = TRUE,
+      raster::image(x = 1:n.occ, y = 1:n.ind, z = t(z), col = mapPalette(10), axes = TRUE,
         xlab = "Year", ylab = "Individual",
         main = 'z matrix of latent states in the CJS model: \nAlive (black) or dead (white) per individual and occasion')
 
       # Plot the observed alive/dead pattern (y, or ch)
-      image(x = 1:n.occ, y = 1:n.ind, z = t(ch), col = mapPalette(10), axes = TRUE,
+      raster::image(x = 1:n.occ, y = 1:n.ind, z = t(ch), col = mapPalette(10), axes = TRUE,
         xlab = "Year", ylab = "Individual",
         main = 'Observed data = capture-history matrix ch in the CJS model: \nDetected (black) or not detected (white) per individual and occasion')
       box()
@@ -100,7 +103,7 @@ simCJS <- function(
       tmp <- z    # copy z into tmp
       tmp[z==1 & ch == 0] <- -1.1  # Mark detection errors as -1
       mapPalette <- colorRampPalette(c("blue", "white", "black"))
-      image(x = 1:n.occ, y = 1:n.ind, z = t(tmp), col = mapPalette(10), axes = TRUE,
+      raster::image(x = 1:n.occ, y = 1:n.ind, z = t(tmp), col = mapPalette(10), axes = TRUE,
       xlab = "Year", ylab = "Individual", main = 'Combopic of z and ch: not in study (white), alive & detected (black), \nalive & undetected (blue) and dead (grey) per individual and occasion')
 
       # Population size trajectory of marked and alive in study area
diff --git a/R/simComm_AHM1_11-2_Simulate_community_data.R b/R/simComm_AHM1_11-2_Simulate_community_data.R
index e10bebb..619f95a 100644
--- a/R/simComm_AHM1_11-2_Simulate_community_data.R
+++ b/R/simComm_AHM1_11-2_Simulate_community_data.R
@@ -102,6 +102,9 @@ simComm <- function(type=c("det/nondet", "counts"), nsites=30, nreps=3, nspecies
 
 if(FALSE) x <- NULL # A kludge to cope with 'curve's odd way of using 'x'
 
+# Load raster package------------------------------------------
+checkNamespace("raster")
+
 # Checks and fixes for input data -----------------------------
 nsites <- round(nsites[1])
 nreps <- round(nreps[1])
@@ -221,19 +224,19 @@ if(type=="det/nondet"){
       mapPalette2 <- colorRampPalette(c("white", "yellow", "orange", "red"))
       # (3) True presence/absence matrix (z) for all species
       # mapPalette was 2 before
-      image(x = 1:nspecies, y = 1:nsites, z = t(z), col = mapPalette1(4), main =
+      raster::image(x = 1:nspecies, y = 1:nsites, z = t(z), col = mapPalette1(4), main =
       paste("True presence/absence (z) matrix\n (finite-sample N species =",
       Ntotal.fs,")"), frame = TRUE, xlim = c(0, nspecies+1),
       ylim = c(0, nsites+1), xlab = "Species", ylab = "Sites")
 
       # (4) Observed detection frequency for all species
-      image(x = 1:nspecies, y = 1:nsites, z = t(y.sum.all), col = mapPalette2(100),
+      raster::image(x = 1:nspecies, y = 1:nsites, z = t(y.sum.all), col = mapPalette2(100),
         main = paste("Observed detection frequencies"),
         xlim = c(0, nspecies+1), ylim = c(0, nsites+1),
         frame = TRUE, xlab = "Species", ylab = "Sites")
 
       # (5) Sites where a species was missed
-      image(x = 1:nspecies, y = 1:nsites, z = t(missed.sites), col = mapPalette1(2),
+      raster::image(x = 1:nspecies, y = 1:nsites, z = t(missed.sites), col = mapPalette1(2),
         main = paste("Matrix of missed presences\n (obs. N species =", Ntotal.obs,")"),
         frame = TRUE, xlim = c(0, nspecies+1), ylim = c(0, nsites+1), xlab = "Species",
         ylab = "Sites")
@@ -353,16 +356,16 @@ if(type=="counts"){
       mapPalette <- colorRampPalette(c("yellow", "orange", "red"))
       # (3) True abundance matrix (log(N+1)) for all species
       # mapPalette was 2 before
-      image(x = 1:nspecies, y = 1:nsites, z = log10(t(N)+1), col = mapPalette(100),
+      raster::image(x = 1:nspecies, y = 1:nsites, z = log10(t(N)+1), col = mapPalette(100),
        main =   paste("True log(abundance) (log10(N)) matrix\n (finite-sample N species =", sum(occurring.in.sample),")"), frame = TRUE, xlim = c(0, nspecies+1),
       zlim = c(0, log10(max(N))), xlab = "Species", ylab = "Sites")
       # (4) Observed maximum counts for all species
-      image(x = 1:nspecies, y = 1:nsites, z = log10(t(ymax.obs)+1), col = mapPalette(100), main = paste("Observed maximum counts (log10 + 1)"),
+      raster::image(x = 1:nspecies, y = 1:nsites, z = log10(t(ymax.obs)+1), col = mapPalette(100), main = paste("Observed maximum counts (log10 + 1)"),
         xlim = c(0, nspecies+1), frame = TRUE, xlab = "Species", ylab = "Sites", zlim = c(0, log10(max(N))))
       # (5) Ratio of max count to true N
       ratio <- ymax.obs/N
       ratio[ratio == "NaN"] <- 1
-      image(x = 1:nspecies, y = 1:nsites, z = t(ratio), col = mapPalette(100),
+      raster::image(x = 1:nspecies, y = 1:nsites, z = t(ratio), col = mapPalette(100),
         main = paste("Ratio of max count to true abundance (N)"),
         xlim = c(0, nspecies+1), frame = TRUE, xlab = "Species", ylab = "Sites",
         zlim = c(0, 1))
diff --git a/R/simDSM.R b/R/simDSM.R
index 494a1f8..70d7e60 100644
--- a/R/simDSM.R
+++ b/R/simDSM.R
@@ -13,6 +13,9 @@
 
 simDSM <- function(X, Ntotal = 400, sigma = 0.65, beta1 = 1.0,
     nsurveys = 2, xlim = c(-0.5, 3.5), ylim = c(-0.5, 4.5), show.plots = TRUE) {
+  
+  # Load raster package
+  checkNamespace("raster")
 
   # Create coordinates rasterized transect
   delta <- 0.2               # 2D bin width
@@ -68,7 +71,7 @@ simDSM <- function(X, Ntotal = 400, sigma = 0.65, beta1 = 1.0,
   if(show.plots) {
     oldpar <- par(mar = c(3,3,3,6)); on.exit(par(oldpar))
     tryPlot <- try( {
-      image(r <- rasterFromXYZ(cbind(gr,x)), col = topo.colors(10))
+      raster::image(r <- raster::rasterFromXYZ(cbind(gr,x)), col = topo.colors(10))
       image_scale(x, col = topo.colors(10))
       lines(X, col = "black", pch = 20, lwd = 3)
       points(sx, sy, pch = 16, col = "black", lwd=1 )
diff --git a/R/simDataDK1_21--.R b/R/simDataDK1_21--.R
index 63348a0..6ebd96e 100644
--- a/R/simDataDK1_21--.R
+++ b/R/simDataDK1_21--.R
@@ -62,6 +62,9 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   # Note that the parameters beta must be chosen such to avoid a too high 'filling' of the discrete approximation of the entire field B
   # ------------------------------------------------------
 
+  # -------------- Load raster package -----------------------
+  checkNamespace("raster")
+
   # -------------- Check and fix input -----------------------
   sqrt.npix <- round(sqrt.npix[1])
   stopifnotLength(alpha, 2)
@@ -83,9 +86,9 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   s.area <- 4
 
   # Approximate the landscape by many small pixels in a raster object
-  s <- temp <- raster(ncol=sqrt.npix, nrow=sqrt.npix, xmn=-1, xmx=1, ymn=-1, ymx=1, crs=NULL)
+  s <- temp <- raster::raster(ncol=sqrt.npix, nrow=sqrt.npix, xmn=-1, xmx=1, ymn=-1, ymx=1, crs=NULL)
   # s will become a Raster Stack, temp is a template to create new layers
-  s.loc <- xyFromCell(temp, 1:ncell(s))    # Coordinates of every pixel in S
+  s.loc <- raster::xyFromCell(temp, 1:raster::ncell(s))    # Coordinates of every pixel in S
 
   # PART A. Ecological process - where are the animals?
   # ---------------------------------------------------
@@ -95,26 +98,26 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   xcov <- standardize(xcov)    # Standardize covariate x
 
   # Fill covariate x into the raster s
-  values(s) <- xcov
+  raster::values(s) <- xcov
   names(s) <- 'x'
 
   # Calculate lambda as a function of X and add to Raster Stack
   X <- cbind(1, xcov)
-  values(temp) <- exp(X %*% beta)   # calculate lambda
+  raster::values(temp) <- exp(X %*% beta)   # calculate lambda
   names(temp) <- 'lambda'
-  s <- addLayer(s, temp)  # add to raster stack
+  s <- raster::addLayer(s, temp)  # add to raster stack
 
   # We use rejection sampling to get random draws from the IPP. The proposal
   #  distribution is uniform, ie, HPP:
   # How 'high' should the proposal density be? It must exceed the IPP everywhere.
-  ( maxlambda <- max(values(s)[,'lambda']) ) # ~ 900 per unit area
+  ( maxlambda <- max(raster::values(s)[,'lambda']) ) # ~ 900 per unit area
   (N.hpp <- suppressWarnings(rpois(1, maxlambda*s.area)))     # Number we need to draw
-  if(is.na(N.hpp) || N.hpp >= ncell(s))
+  if(is.na(N.hpp) || N.hpp >= raster::ncell(s))
     stop("The 'beta' settings result in intensities that are too high\nin the most intense region (more animals than pixels.)", call.=FALSE)
   # Draw from the proposal distribution
-  ind.hpp <- sample(1:ncell(s), size = N.hpp, replace = FALSE)   #  sampling w/o replacement ensures only 1 individual per pixel
+  ind.hpp <- sample(1:raster::ncell(s), size = N.hpp, replace = FALSE)   #  sampling w/o replacement ensures only 1 individual per pixel
   # reject draws depending on lambda (and hence X) to get the IPP draws
-  lambda.hpp <- values(s)[,'lambda'][ind.hpp]  # intensities at those pixels
+  lambda.hpp <- raster::values(s)[,'lambda'][ind.hpp]  # intensities at those pixels
   ind.ipp <- rbinom(N.hpp, 1, lambda.hpp/maxlambda)   # use intensity lambda to determine whether to accept or reject.
   (N.ipp <- sum(ind.ipp))               # about 1800 individuals in IPP
   pixel.id.ipp <- ind.hpp[ind.ipp == 1] # Gives id of every pixel in landscape that has an individual, a vector of numbers, length N.ipp = total population
@@ -125,19 +128,19 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   wcov <- getCovSurface(mWt=c(0.25,0.65), sWt=c(0.25, 0.5),rho=0.1, loc=s.loc)
   wcov <- standardize(wcov)            # Standardize covariate w
   # Add to the raster stack
-  values(temp) <- wcov
+  raster::values(temp) <- wcov
   names(temp) <- 'w'
-  s <- addLayer(s, temp)
+  s <- raster::addLayer(s, temp)
 
   # Compute value of thinning parameter b (= probability of detection)
   #  for each cell as a function of alpha and covariate W
   W <- cbind(1, wcov)   # Design matrix W for thinning
-  values(temp) <- plogis(W %*% alpha) # thinning coefs
+  raster::values(temp) <- plogis(W %*% alpha) # thinning coefs
   names(temp) <- 'pTrue'
-  s <- addLayer(s, temp)
+  s <- raster::addLayer(s, temp)
 
   # ... simulate presence-only data (= detections of individuals as a thinned point process)
-  pTrue.ipp <- values(s)[,'pTrue'][pixel.id.ipp]
+  pTrue.ipp <- raster::values(s)[,'pTrue'][pixel.id.ipp]
   y.ipp1 <- rbinom(N.ipp, size=1, prob=pTrue.ipp) # 1 = detected, 0 = not detected
   # A 1/0 vector, length N.ipp
   pixel.id.det1 <- pixel.id.ipp[y.ipp1 == 1]
@@ -153,12 +156,12 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   #### Part C: simulate replicate count data
   ## ---------------------------------------
   # Get the info on animal location into our Raster Stack as a layer:
-  spop <- dropLayer(s, c('lambda','pTrue')) # clean up, just keep covars
-  z <- rep(0, ncell(spop))
+  spop <- raster::dropLayer(s, c('lambda','pTrue')) # clean up, just keep covars
+  z <- rep(0, raster::ncell(spop))
   z[pixel.id.ipp] <- 1
-  values(temp) <- z
+  raster::values(temp) <- z
   names(temp) <- 'presence'  # 1 if animal in pixel, 0 otherwise (max 1 animal per pixel)
-  spop <- addLayer(spop, temp)
+  spop <- raster::addLayer(spop, temp)
 
   # Form quadrats of side quadrat.size (default 4 -> area 16 -> 625 quadrats)
   quadfact <- c(quadrat.size, quadrat.size)
@@ -166,17 +169,17 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   # mean is ok for x and w, but for 'presence' we need the sum
   abund <- raster::aggregate(raster::subset(spop, 'presence'), fact=quadfact, fun=sum)
   names(abund) <- 'N'
-  squad <- addLayer(squad, abund)
-  squad <- dropLayer(squad, 'presence') # clean up, N/16 not useful
+  squad <- raster::addLayer(squad, abund)
+  squad <- raster::dropLayer(squad, 'presence') # clean up, N/16 not useful
 
   # Simulate replicate counts at every quadrat (aka "site")
-  nsite <- ncell(squad)          # number of sites/quadrats in count design
-  N <- values(squad)[,'N']       # Extract latent abundance at each site
+  nsite <- raster::ncell(squad)  # number of sites/quadrats in count design
+  N <- raster::values(squad)[,'N']       # Extract latent abundance at each site
 
   # Compute values of detection probability for each 16-cell pixel in count survey
   # (Here we use the same covar, w, as the detection-only model, as do Koshkina et al,
   #   but a different covar could be generated.)
-  pcount <- plogis(gamma[1] + gamma[2] * values(squad)[,'w'])
+  pcount <- plogis(gamma[1] + gamma[2] * raster::values(squad)[,'w'])
 
   # Do the nsurveys surveys
   counts <- array(NA, dim = c(nsite, nsurveys))
@@ -184,7 +187,7 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
     counts[,j] <- rbinom(nsite, N, pcount)
   }
 
-  fullCountData <- cbind(quadID=1:nsite, values(squad), counts)
+  fullCountData <- cbind(quadID=1:nsite, raster::values(squad), counts)
 
   # Draw random sample of 'nquadrats' quadrats from the total number, nsite
   selQuad <- sort(sample(1:nsite, nquadrats, replace = FALSE))
@@ -219,7 +222,7 @@ simDataDK1 <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
     mnc <- rowMeans(counts)
     mnc[-selQuad] <- NA
     cnt <- raster::subset(squad, 'N')
-    values(cnt) <- mnc
+    raster::values(cnt) <- mnc
     raster::plot(cnt, colNA='darkgrey',axes = FALSE, box = FALSE, asp=1,
       main = "Mean counts for \neach quadrat surveyed,\ngrey if unsurveyed")
   }    # end show.plot
diff --git a/R/simDataDK_AHM2_10.R b/R/simDataDK_AHM2_10.R
index 1184903..6dde22a 100644
--- a/R/simDataDK_AHM2_10.R
+++ b/R/simDataDK_AHM2_10.R
@@ -62,6 +62,9 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   # Note that the parameters beta must be chosen such to avoid a too high 'filling' of the discrete approximation of the entire field B
   # ------------------------------------------------------
 
+  # -------------- Load raster package -----------------------
+  checkNamespace("raster")
+
   # -------------- Check and fix input -----------------------
   sqrt.npix <- round(sqrt.npix[1])
   stopifnotLength(alpha, 2)
@@ -83,9 +86,9 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   s.area <- 4
 
   # Approximate the landscape by many small pixels in a raster object
-  s <- temp <- raster(ncol=sqrt.npix, nrow=sqrt.npix, xmn=-1, xmx=1, ymn=-1, ymx=1, crs=NULL)
+  s <- temp <- raster::raster(ncol=sqrt.npix, nrow=sqrt.npix, xmn=-1, xmx=1, ymn=-1, ymx=1, crs=NULL)
   # s will become a Raster Stack, temp is a template to create new layers
-  s.loc <- xyFromCell(temp, 1:ncell(s))    # Coordinates of every pixel in S
+  s.loc <- raster::xyFromCell(temp, 1:raster::ncell(s))    # Coordinates of every pixel in S
 
   # PART A. Ecological process - where are the animals?
   # ---------------------------------------------------
@@ -95,7 +98,7 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   xcov <- standardize(xcov)    # Standardize covariate x
 
   # Fill covariate x into the raster s
-  values(s) <- xcov
+  raster::values(s) <- xcov
   names(s) <- 'x'
 
   # Calculate log(lambda) as a function of X
@@ -122,9 +125,9 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   # How many in each pixel?
   n <- tabulate(pixel.id.pop, nbins=sqrt.npix^2)
   # Add to the raster stack
-  values(temp) <- n
+  raster::values(temp) <- n
   names(temp) <- 'n'
-  s <- addLayer(s, temp)
+  s <- raster::addLayer(s, temp)
 
   # PART B. The detection-only observation model
   # --------------------------------------------
@@ -132,9 +135,9 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   wcov <- getCovSurface(mWt=c(0.25,0.65), sWt=c(0.25, 0.5),rho=0.1, loc=s.loc)
   wcov <- standardize(wcov)            # Standardize covariate w
   # Add to the raster stack
-  values(temp) <- wcov
+  raster::values(temp) <- wcov
   names(temp) <- 'w'
-  s <- addLayer(s, temp)
+  s <- raster::addLayer(s, temp)
 
   # Compute value of thinning parameter b (= probability of detection)
   #  for each cell as a function of alpha and covariate W
@@ -162,17 +165,17 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
   # mean is ok for x and w, but for 'n' we need the sum
   abund <- raster::aggregate(raster::subset(s, 'n'), fact=quadfact, fun=sum)
   names(abund) <- 'N'
-  squad <- addLayer(squad, abund)
-  squad <- dropLayer(squad, 'n') # clean up, N/16 not useful
+  squad <- raster::addLayer(squad, abund)
+  squad <- raster::dropLayer(squad, 'n') # clean up, N/16 not useful
 
   # Simulate replicate counts at every quadrat (aka "site")
-  nsite <- ncell(squad)          # number of sites/quadrats in count design
-  N <- values(squad)[,'N']       # Extract latent abundance at each site
+  nsite <- raster::ncell(squad)          # number of sites/quadrats in count design
+  N <- raster::values(squad)[,'N']       # Extract latent abundance at each site
 
   # Compute values of detection probability for each 16-cell pixel in count survey
   # (Here we use the same covar, w, as the detection-only model, as do Koshkina et al,
   #   but a different covar could be generated.)
-  pcount <- plogis(gamma[1] + gamma[2] * values(squad)[,'w'])
+  pcount <- plogis(gamma[1] + gamma[2] * raster::values(squad)[,'w'])
 
   # Do the nsurveys counts
   counts <- array(NA, dim = c(nsite, nsurveys))
@@ -180,7 +183,7 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
     counts[,j] <- rbinom(nsite, N, pcount)
   }
 
-  fullCountData <- cbind(quadID=1:nsite, values(squad), counts)
+  fullCountData <- cbind(quadID=1:nsite, raster::values(squad), counts)
 
   # Draw random sample of 'nquadrats' quadrats from the total number, nsite
   selQuad <- sort(sample(1:nsite, nquadrats, replace = FALSE))
@@ -213,7 +216,7 @@ simDataDK <- function(sqrt.npix = 100, alpha = c(-1,-1), beta = c(6,0.5),
       mnc <- rowMeans(counts)
       mnc[-selQuad] <- NA
       cnt <- raster::subset(squad, 'N')
-      values(cnt) <- mnc
+      raster::values(cnt) <- mnc
       raster::plot(cnt, colNA='darkgrey',axes = FALSE, box = FALSE, asp=1,
         main = "Mean counts for \neach quadrat surveyed,\ngrey if unsurveyed")
     }, silent = TRUE)
diff --git a/R/simDynoccSpatial_AHM2_9.R b/R/simDynoccSpatial_AHM2_9.R
index 8d71e2e..bc9a3c1 100644
--- a/R/simDynoccSpatial_AHM2_9.R
+++ b/R/simDynoccSpatial_AHM2_9.R
@@ -102,6 +102,9 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
   # ask.plot – if TRUE permits to browse through plots (otherwise if FALSE)
 
   if(FALSE) {x <- NULL; rm(x)}  # Stops R CMD check choking on 'curve'.
+  
+  # Load raster package------------------------------------------
+  checkNamespace("raster")
 
   # Checks and fixes for input data -----------------------------
   side <- round(side[1])
@@ -228,19 +231,19 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
 
       # Plot random field covariate XAC
       # rows are in x, columns in y direction
-      image(1:side, 1:side, XAC, col=topo.colors(100), 
+      raster::image(1:side, 1:side, XAC, col=topo.colors(100), 
           main = paste("Gaussian random field XAC with \n neg. exponential correlation (range =", theta.XAC, ")"), 
           xlab = 'x', ylab = 'y')
 
-      image(1:side, 1:side, psi[,,1], col=topo.colors(100),
+      raster::image(1:side, 1:side, psi[,,1], col=topo.colors(100),
           main = paste("Initial occupancy probability"), xlab = 'x', ylab = 'y')
 
-      image(1:side, 1:side, z[,,1], col=c("white", "black"), 
+      raster::image(1:side, 1:side, z[,,1], col=c("white", "black"), 
           main = paste("Initial presence/absence (true system state z):\n black = occupied, white = unoccupied"),
           xlab = 'x', ylab = 'y')
       abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
 
-      image(1:side, 1:side, Xauto[,,1], col=topo.colors(100), 
+      raster::image(1:side, 1:side, Xauto[,,1], col=topo.colors(100), 
           main = "Autocovariate between year 1 and year 2", xlab = 'x', ylab = 'y')
     }, silent = TRUE)
     if(inherits(tryPlot, "try-error")) {
@@ -272,15 +275,15 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
     # --------------------
     if(show.plots) {
       tryPlot <- try( {
-        image(1:side, 1:side, phi[,,k-1], col=topo.colors(100),
+        raster::image(1:side, 1:side, phi[,,k-1], col=topo.colors(100),
               main = paste("Persistence between year", k-1, "and year", k), xlab = 'x', ylab = 'y')
-        image(1:side, 1:side, gamma[,,k-1], col=topo.colors(100),
+        raster::image(1:side, 1:side, gamma[,,k-1], col=topo.colors(100),
               main = paste("Colonization between year", k-1, "and year", k), xlab = 'x', ylab = 'y')
-        image(1:side, 1:side, z[,,k], col=c("white", "black"),
+        raster::image(1:side, 1:side, z[,,k], col=c("white", "black"),
             main = paste('Presence/absence (z) in year', k, ':\n black = occupied, white = unoccupied'),
             xlab = 'x', ylab = 'y')
         abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
-        image(1:side, 1:side, Xauto[,,k], col=topo.colors(100),
+        raster::image(1:side, 1:side, Xauto[,,k], col=topo.colors(100),
             main = paste("Autocovariate between year", k, "and year", k+1), xlab = 'x', ylab = 'y')
       }, silent = TRUE)
       if(inherits(tryPlot, "try-error")) {
@@ -312,7 +315,7 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
     for(j in 1:nsurveys){               # Loop over replicates
       prob <- z[,,k] * p[,,j,k]  # zero out p for unoccupied sites
       y[,,j,k] <- rbinom(M, 1, prob)
-    # image(1:side, 1:side, y[,,j,k], main = paste("Year", k, "and rep", j))
+    # raster::image(1:side, 1:side, y[,,j,k], main = paste("Year", k, "and rep", j))
         # Look at clumped pattern in y
     }
   }
@@ -345,9 +348,9 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
       # Plots comparing true and observed latent states
       par(mfrow = c(2,2), mar = c(5,4,5,2), cex.main = 1.3, cex.lab = 1.5, cex.axis = 1.2)
       for(k in 1:nyears){
-        image(1:side, 1:side, z[,,k], col=c("white", "black"), main = paste('Presence/absence (z) in year', k), xlab = 'x', ylab = 'y')
+        raster::image(1:side, 1:side, z[,,k], col=c("white", "black"), main = paste('Presence/absence (z) in year', k), xlab = 'x', ylab = 'y')
         abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
-        image(1:side, 1:side, zobs[,,k], col=c("white", "black"), main = paste('Ever_detected (zobs) in year', k), xlab = 'x', ylab = 'y')
+        raster::image(1:side, 1:side, zobs[,,k], col=c("white", "black"), main = paste('Ever_detected (zobs) in year', k), xlab = 'x', ylab = 'y')
         abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
       }
     }, silent = TRUE)
diff --git a/R/simExpCorrRF.R b/R/simExpCorrRF.R
index 6f20f8c..70bb909 100644
--- a/R/simExpCorrRF.R
+++ b/R/simExpCorrRF.R
@@ -30,6 +30,9 @@ simExpCorrRF <- function(variance = 1, theta = 1, size = 50, seed = NA, show.plo
 # show.plot: if TRUE, plots of the data will be displayed;
 #  set to FALSE if you are running simulations or use inside of other fct's.
 
+# Load raster package
+checkNamespace("raster")
+
 # Generate correlated random variables in a square
 step <- 1
 x <- seq(1, size, step)
@@ -67,7 +70,7 @@ if(show.plots){
     # Random field
     # image(x, y, field,col=topo.colors(20), main = paste("Gaussian random field with \n negative exponential correlation (theta =", theta, ")"), cex.main = 1)
     par(mar = c(3,2,5,1))
-    raster::plot(rasterFromXYZ(cbind(grid, field)), col=topo.colors(20),
+    raster::plot(raster::rasterFromXYZ(cbind(grid, field)), col=topo.colors(20),
     main = paste("Gaussian random field with \n negative exponential correlation (theta =", theta, ")"), cex.main = 1, legend=FALSE, box=FALSE)
     box()
   }, silent = TRUE)
diff --git a/R/simNmixSpatial_AHM2_9.R b/R/simNmixSpatial_AHM2_9.R
index 07de7da..b8b8604 100644
--- a/R/simNmixSpatial_AHM2_9.R
+++ b/R/simNmixSpatial_AHM2_9.R
@@ -29,6 +29,9 @@ simNmixSpatial <- function(nsurveys = 3, mean.lambda = exp(2), beta = c(2, -2),
 BerneseOberland <- NULL # otherwise "no visible binding for global variable 'BerneseOberland'" when checked
 data(BerneseOberland, envir = environment())
 
+# Load raster package
+checkNamespace("raster")
+
 # Simulate spatial random field
 set.seed(seeds[1])
 s <- simExpCorrRF(variance = variance.RF, theta = theta.RF, show.plots=show.plots)
diff --git a/R/simOccSpatial_AHM2_9.R b/R/simOccSpatial_AHM2_9.R
index 730171f..7a892c0 100644
--- a/R/simOccSpatial_AHM2_9.R
+++ b/R/simOccSpatial_AHM2_9.R
@@ -26,6 +26,9 @@ simOccSpatial <- function(nsurveys = 3, mean.psi = 0.6, beta = c(2, -2), mean.p
 BerneseOberland <- NULL # otherwise "no visible binding for global variable 'BerneseOberland'" when checked
 data(BerneseOberland, envir = environment())
 
+# Load raster package
+checkNamespace("raster")
+
 # Simulate spatial random field
 set.seed(seeds[1])
 s <- simExpCorrRF(variance = variance.RF, theta = theta.RF, show.plots = show.plots)
diff --git a/R/simPPe_AHM2_10.R b/R/simPPe_AHM2_10.R
index 9e1d3d2..0187e45 100644
--- a/R/simPPe_AHM2_10.R
+++ b/R/simPPe_AHM2_10.R
@@ -33,6 +33,9 @@ simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X
   #    (this is the parameter of the gridding process
   #    and determines the size of the quadrats)
 
+  # -------------- Load raster package -----------------------
+  checkNamespace("raster")
+
   # -------------- Check and fix input -----------------------
   buffer.width <- round(buffer.width[1])
   stopifNegative(buffer.width)
@@ -125,7 +128,7 @@ simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X
     tryPlot <- try( {
       # *** Fig. 1: Original point pattern
       # Random field of X with activity-centers overlaid
-      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10),
+      raster::image(raster::rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10),
           main = "Point pattern with\ncore and buffer area",
           xlab = "", ylab = "", axes = FALSE, asp = 1)
       mtext(paste("Mean intensity (lambda) =", round(lambda_pp, 5)), side=1)
@@ -138,7 +141,7 @@ simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X
       # *** Fig. 2: Show abundance and presence/absence in each quadrat on original landscape ***
       # Covariate 1: the Gaussian random field with autocorrelation
       # Reproduce random field with activity centers
-      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Abundance, N",
+      raster::image(raster::rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Abundance, N",
           xlab = "", ylab = "", axes = FALSE, asp = 1)
       mtext(paste0("Mean(N) = ", E_N, ", var(N) = ", round(var(c(Nac)), 2)), side=1)
       polygon(c(buffer.width, size.core+buffer.width, size.core+buffer.width, buffer.width),
@@ -161,7 +164,7 @@ simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X
 
       # Figure 3 for presence/absence of activity centers (= distribution)
       # Reproduce random field with activity centers
-      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Occurrence, z",
+      raster::image(raster::rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Occurrence, z",
           xlab = "", ylab = "", axes = FALSE, asp = 1)
       mtext(paste("Mean(z) =", E_z), side=1)
       polygon(c(buffer.width, size.core+buffer.width, size.core+buffer.width, buffer.width),
diff --git a/R/simSpatialDS_AHM1_9-8-3_Simulate_spatial_DS.R b/R/simSpatialDS_AHM1_9-8-3_Simulate_spatial_DS.R
index cde86a2..b1e732a 100644
--- a/R/simSpatialDS_AHM1_9-8-3_Simulate_spatial_DS.R
+++ b/R/simSpatialDS_AHM1_9-8-3_Simulate_spatial_DS.R
@@ -20,6 +20,9 @@ sim.spatialDS <- function(N=1000, beta = 1, sigma=1, keep.all=FALSE,
 #   sigma: scale of half-normal detection function
 #   keep.all: return the data for all individuals, whether detected or not
 
+# Load raster package------------------------------------------
+checkNamespace("raster")
+
 # Checks and fixes for input data -----------------------------
 N <- round(N[1])
 stopifNegative(sigma, allowZero=FALSE)
@@ -71,7 +74,7 @@ y <- rbinom(N, 1, p)                                           # detected or not
 if(show.plot) {
   op <- par(mar=c(3,3,3,6)) ; on.exit(par(op))
   tryPlot <- try( {
-    image(rasterFromXYZ(cbind(as.matrix(gr),x)), col=topo.colors(10), asp=1, bty='n') # need to convert gr to a matrix
+    raster::image(raster::rasterFromXYZ(cbind(as.matrix(gr),x)), col=topo.colors(10), asp=1, bty='n') # need to convert gr to a matrix
     rect(0, 0, 2*B, 2*B)  # draw box around the image
     # draw.circle(B, B, B)
     points(B, B, pch="+", cex=3)
diff --git a/R/simSpatialDSline_AHM2_11-5.R b/R/simSpatialDSline_AHM2_11-5.R
index f9a16d8..18f5f81 100644
--- a/R/simSpatialDSline_AHM2_11-5.R
+++ b/R/simSpatialDSline_AHM2_11-5.R
@@ -9,6 +9,9 @@ function(N=1000, beta = 1, sigma=0.25, alpha0 = -2, W=1/2, L = 4, perp=FALSE, sh
   #   sigma: scale of half-normal detection function
   #   W : truncation distance = strip half-width.
 
+  # Load raster package
+  checkNamespace("raster")
+
   # Create coordinates rasterized transect
   delta<- 0.1                 # '2D bin width'
   # Following code creates coordinates in order of "raster"
@@ -24,7 +27,7 @@ function(N=1000, beta = 1, sigma=0.25, alpha0 = -2, W=1/2, L = 4, perp=FALSE, sh
   # Create spatially correlated covariate x and plot it
   V <- exp(-e2dist(gr,gr)/1)
   x <- t(chol(V))%*%rnorm(nrow(gr))
-  r <- rasterFromXYZ(cbind(gr,x))
+  r <- raster::rasterFromXYZ(cbind(gr,x))
 
   # Simulate point locations as function of habitat covariate x
   probs <- exp(beta*x)/sum(exp(beta*x)) # probability of point in pixel (sum = 1)
@@ -92,7 +95,7 @@ function(N=1000, beta = 1, sigma=0.25, alpha0 = -2, W=1/2, L = 4, perp=FALSE, sh
   if(show.plots) {
     oldpar <- par(mar=c(3,3,3,6), "mfrow") ; on.exit(par(oldpar))
     tryPlot <- try( {
-      image(r, col=topo.colors(10))
+      raster::image(r, col=topo.colors(10))
       abline(0.5, 0, lwd=2)
       image_scale(x, col=topo.colors(10))
       points(u1, u2, pch=20, col='black', cex = 1)  # plot points ### some pixels have >1 animal
diff --git a/R/simSpatialDSte.R b/R/simSpatialDSte.R
index 92ae7c9..14b1f0d 100644
--- a/R/simSpatialDSte.R
+++ b/R/simSpatialDSte.R
@@ -20,6 +20,9 @@ simSpatialDSte <- function(
     theta=2,        # exponential correlation in the spatial covariate
     show.plots=3) {
 
+  # Load raster package------------------------------------------
+  checkNamespace("raster")
+
   # Checks and fixes for input data -----------------------------
   nsites <- round(nsites[1])
   dim <- round(dim[1])
@@ -133,7 +136,7 @@ simSpatialDSte <- function(
         par(mfrow = c(1,2))
       if(nsurveys > 2)
         par(mfrow = c(2,2))
-      img <- rasterFromXYZ(cbind(gr, x[,i]))
+      img <- raster::rasterFromXYZ(cbind(gr, x[,i]))
       for(j in 1:min(nsurveys, 4)) {
         raster::plot(img, col=rampBYR(255), axes=FALSE, box=FALSE)
         title(main=paste("survey", j), line=0.2)
diff --git a/R/simSpatialHDS.R b/R/simSpatialHDS.R
index b9c63cd..095b6fe 100644
--- a/R/simSpatialHDS.R
+++ b/R/simSpatialHDS.R
@@ -13,6 +13,9 @@ function(lam0 = 4 , sigma= 1.5,  B=3, nsites=100, beta1 = 1, npix = 20, show.plo
 # sigma: scale of half-normal detection function
 # B: circle radius
 
+# Load raster package------------------------------------------
+checkNamespace("raster")
+
 # Checks and fixes for input data -----------------------------
 stopifNegative(lam0, allowZero=FALSE)
 stopifNegative(sigma, allowZero=FALSE)
@@ -63,8 +66,8 @@ for(s in 1:nsites){
 
   if(s <= show.plots) {
     tryPlot <- try( {
-      img <- rasterFromXYZ(cbind(gr,z))
-      image(img, col=topo.colors(10))
+      img <- raster::rasterFromXYZ(cbind(gr,z))
+      raster::image(img, col=topo.colors(10))
       #draw.circle(3,3,B)
       image_scale(z,col=topo.colors(10))
       points(u1,u2,pch=16,col='black')