【发布时间】：2020-09-21 22:27:11
【问题描述】：
我是 R 的新手，并试图了解它是如何为某些物种分布建模完成的。我正在尝试使用 maptools 包中包含的 wrld_simpl 地图创建意大利的面具。当我运行代码时，我在尝试创建 spxy 时收到“.local(obj, ...) 中的错误：坐标中的 NA 值”。我怀疑这个问题可能源于光栅化步骤，但我不确定......我做错了什么？
为代码转储道歉
# Create library
library(raster)
library(dismo)
library(sf)
library(maptools)
library(rgdal)
library(sp)
library(rgeos)
# Call the cleaned nivale csv
nivale <- read.table('C:/Users/David/Documents/nival_1980_GeoDat.csv', header=T, sep = ',')
nivale <- nivale[,2:3]
# Create a simple map of Italy, plot nivale data points
data(wrld_simpl)
plot(wrld_simpl, xlim=c(0,20), ylim=c(40,50), axes=TRUE, col="light yellow")
box()
points(nivale$lon, nivale$lat, col='orange', pch=20, cex=0.75)
points(nivale$lon, nivale$lat, col='black', cex=0.75)
# set CRS of nivale equal to wrld_simpl
coordinates(nivale) <- ~lon+lat
crs(nivale) <- crs(wrld_simpl)
projection(nivale) <- CRS('+proj=longlat +datum=WGS84')
class(nivale)
class(wrld_simpl)
# Sampling Bias Assessment
r <- raster(nivale)
res(r) <- 1
r <- extend(r, extent(r)+1)
nisel <- gridSample(nivale, r, n=1)
p <- rasterToPolygons(r)
plot(p, border='blue')
points(nivale)
points(nisel, cex=1, col='red', pch='x')
# Pseudo-Absences
# Create shapefile from wrld_simpl
italy <- wrld_simpl[is.element(wrld_simpl$NAME, 'Italy'),]
set.seed(1963)
crsi = crs('+proj=longlat +datum=WGS84')
exti = extent(italy)
# Create template for rasterization
rst_temp <- raster(ncols = 1000, nrows = 1000,
crs = crsi,
ext = exti)
# Rasterize italy
rst_italy <- rasterize(italy, rst_temp)
# Random point generation
rand_point <- randomPoints(rst_italy, 250)
#Pseudo-Absence Points
x <- circles(nivale, d=50000, lonlat=TRUE)
pol <- polygons(x)
samp1 <- spsample(pol, 250, type='random', iter=25)
cells <- cellFromXY(rst_italy, samp1)
xy <- xyFromCell(rst_italy, cells)
plot(pol, axes=TRUE)
points(xy, cex=0.75, pch=20, col='blue')
#
# Error - NA values in coordinates
spxy <- SpatialPoints(xy, proj4string=CRS('+proj=longlat +datum=WGS84'))
o <- over(spxy, geometry(x))
xyInside <- xy[!is.na(o), ]
v <- extract(mask, x@polygons, cellnumbers=T)
v <- do.call(rbind, v)
v <- unique(v[,1])
head(v)
m <- italy
m[] <- NA
m[v] <- 1
plot(m, ext=extent(x@polygons)+1)
plot(x@polygons, add=T)
【问题讨论】：