Неправильный расчет площади с геосферой:: areaPolygon?

Question

Я использовал geosphere::areaPolygon() с успехом много раз, но теперь я столкнулся с проблемой.

У меня есть многоугольник sp2, содержащий другой многоугольник sp1. Таким образом, sp2 должен иметь большую площадь, чем sp1. Когда я рассчитываю каждую область с areaPolygon()Я получаю противоположный результат.

areasp1 = 10133977
areasp2 = 9858811

Я использовал gSymdifference, чтобы найти симметричный другой многоугольник sp3, который имеет

areasp3 = 275165.4

sp1: красная пунктирная линия
sp2: черная линия
sp3: зеленая пунктирная линия

Код, который я использовал:

library(sp)
library(geosphere)
library(rgeos)
library(raster)

s1 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 84.9176, 60.40123, 
58.97929, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 28.25964, 
36.89905, 37.72305, 52.58793, 43.47459))

s2 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 54.35377, 58.55841, 
94.95237),lat= c(30.25074, 29.83075, 23.7992, 31.61798, 52.58793, 43.47459))

sp1 <- SpatialPolygons(list(Polygons(list(Polygon(s1)),1)))
crs(sp1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

sp2 <- SpatialPolygons(list(Polygons(list(Polygon(s2)),1)))
crs(sp2) <-CRS ("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

areasp1 <- areaPolygon(sp1)/10^6    # to get area in square km
areasp2 <- areaPolygon(sp2)/10^6

sp3 <- gSymdifference(sp1,sp2,checkvalidity = TRUE)
areasp3 <- areaPolygon(sp3)/10^6 

Все полигоны были проверены на самопересечение или другие проблемы с gIsValid()Таким образом, это не имеет ничего общего с проблемой, упомянутой здесь.

Есть идеи, почему sp1 имеет большую площадь, чем sp2?

Примечание: если вы суммируете areasp2 + areasp3, это почти равно areasp1, Я не знаю, случайно ли это.

Answer 1

Это хорошая иллюстрация того, как вещи не всегда выглядят такими, какими они кажутся, когда вы рассматриваете угловые координаты, как если бы они были плоскими. Ваш график предполагает прямую связь между узлами. Это было бы разумно, если бы координаты были плоскими. Но это не так, и в этом случае все имеет значение. Ниже я перерисовываю график с необработанными данными и после добавления вершин, следующих за кратчайшим расстоянием между узлами. Вы можете видеть, что объект 2 намного меньше, чем кажется.

library(raster)
library(geosphere)

s1 <- cbind(lon= c(130.4327, 129.85127, 121.00775, 84.9176, 60.40123, 58.97929, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 28.25964, 36.89905, 37.72305, 52.58793, 43.47459))
s2 <- cbind(lon= c(130.4327, 129.85127, 121.00775, 54.35377, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 31.61798, 52.58793, 43.47459))
sp1 <- spPolygons(s1, crs="+proj=longlat +datum=WGS84")
sp2 <- spPolygons(s2, crs="+proj=longlat +datum=WGS84")

# add nodes on shortest path
mp1 <- makePoly(sp1, interval=100000)
mp2 <- makePoly(sp2, interval=100000)

# background countries for map
library(maptools)
data(wrld_simpl)
w <- crop(wrld_simpl, extent(mp2) + 40)

plot(w, col='light gray')
points(s2, pch=20, cex=3)
points(s1, pch=20, cex=2, col='gray')

lines(sp2, col='red', lwd=3)
lines(sp1, col='blue', lwd=1)
lines(mp2, col='red', lwd=4, lty=2)
lines(mp1, col='blue', lwd=2, lty=2)

legend("bottomright", c('sp1', 'mp1', 'sp2', 'mp2'), col=rep(c('blue', 'red'), each=2), lwd=rep(c(2,3), each=2), lty=c(1, 2))

Возможно, вы сможете лучше оценить это, спроецировав координаты на плоскую систему отсчета.

library(rgdal)
crs <- CRS("+proj=ortho +lat_0=40 +lon_0=90 +a=6370997 +b=6370997 +units=m")
ww <- spTransform(w, crs)
sp1x <- spTransform(sp1, crs) 
sp2x <- spTransform(sp2, crs) 
mp1x <- spTransform(mp1, crs) 
mp2x <- spTransform(mp2, crs) 

par(mai=c(0,0,0,0))
plot(ww, col='light gray')
lines(sp2x, col='red', lwd=3)
lines(sp1x, col='blue', lwd=1)
lines(mp2x, col='red', lwd=4, lty=2)
lines(mp1x, col='blue', lwd=2, lty=2)

Кратчайший маршрут между Исфаханом и Тайбэем не идет через Непал.