如何栅格化 SpatialPolygonsDataFrame 以便我可以用点（在 R 中）填充 shapefile？答案

【问题标题】：How can I rasterize a SpatialPolygonsDataFrame so that I can fill a shapefile with points (in R)?如何栅格化 SpatialPolygonsDataFrame 以便我可以用点（在 R 中）填充 shapefile？
【发布时间】：2022-08-18 18:29:22
【问题描述】：

我有一大块 R 代码曾经可以工作，但现在不能工作了，我找不到问题所在。代码的目的是用规则间隔的点填充 shapefile。

我的 shapefile 可以在这里访问：https://drive.google.com/drive/folders/1SAbuyIQHevK4fz-0w3TTqpEhz0wKLEII?usp=sharing

如果我开始加载我的 shapefile：

GUA = raster::shapefile(\'Guam3BufferPoly.shp\')

然后我为这个 SpatialPolygonsDataFrame 的坐标参考系统设置了一个变量：

projGUA = crs(GUA)

转换为平面 crs

putm <- spTransform(GUA, projGUA)

创建一个栅格（这是它不起作用的地方）

ext = extent(putm)
r <- raster(ext, res=500)

栅格化多边形并转换为点

r2 <- rasterize(putm, r)
pts <- rasterToPoints(r2, spatial=TRUE)

将点转换为 lon/lat 并绘制结果

pts_lonlat <- spTransform(pts, \"+proj=longlat +datum=WGS84\")
plot(pts_lonlat,pch=\'*\')

栅格 r 为空（破坏了下游的所有代码）。

如果你能帮助我，请告诉我。请善待（这是我第一次在这里发帖，如果我没有正确格式化我的问题，我深表歉意）。谢谢！

标签： r raster shapefile rasterize

【解决方案1】：

欢迎！

基本上，您似乎在一个不适合您执行的步骤的坐标参考系统中工作。

library(raster)

GUA <- shapefile("Guam3BufferPoly.shp")

crs(GUA)
#> Coordinate Reference System:
#> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
#> WKT2 2019 representation:
#> GEOGCRS["GCS_unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["Degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["Degree",0.0174532925199433]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["Degree",0.0174532925199433]]]

似乎 Guam3BufferPoly.prj 提供的 crs 不太符合标准，但我在这里不是 100% 确定。但是，您正在使用 WGS 84 上的纬度/经度坐标。因此，spTransform(GUA, projGUA) 在这里绝对没有效果，因为您再次将您的功能重新投影到相同的 crs。

这导致了一个后续问题：您以度为单位工作并指定 500 个单位的栅格（此处：度，而不是米，可能是预期的）。这正是raster() 所做的，即使它没有任何意义：

ext <- extent(GUA)
ext
#> class      : Extent 
#> xmin       : 144.5962 
#> xmax       : 145 
#> ymin       : 13.19925 
#> ymax       : 13.70405

r <- raster(ext, res = 500)
r
#> class      : RasterLayer 
#> dimensions : 1, 1, 1  (nrow, ncol, ncell)
#> resolution : 500, 500  (x, y)
#> extent     : 144.5962, 644.5962, -486.2959, 13.70405  (xmin, xmax, ymin, ymax)
#> crs        : NA

要解决此问题，您应该使用适合关岛的投影坐标参考系。不确定这有多合适，您可能会知道得更好，但让我试试 WGS 84 / UTM 55 S (EPSG: 32755)：

putm <- spTransform(GUA, CRS("+proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs"))

ext <- extent(putm)
ext
#> class      : Extent 
#> xmin       : 239608.8 
#> xmax       : 283618.4 
#> ymin       : 11460271 
#> ymax       : 11516045

r <- raster(ext, res = 500)
r
#> class      : RasterLayer 
#> dimensions : 112, 88, 9856  (nrow, ncol, ncell)
#> resolution : 500, 500  (x, y)
#> extent     : 239608.8, 283608.8, 11460045, 11516045  (xmin, xmax, ymin, ymax)
#> crs        : NA

r2 <- rasterize(putm, r)
pts <- rasterToPoints(r2, spatial = TRUE)

pts_lonlat <- spTransform(pts, "+proj=longlat +datum=WGS84")
plot(pts_lonlat, pch='*')

【讨论】：