将 shapefile 转换为光栅

Question

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

我从另一个工作正常的示例中提取了代码，但是当我为我的数据尝试它时，我从 plot 函数中什么也得不到。有关 shapefile 和我的代码的链接，请参见下文：

珊瑚礁面临风险：全球综合威胁

# 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)

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

但是以上所有内容都没有真正相关，因为您正在逆向执行此操作。多边形显然来自栅格（看看它们是多么块状）并且栅格在您指向的数据集中可用！

因此，不要进行光栅化，而是执行以下操作：

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。