3

我在栅格化 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)

有什么想法我可能会出错吗？形状文件本身有问题吗？

请和谢谢！

4

1 回答 1

7

您假设多边形在 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。

于 2016-01-30T05:40:16.057 回答