Как вычислить площадь многоугольника в R - PullRequest
Новая
       10

Как вычислить площадь многоугольника в R

0 голосов
/ 18 сентября 2018

У меня проблемы с вычислением в R области импортированного шейп-файла, который имеет многоугольник, состоящий из нескольких частей (один объект, содержащий два отдельных многоугольника).Я заметил, что ArcMap дал мне другое значение для области шейп-файла, чем raster::area.Чтобы выяснить, какая программа дала мне правильную область, я разбил шейп-файл на отдельные части и пересчитал площадь двух отдельных многоугольников: 

library(raster)

> single_part <- shapefile("../Desktop/test/test_sp.shp")
> area(single_part)
[1] 575924.0 433409.8
> sum(area(single_part))
[1] 1009334
> 
> multi_part <-  shapefile("../Desktop/test/test_mp.shp")
> area(multi_part)
[1] 1018390

Теперь я понимаю, что знаю об этой проблеме, мне следуетвсегда разбивать классы объектов полигонов на отдельные части, но кто-нибудь знает, как raster::area вычисляет площадь многоугольников?Я также пытался использовать rgeos::gArea, но получил тот же результат.Есть ли способ расчета площади многоугольников в R?

Мне бы очень хотелось знать, потому что они довольно распространены, и я пытаюсь перейти от всех моих анализов в ArcMap к R.

В случае, если это полезно, вот изображениешейп-файла: многокомпонентный поли-шейп-файл

РЕДАКТИРОВАТЬ ДОБАВЛЕНО 21/9/2018 -------------------------------------------------------

Вот ссылка на шейп-файл test_mp.shp

Из того, что я могу сказать, кажется, что проблема связана с тем, как R (против ArcMap) интерпретирует дыры.См. Разницу между дисплеем ArcMap и R дисплеем .По какой-то причине R заполняет эти отверстия как часть шейп-файла, что должно быть причиной того, что я получаю различные расчеты для этой области.Что-то не так с шейп-файлом или как я его импортирую?

1 Ответ

0 голосов
/ 18 сентября 2018

Очевидно, что ваш объект с именем 'multi_part' имеет только один (multi?) Многоугольник, так как area возвращает одно значение.Здесь я проиллюстрирую, как исследовать то, что вам нужно: 

library(raster)
d <- getData('GADM', country='Isle of Man', level=0)

area(d)
[1] 579672897

Разделить на 4 полигона (острова) 

dd <- disaggregate(d)
a <- area(dd)
a
[1]     19424.12   2705442.41     25629.79 576922400.90
sum(a)
[1] 579672897

Одна и та же область, и нет никаких причин, по которым они будут другими,За исключением, возможно, если есть путаница с отверстиями многоугольника.Трудно комментировать без ваших данных.

Вы можете записать эти объекты на диск (см. Ниже) и посмотреть, что ArcGIS предоставляет вам как область (но обратите внимание, что в этом примере используются координаты lon / lat, я не уверен, что ArcGIS может вычислить области по ним).

shapefile(d, "man.shp")

Вот случай с и без отверстия: 

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))

# two (overlapping) polygons (no hole)   
pol1 <- spPolygons(p1, p2, crs="+proj=utm +zone=1 +datum=WGS84")
# single polygon with hole
pol2 <- spPolygons(list(p1, p2), crs="+proj=utm +zone=1 +datum=WGS84")

a <- area(pol1) / 10e+9
b <- area(pol2) / 10e+9

a
#[1] 10925   800
sum(a)
#[1] 11725
a[1]-a[2]
#[1] 10125

b соответствует a[1] - a[2], как и ожидалось 

b
#[1] 10125

Я получаю точноте же результаты с ArcGIS, используя «вычислить геометрию» для поля в таблицах атрибутов.

...