首页 > 解决方案 > 简化 R 中的二进制光栅化

问题描述

我有一些非常小的国家级多边形和点形状文件,我想在 R 中进行栅格化。最终产品应该是一个全局二进制栅格(表明网格单元中心是否被多边形覆盖/点是否位于单元内) . 我的方法是遍历 shapefile 并对每个 shapefile 执行以下操作:

  # load shapefile  
  shp       = sf::read_sf(shapefile_path)
  
  # create a global raster template with resolution 0.0083
  ext       = extent(-180.0042, 180.0042, -65.00417, 75.00417)
  gridsize  = 0.008333333
  r         = raster(ext, res = gridsize)
  
  # rasterize polygon or point shapefile to raster
  rr        = rasterize(shp, r, background = 0) #all grid cells that are not covered get 0

  # convert to binary raster
  values(rr)[values(rr)>0] = 1

在这里,rr是一个栅格文件,其中的多边形/点shp被编码为 1,所有其他网格单元被编码为 0。之后,我将总和rr得到一个包含所有多边形/点的全局二进制栅格文件。

最后两个步骤非常缓慢。此外,当我尝试用 替换所有正值时,我会遇到 RAM 问题,rr因为1由于分辨率高,单元格计数非常大。我想知道是否有可能为我想要实现的目标提出一个更智能的解决方案。

我已经找到了fasterize一个可以快速实现的包,rasterize它运行良好。我认为,如果有人有rasterize直接返回二进制栅格的解决方案,那将有很大帮助。

标签: rgisrastershapefilesf

解决方案


这就是您可以使用raster. 请注意这个value=1论点,以及我更改了您对范围的规范——因为您所做的可能不正确。

library(raster)
v <- shapefile(shapefile_path)
ext <- extent(-180, 180, -65, 75)
r <- raster(ext, res = 1/120)
rr <- rasterize(v, r, value=1, background = 0)

不需要你的最后一步,但你可以做到

rr <- clamp(rr, 0, 1)
# or 
rr <- rr > 0
# or
rr <- reclassify(rr, cbind(1, Inf, 1))

raster::calc对于像这样的简单算术不是很有效

在一个步骤中光栅化所有矢量数据应该快得多,而不是在一个循环中,特别是对于像这样的大型光栅(程序可能需要为每次迭代编写一个临时文件)。

用示例数据说明这个解决方案

library(raster)
cds1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60))
cds2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55))
cds3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45))
v <- spLines(cds1, cds2, cds3)   
r <- raster(ncols=90, nrows=45)

r <- rasterize(v, r, field=1)

为了加快速度,您可以使用terra(光栅的替代品)

library(raster)
f <- system.file("ex/lux.shp", package="terra")
v <- as.lines(vect(f))
r <- rast(v, ncol=75, nrow=100)
x <- rasterize(v, r, field=1)  

推荐阅读