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 2x2 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.