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

[Karl Ove Hufthammer]

Matthew Vavrek wrote:

> have you looked at the PBSmapping package? It has a calcArea() function
> that might do what you need it to do. However, you need to get your data
> converted to a PolySet rather than a SpatialPolygons object.

Thanks for the suggestion. It seems that PBSmapping uses an UTM 
transformation based on the location of the (entire) PolySet.
The accuracy is thus not good enough (e.g., if you want to see
the distortion caused *by* a standard UTM projection), but it gave
me an idea:

For each polygon / grid cell, create an projection tailored to this 
exact polygon, for example a non-standard transverse Mercator projection.
Use this to calculate the area of the polygon. If the polygon isn’t too 
big, this should work well.

Here’s some example code, that aims to visually present the amount of
area distortion caused by an UTM zone 33 projection. Note the use
of +k=1.0 to make the distortion in the transverse Mercator projection
be ~0 in the *middle* of each polygon.

#########################################################################
# Various packages used
library(rgdal)
library(mapdata)
library(maptools)

# Bounding box
xl=c(4,30)
yl=c(57,72)

unproj=CRS("+init=epsg:4326") # Unprojected (longlat)
proj=CRS("+init=epsg:32633")  # UTM zone 33

# Fetch a simple map for overlaying
wmap=map("worldHires", c("Norway","Sweden","Denmark"), xlim=xl, ylim=yl, fill=TRUE, plot=FALSE, resolution=0)
wmap.sp=map2SpatialPolygons(wmap, wmap$names, proj4string=unproj)
wmap.sp.proj=spTransform(wmap.sp, proj)

# Create a grid 
bb=bbox(wmap.sp.proj)
grid.lowerleft=floor(c(bb[1,1],bb[2,1]))
cell.s=50*c(1e3,1e3)
cell.d=ceiling((c(diff(bb[1,]),diff(bb[2,])))/cell.s)
cells=as.SpatialPolygons.GridTopology(GridTopology(grid.lowerleft, cell.s, cell.d))
proj4string(cells)=proj

# Transform the grid cells to longlat, so that we can extract their
# centroid longlat coordinates for use in determining an appropriate projection
cells.ll=spTransform(cells, unproj)
areas=numeric(length(cells.ll)) # Vector to hold the cell areas
for(i in seq_along(cells.ll)) {
  cell=cells.ll[i,]
  pos=coordinates(cell)
  tproj=CRS(paste("+proj=tmerc +lat_0=", pos[2]," +lon_0=",pos[1], " +k=1.0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0",sep=""))
  areas[i]=spTransform(cell, tproj)@polygons[[1]]@area/1e6
}

# Convert the grid cells to a SpatialPolygonsDataFrame
# holding the areas of each cell
cells=as(cells, "SpatialPolygonsDataFrame")
cells$dummy=areas
cell.area.exp=prod(cell.s)/1e6
maxdiff=ceiling(max(abs(areas-cell.area.exp)))
spplot(cells, sp.layout = list(list("sp.polygons", wmap.sp.proj,lw=2)),
              at=seq(cell.area.exp-maxdiff,cell.area.exp+maxdiff, length=18))
#########################################################################

The result does at least look *reasonable*, with the least amount of 
distortion not in the centre of the UTM projection, but in two bands
a few degrees east/west of the centre.

I wouldn’t trust the results too much, though, as other choices of
projections (e.g., ‘omerc’ instead of ‘tmerc’) gives other results.
I’m not sure which projection to use to ensure optimal results. 

-- 
Karl Ove Hufthammer