Calculating areas on Earth’s surface

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&lt;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