Wrong area calculation with geosphere:: areaPolygon?

Question

I have used geosphere::areaPolygon() with success many times but now I have encountered an issue.

I have a polygon sp2 containing another polygon sp1. So sp2 should have larger area than sp1. When I calculate each area with areaPolygon(), I get the opposite result.

areasp1 = 10133977
areasp2 = 9858811

I used gSymdifference to find the symmetrical different polygon sp3, which has

areasp3 = 275165.4

sp1 : red dashed line
sp2 : black line
sp3 : green dotted line

Code I used:

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 

All polygons were tested for self-intersections or other issues with gIsValid(), so it has nothing to do with the problem mentioned here.

Any idea on why sp1 has larger area than sp2?

Note: If you sum areasp2 + areasp3, it almost equals areasp1. I don't know if this is by chance.

Answer 1

This is nice illustration of how things aren't always what they seem when you treat angular coordinates as if they were planar. Your plot assumes straight line connections between nodes. That would be reasonable if the coordinates were planar. But they are not, and that makes all the difference in this case. Below I redraw the plot with the raw data and after adding vertices following the shortest distance between the nodes. You can see that object 2 is much smaller than it seems.

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))

You can perhaps better appreciate this by projecting the coordinates to a planar reference system.

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)

The shortest route between Isfahan and Taipei does not go across Nepal.