[R-sig-Geo] Calculating areas on Earth’s surface
Karl Ove Hufthammer
karl at huftis.org
Thu Jul 21 13:18:44 CEST 2011
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
More information about the R-sig-Geo
mailing list