I have been trying to conduct a MaxEnt test in R and have run into an issue that I can’t solve. I would appreciate it if someone could help me. My environmental raster data and presence data are loading correctly. My pairwise correlation is showing accurate results. The issue is probably something straightforward, but the answer is eluding me. Here is my complete code and the error message I am receiving:

# Introduction ————————————————————

 

# Pteoris volitans occurrences from GBIF (coordinates) along with Bio-Oracle SST 5.85

# Worst Case Scenario SSP5-8.5.

 

# Load Packages ———————————————————–

 

# Load necessary packages

library(installr)

library(rgbif)

library(rJava)

library(occ)

library(sp)

library(raster)

library(ggplot2)

library(lattice)

library(rasterVis)

library(dismo)

library(GGally)

library(tidyverse)

library(sf)

library(terra)

library(SDMtune)

library(biomod2)

library(zeallot)

library(gridExtra)

library(farver)

 

# Load GBIF Data on P. volitans ——————————————-

 

# Perform GBIF occurrence search

gbif_p.volitans_data <- rgbif::occ_data(scientificName = “Pterois volitans”, hasCoordinate = TRUE, limit = 13000)

myspecies_coords <- gbif_p.volitans_data$data[, c(“decimalLongitude”, “decimalLatitude”, “individualCount”, “occurrenceStatus”, “coordinateUncertaintyInMeters”, “institutionCode”, “references”)]

 

# Convert the list to a data frame

myspecies_coords <- data.frame(decimalLatitude = myspecies_coords$decimalLatitude, decimalLongitude = myspecies_coords$decimalLongitude)

 

# Change name of myspecies_coords

p_volitans_coords <- myspecies_coords

 

# Transform GBIF Data to .csv ———————————————-

 

# Write data frame to CSV

write.csv(myspecies_coords, file = “/Users/micahgisclair/Documents/UCC Marine Biology MSc/Dissertation/R_code/GBIF_and_BioOracle/GBIF_and_Bio_Oracle_R_SSP585/Pterois_volitans_Occurence/myspecies_coords.csv”, row.names = FALSE)

 

# Extract Coordinates from GBIF Data ————————————–

 

# Extract latitude and longitude from GBIF Pterois volitans data

decimallatitude <- gbif_p.volitans_data$data$decimalLatitude

decimallongitude <- gbif_p.volitans_data$data$decimalLongitude

 

# Path to Bio-Oracle SST SSP5-8.5 Layer ————————————–

 

nc_path1 <- “/Users/micahgisclair/Documents/UCC Marine Biology MSc/Dissertation/R_code/GBIF_and_BioOracle/GBIF_and_Bio_Oracle_R_SSP585/data/sst/SST_ssp585_2020_2100_depthsurf_773f_eb2c_59ac_U1713367461606.nc”

 

# Read Bio-Oracle layer SST SSP5-5.85—————————————-

 

bo_SST_SSP585 <- raster(nc_path1)

 

# Saving Bio-Oracle SST SSP5-5.85 Raster Tiff ——————————————-

 

# Write Bio-Oracle SST SSP5-5.85 layer as GeoTIFF and .csv

# Save the Data as a Raster GeoTIFF

writeRaster(bo_SST_SSP585, “./data/sst/bo_SST_SSP585_terra.tif”, overwrite = TRUE)

 

# Plot SST raster————————————————————

 

plot(bo_SST_SSP585, main = “Sea Surface Temperature”)

 

# Path to Bio-Oracle SWS SSP5-8.5 Layer ————————————-

 

nc_path2 <- “/Users/micahgisclair/Documents/UCC Marine Biology MSc/Dissertation/R_code/GBIF_and_BioOracle/GBIF_and_Bio_Oracle_R_SSP585/data/sws/sws_ssp585_2020_2100_depthsurf_b657_a572_c6b1_U1715786576023.nc”

 

# Read Bio-Oracle SWS SSP5-8.5 Layer—————————————–

 

bo_SWS_SSP585 <- raster(nc_path2)

 

# Saving Bio-Oracle SWS SSP5-8.5 Raster Tiff ——————————

 

writeRaster(bo_SWS_SSP585, “./data/sws/bo_SWS_SSP585_terra.tif”, overwrite=TRUE)

 

# Plot SWS raster————————————————————

 

plot(bo_SWS_SSP585, main = “Sea Water Speed”)

 

# Path to Bio-Oracle SO SSP5-8.5 Layer ————————————–

 

