从形状文件计算覆盖百分比

Question

我通常不使用shapefile，所以我在这里有点失落。我有两个shapefile每个都有多个对象。首先是一组32个多边形（每个是一个图）。第二个shapefile具有> 10,000个对象，这些对象代表每个绘图内不同大小的植被簇。我试图弄清楚。

1）如何计算每个地点的总植被覆盖率？

2）每块地块每个植被覆盖面积的百分比小于5米？

这就是我的数据在ArcGIS中对于单个图的外观。

Answer 1

下面的代码会做你想要什么，我想。

注意：这使用存储在shapefile多边形中的区域信息（如下所述）。它确实而不是在植被shapefile数据部分使用Area列。在大多数情况下，您的Area与存储在shapefile中的区域相同，但在某些情况下，您的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”对象。这些对象基本上有两部分：一个多边形部分，它具有绘制每个多边形所需的坐标;以及一个数据部分，其中包含每个多边形的数据（所以每个多边形只有一行）。如果shape文件导入为，例如，map，

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

然后多边形和数据部分可作为map@polygon和map@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是该地块植被覆盖的总面积。因为对于两种形状文件，投影都是：

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

单位是米，面积是msq。要确定有多少植被，从区域来< 5 MSQ，我们总用面积< 5植被区和结果合并起来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