Процент перекрытия между SpatialPolygonsDataFrame - PullRequest
Новая
       46

Процент перекрытия между SpatialPolygonsDataFrame

0 голосов
/ 23 января 2019

Мне нужно рассчитать процент перекрытия между несколькими объектами типа SpatialPolygonsDataFrame в R. У меня они есть в виде шейп-файлов и загрузили выделение здесь . 

## 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 шейп-файлов для этого. Поэтому мне нужен более удобный метод. Мое идеальное решение выложено здесь . Но я борюсь с тем, чтобы адаптировать следующую строку кода из ответа на мои файлы: 

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

Как мне составить список моих файлов, присваивая переменные "a", "b", "c" и т. Д., Которые необходимы для следующих вычислений? Все, что мне нужно сделать, это создать объект типа poly из ответа с моими входными SpatialPolygonsDataFrame - или SpatialPolygons -файлами. Кто-нибудь знает, как это сделать? Большое спасибо!

1 Ответ

0 голосов
/ 24 января 2019

Вот адаптация этого ответа 

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)

И возьми оттуда.

...