nc_path3 <- “/Users/micahgisclair/Documents/UCC Marine Biology MSc/Dissertation/R_code/GBIF_and_BioOracle/GBIF_and_Bio_Oracle_R_SSP585/data/so/so_ssp585_2020_2100_depthsurf_9ed4_e29f_e3fe_U1715785506682.nc”

 

# Read Bio-Oracle SO SSP5-8.5 Layer——————————————

 

bo_SO_SSP585 <- raster(nc_path3)

 

# Saving Bio-Oracle SO SSP5-8.5 Raster Tiff ——————————

 

# Write Bio-Oracle SO SSP5-5.85 layer as GeoTIFF

writeRaster(bo_SO_SSP585, “./data/so/bo_SO_SSP585_terra.tif”, overwrite =TRUE)

 

# Plot SO raster———————————————————-

 

plot(bo_SO_SSP585, main = “Salinity”)

 

# Path to Bio-Oracle Phyc SSP5-8.5 Layer ———————————-

 

nc_path4 <- “/Users/micahgisclair/Documents/UCC Marine Biology MSc/Dissertation/R_code/GBIF_and_BioOracle/GBIF_and_Bio_Oracle_R_SSP585/data/phyc/Primary Productivity_ssp585_2020_2100_depthsurf_7f72_0431_8307_U1713367491076.nc”

 

# Read Bio-Oracle Phyc SSP5-8.5 Layer————————————–

 

bo_phyc_SSP585 <- raster(nc_path4)

 

# Saving Bio-Oracle Phyc SSP5-8.5 Raster Tiff and .csv ——————————

 

# Write Bio-Oracle Phyc SSP5-5.85 layer as GeoTIFF

writeRaster(bo_phyc_SSP585, “./data/phyc/bo_phyc_SSP585_terra.tif”, overwrite =TRUE)

 

# Plot Phyc raster——————————————————————-

 

plot(bo_phyc_SSP585, main = “Primary Productivity”)

 

# Load layer data to convert into polygons within a grid —————–

 

# Read geotiff files

# Load GeoTIFF files

bo_SST_SSP585 <- raster(“./data/sst/bo_SST_SSP585_terra.tif”)

bo_SO_SSP585 <- raster(“./data/so/bo_SO_SSP585_terra.tif”)

bo_SWS_SSP585<- raster(“./data/sws/bo_SWS_SSP585_terra.tif”)

bo_phyc_SSP585 <- raster(“./data/phyc/bo_phyc_SSP585_terra.tif”)

 

# Load grid data

grid <- raster::rasterToPolygons(raster::raster(xmn = -180, xmx = 180, ymn = -90, ymx = 90, resolution = 1))

 

# Plot grid

plot(grid)

 

# Use sf (simple features) to convert grid and give unique ID to polygons———————————————

 

# Create sf grid from current grid information

sf_grid <- st_as_sf(grid, wkt = “geometry”)

sf_grid$ID <- 1:nrow(grid)

points_sf <- st_as_sf(sf_grid, coords = c(“decimalLongitude”,”decimalLatitude”), crs = 4326)

sf_grid$ID <- seq_len(nrow(sf_grid))

join_data <- st_join(points_sf, sf_grid, join = st_within)

 

# Reduce down to 1 point per cell

oppc <- join_data %>% group_by(ID.y) %>% sample_n(1) %>% ungroup()

 

# Convert to Spatial object

sp_oppc <- as(oppc, “Spatial”)

sp_oppc <- as.data.frame(sp_oppc)

 

# Comprising the final_df for pairwise correlation test———————————————————–

 

# Plot each environmental raster layer

par(mfrow = c(2, 2))  # Set up a 2×2 grid for plots

plot(bo_SST_SSP585, main = “sst”)

plot(bo_SO_SSP585, main = “so”)

plot(bo_SWS_SSP585, main = “sws”)

plot(bo_phyc_SSP585, main = “phyc”)

par(mfrow = c(1, 1))  # Reset the plotting layout to default

 

# Create a dataframe with species occurrences

# Calculate the number of repetitions needed for each layer

repetitions <- ceiling(12996 / 4)

 

# Create a dataframe with species occurrences

final_df <- data.frame(

  species = rep(“Pterois volitans”, repetitions * 4),  # Repeat the species name repetitions * 4 times

  PA = rep(1, repetitions * 4),  # Repeat the value 1 for PA repetitions * 4 times

  layer = rep(c(“bo_SST_SSP585”, “bo_SO_SSP585”, “bo_SWS_SSP585”, “bo_phyc_SSP585”), each = repetitions),  # Repeat each layer name repetitions times

  ID = rep(1:repetitions, 4),  # Sequential IDs repeated 4 times

  decimalLongitude = runif(repetitions * 4, min = -180, max = 180),  # Random longitudes

  decimalLatitude = runif(repetitions * 4, min = -90, max = 90)  # Random latitudes

)

 

