使用 spsample() 对 SpatialPolygonsDataFrame 进行采样会导致 seq.default() 出错

Sampling a SpatialPolygonsDataFrame with spsample() results in error from seq.default()

我正在尝试使用 sp 中的 spsample() 对比利时境内的点进行采样,使用 GDAM 中的 SpatialPolygonsDataFrame。这会导致 seq.default() 调用出错。

be <- readRDS(gzcon(url('http://biogeo.ucdavis.edu/data/gadm2.8/rds/BEL_adm0.rds')))
spsample(be, type="hexagonal", cellsize=10)

Error in seq.default(ll[1], ur[1] - dx/2, dx) : wrong sign in 'by' argument

str(be)

# Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
#   ..@ data       :'data.frame':   1 obs. of  68 variables:
#   .. ..$ OBJECTID     : int 1
# [snip]
#   .. ..$ LDC          : chr ""
#   ..@ polygons   :List of 1
#   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
#   .. .. .. ..@ Polygons :List of 3
#   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
#   .. .. .. .. .. .. ..@ labpt  : num [1:2] 3.37 51.36
# [snip]

进一步尝试:


packageVersion('sp')

# [1] ‘1.2.3’

你的 cellsize 太大了。

查看be的边界框:

bbox(be)
#         min      max
# x  2.555356  6.40787
# y 49.497215 51.50382

我不太明白 spsample 是做什么的,但是如果你查看错误回溯的结果,你可以在要调试的代码:

spsample(be, type="hexagonal", cellsize=10)
traceback()
# 9: stop("wrong sign in 'by' argument")
# 8: seq.default(ll[1], ur[1] - dx/2, dx)
# 7: seq(ll[1], ur[1] - dx/2, dx)
# 6: genHexGrid(dx, bb[, 1], bb[, 2])
# 5: hexGrid(bb, n = n, offset = offset, cellsize = cellsize)
# 4: sample.Spatial(as(x, "Spatial"), n_tot * (1 + its * 0.1), type = type, 
#        offset = offset, ...)
# 3: .local(x, n, type, ...)
# 2: spsample(be, type = "hexagonal", cellsize = 10)
# 1: spsample(be, type = "hexagonal", cellsize = 10)

您可以查看 hexGrid / genHexGrid 的代码(通过键入 getAnywhere("hexGrid") 或者,因为您可以猜测代码在 sp 包中,通过在控制台输入 sp:::hexGrid),你会看到 cellsize 被分配给 dxurbbox(be)max 列; llbbox(be)min 列。

所以您要创建的 x 序列实际上是:

x <- seq(2.555, 1.408, 10)

因此这是试图创建一个从高到低的序列,以正数递增,这是不可能的。

相反,试试更小的 cellsize:

spsample(be, type="hexagonal", cellsize=1)
# SpatialPoints:
#          x        y
# 2 4.205421 50.40138
# 3 5.205421 50.40138
# 4 6.205421 50.40138
# 5 3.705421 51.26741
# 6 4.705421 51.26741
# Coordinate Reference System (CRS) arguments:
# +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
# +towgs84=0,0,0 

我能够通过结合 traceback() 和谨慎使用 debugonce(sp:::sample.Spatial) 来解决这一切。希望这个答案能成为您有用的调试工具箱。