所以我一直在尝试制作一个未来的 Salsola Tragus 物种的投影模型,但我在接近尾声时遇到了一个障碍
bc <- bioclim(presvals[,c('wc2.1_10m_bioc_CNRM-CM6-1_ssp370_2021-2040','-wc2.1_10m_prec_CNRM-CM6-1_ssp370_2021-2040','wc2.1_10m_tmax_CNRM-CM6-1_ssp370_2021-2040','wc2.1_10m_tmin_CNRM-CM6-1_ssp370_2021-2040')])
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'x' in selecting a method for function 'bioclim': subscript out of bounds
我一直在搞乱配对功能,但我无法弄清楚其他任何事情。
下面是我一直在运行的代码
setwd("~/gfds/Salsola Files")
#loading in the needed library
library(dismo)
library(maptools)
library(raster)
library(sp)
library(rgeos)
library(caret)
library(rJava)
library(gdalUtils)
library(rgdal)
#---misc variables manually defined-----------------------------------------------------------------
e <- extent (-180,180, -80, 80) #Set extent of area (long_min, long_max, lat_min, lat_max)
n_abs <- 1000
datafile <- "Salsola Tragus Data.csv"
data <- read.csv(datafile,sep = ",", quote = "\"",dec = ".", fill = TRUE, comment.char = "")
data <- data[,c("Scientific.Names","Longitude","Latitude")]
salsola.coords <- data[,c("Longitude", "Latitude")]
data.sdf<-data #spacial data frame for salsola data
coordinates(data.sdf) <- ~Longitude+Latitude
projection(data.sdf) <- CRS('+proj=longlat +datum=WGS84')
#---Load worldclim data
#---Read in the Tif Files -------------------------------------------------------------
#Load bioclim stack
future1<-'wc2.1_10m_bioc_CNRM-CM6-1_ssp370_2021-2040.tif'
future2<-'wc2.1_10m_prec_CNRM-CM6-1_ssp370_2021-2040.tif'
future3<-'wc2.1_10m_tmax_CNRM-CM6-1_ssp370_2021-2040.tif'
future4<-'wc2.1_10m_tmin_CNRM-CM6-1_ssp370_2021-2040.tif'
#----------------creating the rasters and stacking them on top of eachother ----------------
ras1<-raster(future1)
ras2<-raster(future2)
ras3<-raster(future3)
ras4<-raster (future4)
Rastack<-stack(ras1,ras2,ras3,ras4) #stacking up all the rasters so that htey can replaced the layered one
worldclim.extent <-stack(crop(Rastack,e)) # adding hte extent to the stack so that hte raster is trimmed to the right size
#---pseudo absence simulation------------------------------------------------------------------------------------------------------------------
#Define a mask from worldclim data
mask_land <- raster(worldclim.extent@layers[[1]])
# Create circles with a radius of 100 km
x <- circles(data.sdf, d=10000, lonlat=TRUE)
presvals <- extract(worldclim.extent, salsola.coords) # start here --------------------------------
# setting random seed to always create the same
# random set of points for this example
set.seed(0)
stvals <- extract(worldclim.extent, randomPoints(worldclim.extent, n_abs)) #simulated absence values
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(stvals)))#presence & absence binary array
sdmdata <- data.frame(cbind(pb, rbind(presvals, stvals)))#dataframe of presence and absence data
#Plot using "pairs" function--------------------------------------------------------------------------------------------------------------------------------
pairs(sdmdata[,2:5], cex=0.1)
#Save outputs to be used in next section
saveRDS(sdmdata, "s_tragus.Rds")
saveRDS(presvals, "s_tragus_pvals.Rds")
#---Model Fitting
#GLM: Linear Models
m2 = glm(pb ~ ., data=sdmdata)
#MaxEnt model
#get maxent() #Run this prior to attempting to load dismo
#Create training data for maxent------------------------------------------------------------------------------------------------------------------------------
group<-kfold(salsola.coords,5)
pres_test <- salsola.coords[group==1, ]
pres_train <-salsola.coords[group !=1,]
#Run MaxEnt: I am not totally sure how "factors" works, ------------------------------------------------------------------------------------------------------
#but it tells MaxEnt which layers we want to consider "categorical". I've set
#this series of model creation steps up to try each of our chosen variables,
#as well as one which doesn't specify a factor (max).
#FYI, MaxEnt has a "removeDuplicates" Boolean option which will remove points
#which fall in the same grid cell. That might be really useful if your abiotic
#rasters have varying resolution.
maxcom <- maxent(worldclim.extent,pres_train)#####THIS SHOULD BE PRES_TRAIN####
maxcom2 <- maxent(worldclim.extent,pres_train,factors='wc2.1_10m_bioc_CNRM-CM6-1_ssp370_2021-2040')
maxcom3 <- maxent(worldclim.extent,pres_train,factors='wc2.1_10m_prec_CNRM-CM6-1_ssp370_2021-2040')
maxcom4 <- maxent(worldclim.extent,pres_train,factors='wc2.1_10m_tmax_CNRM-CM6-1_ssp370_2021-2040')
maxcom5 <- maxent(worldclim.extent,pres_train,factors='wc2.1_10m_tmin_CNRM-CM6-1_ssp370_2021-2040')
#bioclim model-------------------------------------------------------------------------------------------------------------------------------------------
bc <- bioclim(presvals[,c('wc2.1_10m_bioc_CNRM-CM6-1_ssp370_2021-2040','-wc2.1_10m_prec_CNRM-CM6-1_ssp370_2021-2040','wc2.1_10m_tmax_CNRM-CM6-1_ssp370_2021-2040','wc2.1_10m_tmin_CNRM-CM6-1_ssp370_2021-2040')])
任何能指引我正确方向的东西都会有很大帮助