从 shapefile 计算百分比覆盖率

Posted

技术标签:

【中文标题】从 shapefile 计算百分比覆盖率【英文标题】:Calculating percent cover from shapefiles 【发布时间】:2014-01-09 15:53:22 【问题描述】:

我通常不使用 shapefile,所以在这里我有点迷失了。我有两个 shapefile,每个文件都有多个对象。第一个是一组 32 个多边形(每个都是一个图)。第二个 shapefile 有超过 10,000 个对象,它们代表每个地块内不同大小的植被簇。我想弄清楚。

1) 我如何计算每个站点内总植被覆盖率的百分比?

2) 每块地块面积小于5米的植被占百分之几?

这就是我的数据在 ArcGIS 中单个绘图的样子。

【问题讨论】:

通常,对于多边形文件,该区域嵌入在 shapefile 本身中。请将您的 shapefile 放到网上(Dropbox?)并发布链接作为对您问题的编辑。 此链接将带您访问我的保管箱中的 shapefile:dropbox.com/sh/wyokxximppexyb3/p7VC-pfF2E 另外三个问题(对不起......):(1)当我这样做时,我得到的植被覆盖率非常低(通常 @jlhoward; 0.2% 的覆盖率不加起来。单位均应以米为单位,面积均应以平方米为单位。面积是指每个植被对象的面积,以平方米为单位。 FID 绘图与绘图 shapefile 中的绘图 ID 不对应。感谢您在此过程中提供的帮助.. 好的。在两个 shapefile 的数据部分都有一列 LocCode。似乎每个绘图多边形都有一个LocCode。这是将植被对象与绘图多边形相关联的字段吗?当我使用它而不是 FID_plots 时,我的覆盖率在 13% 到 66% 之间。听起来对吗??顺便说一句:所有这些来回的原因是我不想发布一个愚蠢的答案。 【参考方案1】:

我想,下面的代码会做你想做的事。

注意:这使用存储在 shapefile 多边形中的区域信息(如下所述)。它确实使用植被 shapefile 数据部分中的 Area 列。在大多数情况下,您的 Area 与存储在 shapefile 中的区域相同,但在某些情况下,您的 Area 更大。由于我不知道您的 Area 数据来自哪里,因此使用与 shapefile 多边形一起存储的信息似乎更安全。

library(rgdal)
library(ggplot2)

setwd("<directory containing all your shapefiles>")
plt.map <- readOGR(dsn=".",layer="plots")
veg.map <- readOGR(dsn=".",layer="veg_in_plots")
# associate LocCode with polygon IDs
plt.data <- cbind(id=rownames(plt.map@data), LocCode=plt.map@data$LocCode)
veg.data <- cbind(id=rownames(veg.map@data), LocCode=veg.map@data$LocCode)
# function to extract area from polygon data
get.area <- function(polygon) 
  row <- data.frame(id=polygon@ID, area=polygon@area, stringsAsFactors=F)
  return(row)

# area of each plot polygon
plt.areas <- do.call(rbind,lapply(plt.map@polygons, get.area))
plt.data  <- merge(plt.data,plt.areas, by="id")  # append area column to plt.data
# area of each vegetation polygon
veg.areas <- do.call(rbind,lapply(veg.map@polygons, get.area))
veg.data  <- merge(veg.data,veg.areas, by="id")  # append area column to veg.data
# total area of vegetation polygons by LocCode
veg.smry  <- aggregate(area~LocCode,data=veg.data,sum)
smry      <- merge(plt.data,veg.smry,by="LocCode")
smry$coverage <- with(smry,100*area.y/area.x)    # coverage percentage
# total area for vegetation object with A < 5 msq
veg.lt5   <- aggregate(area~LocCode,data=veg.data[veg.data$area<5,],sum)
smry      <- merge(smry, veg.lt5, by="LocCode") 
# fraction of covered area coming from veg. obj. with A < 5 msq
smry$pct.lt5 <- with(smry, 100*area/area.y)

产生这个:

# head(smry)
#   LocCode id   area.x   area.y coverage      area  pct.lt5
# 1       1  3 1165.916 259.2306 22.23408  60.98971 23.52720
# 2      10 11 1242.770 366.3222 29.47626  88.21827 24.08216
# 3      11 12 1181.366 213.2105 18.04779 129.21612 60.60496
# 4      12 13 1265.352 577.6037 45.64767 236.83946 41.00380
# 5      13 14 1230.662 226.2686 18.38593  48.09509 21.25575
# 6      14 15 1274.538 252.0577 19.77640  46.94874 18.62619

说明:

可以使用rgdal 包中的readOGR(...) 将形状文件导入R。导入多边形 shapefile 时,结果是“SpatialPolygonDataFrame”对象。这些对象基本上有两个部分:一个多边形部分,其中包含绘制每个多边形所需的坐标,以及一个数据部分,其中包含每个多边形的数据(因此,每个多边形一行)。如果 shapefile 被导入为,例如,map

map <- readOGR(dsn=".",layer="myShapeFile")

然后可以以map@polygonmap@data 访问多边形和数据部分。事实证明,多边形区域存储在多边形部分中。为了获取区域,我们定义了一个函数get.area(...),它从多边形中提取区域和多边形 ID。然后我们使用lapply(...) 为所有多边形调用该函数,并使用rbind(...) 将所有返回值逐行绑定在一起:

plt.areas <- do.call(rbind,lapply(plt.map@polygons, get.area))
veg.areas <- do.call(rbind,lapply(veg.map@polygons, get.area))

现在我们需要将植被区域与绘图多边形相关联。这是通过列 LocCode 完成的,该列存在于每个 shapefile 的数据部分中。因此,我们首先将多边形 ID 与地块和植被区域的 LocCode 相关联:

plt.data <- cbind(id=rownames(plt.map@data), LocCode=plt.map@data$LocCode)
veg.data <- cbind(id=rownames(veg.map@data), LocCode=veg.map@data$LocCode)

然后我们根据多边形ID追加面积列:

plt.data  <- merge(plt.data,plt.areas, by="id")  # append area column to plt.data
veg.data  <- merge(veg.data,veg.areas, by="id")  # append area column to veg.data

然后我们需要通过LocCode对植被面积求和:

veg.smry  <- aggregate(area~LocCode,data=veg.data,sum)

最后将其与绘图多边形区域合并:

smry      <- merge(plt.data,veg.smry,by="LocCode")

smry 数据框中,area.x 是地块的面积,area.y 是该地块中植被覆盖的总面积。因为,对于这两个 shapefile,投影是:

 +proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

单位以米为单位,面积以 msq 为单位。为了确定有多少植被来自 smry 合并:

veg.lt5   <- aggregate(area~LocCode,data=veg.data[veg.data$area<5,],sum)
smry      <- merge(smry, veg.lt5, by="LocCode") 

最后,利用我们拥有的数据,可以直接为每个绘图区域渲染地图:

cols   <- c("id","LocCode")
plt.df <- fortify(plt.map)
plt.df <- merge(plt.df,plt.data[cols],by="id")
veg.df <- fortify(veg.map)
veg.df <- merge(veg.df,veg.data[cols],by="id")
ggp <- ggplot(plt.df, aes(x=long, y=lat, group=group))
ggp <- ggp + geom_path()
ggp <- ggp + geom_polygon(data=veg.df, fill="green")
ggp <- ggp + facet_wrap(~LocCode,scales="free")
ggp <- ggp + theme(axis.text=element_blank())
ggp <- ggp + labs(x="",y="")
ggp <- ggp + coord_fixed()
ggp

【讨论】:

以上是关于从 shapefile 计算百分比覆盖率的主要内容,如果未能解决你的问题,请参考以下文章

如何计算被另一个多边形覆盖的多边形的百分比

使用 Jenkins 从 Karma.js 获得百分比覆盖率

从以前的记录计算百分比

数据仓库事实表的设计

从多列值计算百分比

使用 FxCop 获得百分比代码覆盖率