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@polygons, function(x) x@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

_______________________________________________

R-sig-Geo mailing list

[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo