[R-sig-Geo] Calculating areas on Earth’s surface

[Karl Ove Hufthammer]

Calculating areas on Earth’s surface

I’m interested in the distortion in apparent area from various projections,
and need to calculate the *real* area of various polygons (or grid cells).
In an earlier post on this list, it was suggested to use a equal-area
projection. This partly works, but the accuracy is not good enough for
measuring the distortion for my purpose (I guess an accuracy < 0.05% for
moderately large areas would be OK).
 
Here’s a simple example:

library(sp)
library(rgdal)

unproj=CRS("+init=epsg:4326")                  # Unprojected, i.e. longitudes and latitudes
xy=GridTopology(c(0,55),c(8, 4), c(5,5))       # Small grid (covering Norway)
xy=as.SpatialPolygons.GridTopology(xy, unproj) # Grid as polygons

# Calculate the area of a polygons using the 'proj' projection
calcArea=function(xy, proj) {
  xy.trans=spTransform(xy, proj)
  sapply( xy.trans at polygons, function(x) x at area)/1e6
}

# Various equal-area projections
proj1=CRS("+proj=moll +lat_0=65 +lon_0=10") # Mollweide
proj2=CRS("+proj=sinu +lat_0=65 +lon_0=10") # Sinusoidal
proj3=CRS("+proj=tcea +lat_0=65 +lon_0=10") # Transverse Cylindrical
proj4=CRS("+init=epsg:32633")               # UTM 33N (*not* equal area)

# Areas according to the different projections all differ
a1=calcArea(xy, proj1)
a2=calcArea(xy, proj2)
a3=calcArea(xy, proj3)
a4=calcArea(xy, proj4)

# All areas calculated using Mollweide are small than the areas
# calculated using the sinusoidal projection:
all(a1<a2)

Since the areas from the various projections all differ, I’m not sure which 
one is most accurate. Does anyone know of a ‘real’ area function available 
somewhere in R, which can calculate the area of simple polygons on the WGS 84
ellipsoid (or perhaps even a spherical approximation would be good enough).
It would be very useful for drawing maps showing the distortion in area for
various projections.

-- 
Karl Ove Hufthammer