【问题标题】:Rasterize SpatialPolygons with holes in R在 R 中用孔栅格化 SpatialPolygons
【发布时间】:2018-08-10 21:31:01
【问题描述】:

我正在尝试栅格化 SpatialPolygons 对象,以获取与栅格的每个单元格重叠的多边形数量。我的多边形有洞。我遇到的问题是,如果两个孔重叠,则计数中包含此重叠。

我创建了一个玩具示例:

library(sp)
library(raster)

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-180,-20), c(-130,55), c(0,60), c(40,5), c(15,-45), c(-180,-20))
p3 <- list(p3, hole)

pols <- spPolygons(p1, p2, p3)

pols@polygons[[1]]@Polygons[[2]]@hole
pols@polygons[[3]]@Polygons[[2]]@hole


r1 <- raster(ncol=180, nrow=180)
r <- rasterize(pols, r1, fun="count")
plot(pols, col=rgb(1,0,0,0.5))
plot(r)

我认为这可能是导致问题 use R to rasterize shapefile with holes? 相关问题的原因,但没有最终解决方案。问题rasterize ESRI shapefile with holes but FALSE hole slots 的问题是孔的格式不正确,但这不是我的情况。如果您查看这些孔,它们都是 TRUE,请参阅:

pols@polygons[[1]]@Polygons[[2]]@hole
pols@polygons[[3]]@Polygons[[2]]@hole

任何帮助将不胜感激!

【问题讨论】:

    标签: r r-raster sp


    【解决方案1】:

    使用包fasterize 进行光栅化(除了更快)似乎可以避免这个问题:

    library(fasterize)
    library(sf)
    
    r1 <- raster(ncol=180, nrow=180)
    r <- fasterize::fasterize(sf::st_as_sf(pols),
                              r1, fun = "count")
    plot(r, col = rainbow(3))
    

    HTH!

    【讨论】:

    • 太棒了!但是,出于好奇,为什么fasterize::fasterize(sf::st_as_sf(pols), r1, fun = "count") 而不仅仅是fasterize(st_as_sf(pols), r1, fun="count")
    • 只是因为我通常在调用包的函数前加上包名。这是为了避免函数阴影并提醒每个函数属于哪个包。它不会改变任何东西。
    【解决方案2】:

    这是一种解决方法。首先栅格化每个多边形,然后将它们相加。通过这样做,孔将为零。

    library(sp)
    library(raster)
    
    # Create a list of Spatial Polygons
    pols_list <- lapply(list(p1, p2, p3), spPolygons)
    
    # Create a blank raster
    r1 <- raster(ncol = 180, nrow = 180)
    
    # Rasterize each polygon
    r_list <- lapply(pols_list, rasterize, y = r1, fun = "count")
    
    # Create a raster stack
    s <- stack(r_list)
    
    # Calculate the sum
    r2 <- sum(s, na.rm = TRUE)
    
    # Plot r2
    plot(r2)
    

    【讨论】:

      猜你喜欢
      • 2017-05-22
      • 2017-05-02
      • 2017-09-17
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2016-04-13
      相关资源
      最近更新 更多