# Set coordinates

coordinates(final_df) <- c(“decimalLongitude”, “decimalLatitude”)

 

# Create a list to store extracted values for each variable

extracted_values <- list()

 

# Extract values at occurrences for each variable and store them in the list

extracted_values$sst <- raster::extract(bo_SST_SSP585, final_df)

extracted_values$so <- raster::extract(bo_SO_SSP585, final_df)

extracted_values$sws <- raster::extract(bo_SWS_SSP585, final_df)

extracted_values$phyc <- raster::extract(bo_phyc_SSP585, final_df)

 

# Add extracted values to final_df

final_df <- cbind(final_df, extracted_values)

 

# Rename the columns of the SpatialPointsDataFrame

names(final_df@data) <- c(“species”, “PA”, “layer”, “ID”, “sst”, “so”, “sws”, “phyc”)

 

# Convert final_df to a data frame

final_df <- as.data.frame(final_df)

 

# Subset final_df for correlation analysis

correlationdata <- final_df[, c(“sst”, “so”, “sws”, “phyc”)]

 

# Create predictors stack

predictors <- stack(bo_SST_SSP585, bo_SO_SSP585, bo_SWS_SSP585, bo_phyc_SSP585)

names(predictors) <- c(“bo_SST_SSP585”, “bo_SO_SSP585”, “bo_SWS_SSP585”, “bo_phyc_SSP585”)

plot(predictors)

 

# Filter presence and absence

presence <- subset(final_df, PA == 1)

absence <- subset(final_df, PA == 0)

 

# Convert presence and absence data frames to a two-column matrix of coordinates

presence_coords <- cbind(presence$decimalLongitude, presence$decimalLatitude)

absence_coords <- cbind(absence$decimalLongitude, absence$decimalLatitude)

 

# Convert predictors to SpatRaster object

predictors_spat <- terra::rast(predictors)

print(predictors_spat)

 

# Utilize ggllay (ggpairs) to make a pairwise correlations to test (30%) presence/absence data—————————————–

 

# Randomly sample points (same number as our observed points)

background <- terra::spatSample(x = predictors_spat, size = 12996,  # generate pseudo-absence points equal to the number of observed points

                                values = FALSE,  # don’t need values

                                na.rm = TRUE,   # don’t sample from NA values

                                xy = TRUE)      # just need coordinates

 

# Convert background to SpatRaster object

background_terra <- terra::vect(background)

 

# Use prepareSWD from the ‘dismo’ package with ‘terra’ objects

predictors_terra <- prepareSWD(

  species = “Pterois volitans”,

  p = presence_coords,

  a = background,

  env = predictors_spat

)

 

# Split data into training, validation, and testing sets

datasets <- trainValTest(predictors_terra, test = 0.3, only_presence = FALSE, seed = 25)

train <- datasets[[1]]

test <- datasets[[2]]

 

# Subset final_df for correlation analysis

correlationdata <- final_df[, c(“sst”, “so”, “sws”, “phyc”)]

 

# Remove NAs from the correlation data frame

correlationdata_df <- na.omit(as.data.frame(correlationdata))

 

# Plot pairwise correlations

GGally::ggpairs(correlationdata_df) + theme_bw()

 

# Perform MaxEnt test —————————————————–

 

# Check the number of observations in the data slot (coordinates of environmental layers)

# Assuming train is your SWD object

coords_df <- train@coords  # Extracting the coordinates data frame

 

# Creating the data_obs object with X and Y columns

data_obs <- coords_df[, c(“X”, “Y”)]

 

# Displaying the first few rows of data_obs to verify

head(data_obs)

 

# Check the length of the pa slot (presence/absence data)

pa_obs <- length(train@pa)

print(pa_obs)

 

# Perform MaxEnt test

model_maxent_SSP585 <- train(data = train, method = “Maxent”, fc = “lqph”, reg = 1)

 

# Perform MaxEnt test

model_maxent_SSP585 <- train(data = train, method = “Maxent”, fc = “lqph”, reg = 1)

Error in data.frame(…, check.names = FALSE) : arguments imply differing number of rows: 1, 2, 0

I believe it may be my data_obs and pa_obs results. They both show just one column, one row, and the total presence data.