8

我有一个简单的栅格（使用 R-package: raster 创建）。使用函数“rasterToPolygons”我得到包含值“1”的所有栅格单元的多边形：

library(raster) 
dat = list()
dat$x = seq(1.5, by = 10, len = 10)
dat$y = seq(3.5, by = 10, len = 15)
dat$z = matrix(sample(c(0,1), size = 10*15, replace = T), 10, 15)

r=raster(dat);plot(r)

r_poly = rasterToPolygons(r, fun = function(r) {r == 1}, dissolve = F)
plot(r_poly, add = T)

我不使用“dissolve = T”来避免所有多边形合并成一个大多边形。相反，我希望获得一个新的 SpatialPolygonsDataFrame，其中合并了所有共享边或点的多边形。明确分离的多边形应可识别为单独的多边形。基于新的 SpatialPolygonsDataFrame 我想分析组合多边形的大小如下：

b = extract(r,r_poly_new) # "r_poly_new" contains the combined polygons
str(b)                    # list of clearly separated polygons
tab = lapply(b,table)      
tab

我的问题是双重的：1）如何组合共享边或点的多边形？2）如何将这些信息转换为允许分析组合多边形区域的格式？非常感谢您的反馈。

4

3 回答 3

10

您可以首先使用raster::clump()来识别连接的栅格单元的集群，然后应用于rasterToPolygons()“多边形化”这些单元。（但请注意，每个团块的面积可以直接从 计算RasterLayer而不将其转换为 a SpatialPolygonsDataFrame，如下所示）：

library(rgeos) ## For the function gArea

## Clump and polygonize
Rclus <- clump(r)
SPclus <- rasterToPolygons(Rclus, dissolve=TRUE)

## Check that this works
plot(SPclus, col = seq_along(SPclus))

## Get cluster areas from RasterLayer object
transform(data.frame(freq(Rclus)), 
          area = count*prod(res(Rclus)))

## Get cluster areas from SpatialPolygons object
transform(data.frame(SPclus), 
          area = gArea(SPclus, byid=TRUE))

在此处输入图像描述

于 2013-12-18T16:52:31.657 回答
3

该rgeos软件包有许多多边形操作工具。gUnion将结合在一起接触多边形：

require(rgeos)
uni <- gUnion( r_poly , r_poly )
plot( uni , col = 2 )

在此处输入图像描述

于 2013-12-18T13:24:04.250 回答
2

rasterToPolygons()是一个计算上非常昂贵的操作，因此，假设 CRS 是平面的，我会选择：

m <- clump(r)
f <- freq(m)
f[,2] <- f[,2] * xres(r) * yres(r)

对于经度/纬度，我会使用：

a <- area(r)
zonal(a, m, 'sum')
于 2013-12-18T17:41:03.853 回答