SpatialPolygonsDataFrame 之间的重叠百分比

Question

我需要计算 R 中 SpatialPolygonsDataFrame 类型的多个对象之间的重叠百分比。我将它们作为 shapefile 并上传了一个选择 here。

## load libraries
library(raster)
library(rgdal)
library(rgeos)
library(maptools)
library(sp)
library(FRK)

setwd("...")
Polar1 <- shapefile("Din_biome_Rock_Ice.shp")
Polar2 <- shapefile("FAO_GEZ_Polar.shp")
Polar3 <- shapefile("KG_Beck_EF.shp")

看看：

> Polar1
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -179.9999, 180, -89.89197, 82.3525  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : dummy 
value       :     0

我可以将它们转换为 SpatialPolygons:

T_Polar1 <- as(Polar1, "SpatialPolygons")

如我们所见...

> T_Polar1
class       : SpatialPolygons
features    : 1 
extent      : -179.9999, 180, -89.89197, 82.3525  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

无论哪种方式，我都可以像那样简单地计算百分比：

intersection_Polar1_Polar2 <- gIntersection(Polar1, Polar2)
round(area(intersection_Polar1_Polar2)/area(Polar1)*100, 2)

但是我有 200 多个 shapefile 可以这样做。因此，我需要一个更方便的方法。我理想的解决方案已发布 here。但是我正在努力的是将以下代码行从我的文件的答案中改编出来：

poly <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "a"),Polygons(list(Polygon(p2)), "b"),Polygons(list(Polygon(p3)), "c")),1L:3L)

我怎样才能列出我的文件，分配以下计算所需的变量 "a"、"b"、"c" 等等？我需要做的就是用我输入的 SpatialPolygonsDataFrame- 或 SpatialPolygons- 文件从答案创建一个像 poly 这样的对象。有人知道怎么做这个吗？非常感谢！

Answer 1

这里是对那个答案的改编

library(raster)
library(sp)
## example data:
p1 <- structure(c(0, 0, 0.4, 0.4, 0, 0.6, 0.6, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
p2 <- structure(c(0.2, 0.2, 0.6, 0.6, 0, 0.4, 0.4, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
p3 <- structure(c(0, 0, 0.8, 0.8, 0, 0.8, 0.8, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
poly <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "a"),Polygons(list(Polygon(p2)), "b"),Polygons(list(Polygon(p3)), "c")),1L:3L)
plot(poly)
crs(poly) <- "+proj=utm +zone=1, +datum=WGS84"

我了解到您有单独的 SpatialPolygonDataFrame 对象，每个对象都有一个（多个）多边形。您可以将这些存储在列表中。使用示例数据：

ss <- list(poly[1,], poly[2,], poly[3,])

然后做这样的事情：

n <- length(ss)
overlap <- matrix(0, nrow=n, ncol=n)
diag(overlap) <- 1

for (i in 1:n) {
    ss[[i]]$area <- area(ss[[i]])
}


for (i in 1:(n-1)) {
    for (j in (i+1):n) {
        x <- intersect(ss[[i]], ss[[j]])
        if (!is.null(i)) {
            a <- area(x)
            overlap[i,j] <- a / ss[[i]]$area
            overlap[j,i] <- a / ss[[j]]$area
        } 
    }
}

根据您的数据，您可以制作一个多边形列表，如下所示：

ff <- c("Din_biome_Rock_Ice.shp", "FAO_GEZ_Polar.shp", "KG_Beck_EF.shp")
ss <- lapply(ff, shapefile)

然后从那里拿走。

SpatialPolygonsDataFrame 之间的重叠百分比

Percentage of overlap between SpatialPolygonsDataFrame

r

polygon

shapefile

overlap

intersect