compute surface area for lat-long grids.

areaLL.R

areaLL <- function(r){
    ## cell areas in km^2
    areas = area(r)

    ## area of non-NA in km^2
    ## this is the flat area
    a = sum(values(area(r)*(!is.na(r))))
    centre = cellFromRowCol(r, nrow(r)/2, ncol(r)/2)
    ## compute w/h of middle cell 
    a1 = areas[centre]
    csize = sqrt(a1*1000*1000)
    corners = xyFromCell(r,cellFromRowColCombine(r,(0:1)+nrow(r)/2,(0:1)+ncol(r)/2))
    wh=1000*spDists(corners[1,,drop=FALSE], corners[2:3,],longlat=TRUE)
   
    ## now find surface area assuming all cells like the central cell@
    m = raster:::as.matrix(r)
    sa = surfaceArea(m, cellx=wh[1], celly=wh[2], byCell=TRUE)
    uncorrected = sum(sa,na.rm=TRUE)/(1000*1000)

    ## correct area by scaling cell surface areas by relative cell areas:
    corrected = sum(sa*(raster:::as.matrix(areas)/a1),na.rm=TRUE)/(1000*1000)
    list(flat=a, uncorrected=uncorrected, corrected=corrected)
}