I have created a lot of maps before without an issue and I am wondering if this is because of a simple update requirement or something similar. When I use ggplot + geom_stars for creating occurrence maps with spOccupancy prediction, I get a warning message and the map drawn does not capture details of the map. "Warning message: Removed 38051 rows containing missing values or values outside the scale range (geom_raster()). "

Prediction is for 30910 non-sampled locations and I am confused where the 38051 is coming from. Projection used is NAD83, when I change it to UTM, the map doesn't even load and the error says "Error: cannot allocate vector size of xxx " which I know is due to memory but I do have more memory in the computer than it says it needs. Therefore, I am stuck. All other codes work and the analysis is complete, I just can't show it in the map. I have restarted a couple of times and cleared memory but nothing works. Any suggestions?

Below is the code I used and files are attached:

library(spOccupancy)
library(MCMCvis)
library(ggplot2)
library(sf)
library(stars)

setwd('')
y<-read.csv("SG_Spring_finaldata.csv", header=TRUE)

Plot (site) codes.

plot.codes <- sort(unique(y$FID))
# Maximum number of replicates at a site
K <- 2
# Number of sites
J <- length(unique(y$FID))
# Array for detection-non detection data. 
ymat <- array(y$Join_Count, dim = c(J, K),)

# Look at the structure of our array y
str(ymat)

Site covariate and coordinates

sc <- read.csv("occ_covs.csv", header=TRUE)
coords<-read.csv("coords.csv", header=TRUE)
coords1<-as.data.frame(coords)
coords<-as.matrix(coords1)
str(coords)
str(sc)

Format detection covs

library(dplyr)
PH <- y %>%
  group_by(YEAR,FID) %>%
  summarize(PH = unique(Perren_herb) 
  ) %>%
  ungroup() %>%
  glimpse()

SB <- y %>%
  group_by(YEAR,FID) %>%
  summarize(SB = unique(Sagebrush) 
  ) %>%
  ungroup() %>%
  glimpse()

PH<- matrix(PH$PH, nrow = J, ncol = K)
SB<-matrix(SB$SB, nrow = J, ncol = K)

det.covs <- list(PH=PH,
                 SB=SB)

str(det.covs)

data <- list(y=ymat,
             occ.covs=sc,
             det.covs=det.covs,
             coords=coords)
str(data)

Model

output <- spPGOcc(occ.formula = ~ scale(Rugged)+ scale(Elevation),
                  det.formula = ~ scale(PH)+ scale(SB),
                  data = data,
                  n.batch = 400,
                  batch.length = 25,
                  NNGP = TRUE, 
                  n.neighbors = 5,
                  n.thin = 10, 
                  n.burn = 3000, 
                  n.chains = 5,
                  n.report = 100)

ppc.output <- ppcOcc(output, fit.stat = 'freeman-tukey', group = 1)
summary(ppc.output)
waicOcc(output)

predicting throughout the landscape

data2<-read.csv("base_data2.csv", header=TRUE)
coords2<-read.csv("base_coords500.csv", header=TRUE)
coords3<-as.data.frame(coords2)
coords.0<-as.matrix(coords3)
str(coords.0)

elevation.0 <- (data2[, 'Elevation'] - mean(data$occ.covs$Elevation)) / 
  sd(data$occ.covs$Elevation)
rug.0 <- (data2[, 'Rugged'] - mean(data$occ.covs$Rugged)) / 
  sd(data$occ.covs$Rugged)

X.0 <- cbind(1, rug.0, elevation.0)
out.pred <- predict(output, X.0, coords.0)
str(out.pred)

for occurence

psi.0.mean <- apply(out.pred$psi.0.samples, 2, mean)
# Occupancy probability standard deviations
psi.0.sd <- apply(out.pred$psi.0.samples, 2, sd)
# Spatial process mean and sd
w.0.mean <- apply(out.pred$w.0.samples, 2, mean)
w.0.sd <- apply(out.pred$w.0.samples, 2, sd)

plot.df <- data.frame(psi.mean = psi.0.mean,
                      psi.sd = psi.0.sd,
                      w.mean = w.0.mean, 
                      w.sd = w.0.sd,
                      x = coords.0[, 1],
                      y = coords.0[, 2])

str(plot.df)
pred.stars <- st_as_stars(plot.df, dims = c('x', 'y'))

plotting occurence

library(viridis)

ggplot() +
  geom_stars(data = pred.stars, aes(x = x, y = y, fill = psi.mean),interpolate = TRUE) +
  scale_fill_viridis("", limits = c(0, 1),
                     na.value = NA) +
  theme_bw(base_size = 18) + 
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank()) +
  labs(x = "Easting", y = "Northing", title = 'Spring')

base_coords500.csv base_data2.csv coords.csv occ_covs.csv SG_Spring_finaldata.csv

biodiverse / spOccupancy

Error while creating maps regarding removing rows, not producing clear maps #48