R中如何将Shapefile转换为栅格数据？

Question

我从worldwildlife.org下载了一份有关世界陆地生态区域的shapefile文件。该文件可在以下链接中加载：http://worldwildlife.org/publications/terrestrial-ecoregions-of-the-world。

这是一个标准的shapefile文件，我想对它进行两种操作： 第一步：将本地目录中的shapefile裁剪到东北美洲的范围（ext = extent(-95,-50,24,63)）。

# Read shapefile using package "maptools"
eco_shp <- readShapeLines("F:/01_2013/Ecoregions/Global/wwf_terr_ecos.shp", 
                          proj4string=CRS("+proj=utm +zone=33 +datum=WGS84")) 


# Set the desired extent for the final raster using package "raster" 
ext <- extent(-95, -50, 24, 63)

我确定我需要使用"raster"包中的栅格化函数，但我仍然无法正确地使用它。希望能得到关于如何使用它的任何建议。

- I Del Toro

---

Mnel：我确实需要将它光栅化为.asc文件，以与我的其他环境图层匹配。 - I Del Toro

为什么不直接将多边形叠加到您希望匹配的SpatialGridDataFrame上（使用overlay）？ - mnel

我的模型要求输入的数据框必须都是asc格式，并且具有相同的空间范围。如果我只想可视化它，那么"overlay"就可以工作了。 :) - I Del Toro

R中的overlay函数将提取网格点处的值，然后可以保存为ASCII文件。 - mnel

很抱歉，我在使用R进行GIS工作方面还比较新。您的意思是创建一个具有所需范围的空栅格，然后叠加形状文件并以此方式提取点的值吗？ - I Del Toro

Answer 1

你的想法是正确的，对于空间栅格数据，你应该使用raster（而不是sp空间栅格类）。在读取、写入和其他操作空间矢量数据时，你也应该使用rgdal（而不是maptools）。

这将帮助你入门：

library(rgdal)
library(raster)

## Read in the ecoregion shapefile (located in R's current working directory)
teow <- readOGR(dsn = "official_teow/official", layer = "wwf_terr_ecos")

## Set up a raster "template" to use in rasterize()
ext <-  extent (-95, -50, 24, 63)
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
n <- 5
r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)

## Rasterize the shapefile
rr <-rasterize(teow, r)

## A couple of outputs
writeRaster(rr, "teow.asc")
plot(rr)

- Josh O'Brien

---

很高兴听到这个消息。请注意，在实践中，您可能希望使用自己的栅格图层作为模板栅格（rasterize(teow, r) 中的 r），并且可能需要进行一些调整以使 proj4strings 匹配（尽管 raster 和 rgdal 都非常擅长处理投影元数据）。 - Josh O'Brien

需要注意的是，对于大型数据集，rasterize() 函数似乎非常不稳定和低效。最近我尝试以6英尺分辨率栅格化了西雅图市的覆盖范围（输出为14 MB 的 GeoTIFF 文件，大小为8095 x 14819），R 运行了7个线程约三个小时才“完成”，但没有任何输出和错误信息。使用 R 生成所需范围和分辨率的空白 GeoTIFF 栅格，然后通过 gdal（在 QGIS 的帮助下）运行 rasterize 操作，在单个线程中不到半小时即可完成。 - Miles Erickson

确实。对于除最小光栅外的任何光栅，我倾向于使用 gdalUtils::gdal_rasterize 的包装器。对于 1000x1000 的光栅，我可以获得 6 倍的速度提升，对于 2000x2000 的光栅，我可以获得 12 倍的速度提升（从 70 秒到 6 秒）。 - Josh O'Brien

@BhargavRao 我确实做了。谢谢！（已修复。） - Josh O'Brien

对于一个79862x67313的栅格，gdal（python）需要大约100秒。R raster包需要一整天时间并且通常会崩溃。区别在于gdal是C应用程序并且经过了优化...而raster->rasterize的代码是用R编写的，并且使用了太多的内存。 - user1269942