将 shapefile 转换为栅格

Question

我在栅格化 shapefile 以在 0.5*0.5 网格上生成点时遇到问题。 shapefile 表示全球珊瑚礁对综合威胁的风险级别分类（低-0、中-100、高-1000、极高-1500）。

我从另一个工作正常的示例中提取了代码，但是当我为我的数据尝试它时，我从 plot 函数中没有得到任何结果。请参阅下面的 link 到 shapefile 和我的代码：

Reefs At Risk: Global Integreated Threats

# Read shapefile into R
library(rgdal)
library(raster)    

int.threat.2030 <- readOGR(dsn = "Global_Threats/Integrated_Future", 
                           layer = "rf_int_2030_poly")

## Set up a raster "template" for a 0.5 degree grid
ext <- extent(-110, -50, 0, 35)
gridsize <- 0.5
r <- raster(ext, res=gridsize)

## Rasterize the shapefile
rr <- rasterize(int.threat.2030, r)

## Plot raster
plot(rr)

我有什么地方可能出错的想法吗？是 shapefile 本身的问题吗？

拜托谢谢！

Answer 1

您假设多边形在 lon/lat 坐标中，但它们不是：

library(raster)
library(rgdal)
p <- shapefile('Global_Threats/Integrated_Future/rf_int_2030_poly.shp')
p

#class       : SpatialPolygonsDataFrame 
#features    : 63628 
#extent      : -18663508, 14601492, -3365385, 3410115  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=cea +lon_0=-160 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 3
#names       :    ID, THREAT, THREAT_TXT 
#min values  :     1,      0,   Critical 
#max values  : 63628,   2000,  Very High

您可以更改投影

pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

然后做类似的事情：

ext <- floor(extent(pgeo))
rr <- raster(ext, res=0.5)
rr <- rasterize(pgeo, rr, field=1)

或者保留原始 CRS 并执行以下操作：

ext <- extent(p)
r <- raster(ext, res=50000)  
r <- rasterize(p, r, field=1)
plot(r)

请注意，您正在将非常小的多边形栅格化为大栅格单元。如果多边形覆盖单元格的中心（即假设多边形覆盖多个单元格的情况），则多边形被视为 'inside'。因此，对于这些数据，您需要使用更高的分辨率（然后可能汇总结果）。或者，您可以栅格化多边形质心。

但是上面的 none 确实是相关的，因为你正在向后做这一切。多边形显然是从栅格派生的（看看它们有多块），并且栅格在您指向的数据集中可用！

所以不要光栅化，而是：

x <- raster('Global_Threats/Integrated_Future/rf_int_2030')
x
#class       : RasterLayer 
#dimensions  : 25456, 80150, 2040298400  (nrow, ncol, ncell)
#resolution  : 500, 500  (x, y)
#extent      : -20037508, 20037492, -6363885, 6364115  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 
#data source : C:\temp\Global_Threats\Integrated_Future\rf_int_2030 
#names       : rf_int_2030 
#values      : 0, 2000  (min, max)
#attributes  :
#   ID  COUNT THREAT_TXT
#    0  80971        Low
#  100 343535     Medium
# 1000 322231       High
# 1500 168518  Very High
# 2000  83598   Critical

这里绘制了巴拉望岛的一部分：

e <- extent(c(-8990636, -8929268, 1182946, 1256938))
plot(x, ext=e)
plot(p, add=TRUE)

如果您需要较低的分辨率，请参阅 raster::aggregate。对于不同的坐标参考系，请参见 raster::projectRaster。

将 shapefile 转换为栅格

Converting shapefile to raster

r

shapefile

rasterize