如何从owin格式二值掩码中计算聚类区域面积？

问题背景 #

手里有个owin格式的二值掩码，包含5个不同大小的聚类点，需要计算每个聚类的面积。原以为spatstat包无对应功能，便转成Raster图层用igraph的clump函数识别聚类，但计算面积时遇到警告和错误。

尝试过的操作与问题 #

1. 转换图层并识别聚类
   执行如下代码后，能画出分色的5个聚类，但调用area(clusters)时收到警告：

   library(igraph)
   raster_mask <- raster(as.im(dd))
   clusters <- clump(raster_mask)
   

   This function is only useful for Raster* objects with a longitude/latitude coordinates

2. 设置CRS后警告依旧
   查看crs(clusters)返回NA，设置CRS：

   crs(clusters) <- "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs +type=crs"
   

   此时图层信息如下：

   > print(areas)
   class      : RasterLayer 
   dimensions : 135, 129, 17415  (nrow, ncol, ncell)
   resolution : 0.9790567, 1.224028  (x, y)
   extent     : 490715.5, 490841.8, 5429337, 5429502  (xmin, xmax, ymin, ymax)
   crs        : +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs +type=crs 
   source     : memory
   names      : layer 
   values     : 1.198393, 1.198393  (min, max)
   

   但调用area(clusters)仍弹出相同警告。

3. 转换经纬度CRS时出错
   尝试在聚类前转换CRS：

   crs(raster_mask) <- "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs +type=crs"
   raster_mask <- projectExtent(raster_mask, "+proj=longlat +datum=WGS84")
   plot(raster_mask)
   

   报错：

   Error in .plotraster2(x, col = col, maxpixels = maxpixels, add = add, :
   no values associated with this RasterLayer

解决方案 #

方案1：直接用spatstat原生方法（推荐） #

spatstat其实自带识别聚类并计算面积的功能，无需转Raster图层：

library(spatstat)
# 识别owin中的聚类，返回多个子owin组成的列表
clust_owin <- connected(dd)
# 计算每个聚类的面积
cluster_areas <- sapply(clust_owin, area)

结果的单位与原owin的坐标单位一致（你的数据是UTM米单位，面积即为平方米）。

方案2：基于Raster图层的解决方法 #

因为你的CRS是UTM平面坐标系，area()函数是针对经纬度椭球计算面积的，所以会触发警告。直接通过像素数计算面积更准确：

# 计算单个像素的面积
pixel_area <- res(clusters)[1] * res(clusters)[2]
# 统计每个聚类的像素数，再乘以单个像素面积
cluster_stats <- freq(clusters)
cluster_stats$area <- cluster_stats$count * pixel_area

最终cluster_stats会包含每个聚类的ID、像素数和对应的面积。

最小可复现示例制作方法 #

用spatstat生成模拟的owin格式二值掩码，即可复现场景：

library(spatstat)
# 创建矩形窗口
win <- owin(xrange = c(0, 100), yrange = c(0, 100))
# 生成含5个聚类的二值掩码
set.seed(123)
pp <- rpoispp(5, win = win)
mask <- as.mask(pp, dimyx = c(135, 129))
# mask即为带5个聚类的owin格式二值掩码，可用于测试

内容的提问来源于stack exchange，提问作者Qiyuan