Skip to content

Calculate area of raster cells

5 messages · Tim Erbrecht, Paul Hiemstra, Roger Bivand +1 more

#
Tim Erbrecht wrote:
with the meuse dataset:

library(sp)
data(meuse.grid)
gridded(meuse.grid) = ~x+y
fullgrid(meuse.grid) = TRUE
bla = gridparameters(meuse.grid)
area = bla$cellsize[1] * bla$cellsize[2]
area

cheers,
Paul
#
On Wed, 28 Apr 2010, Paul Hiemstra wrote:

            
These are projected, not geographical coordinates, unfortunately. The case 
for calculating areas on an ellipsoid is much harder, I'm afraid.

Roger

  
    
#
Tim,

if your SpatialGridDataFrame is called 'spgdf', you can do

library(raster)
r = raster(spgdf)
a = area(r)
plot(a)
spgdf2 = as(a, 'SpatialGridDataFrame')

That gives a reasonable estimate in most cases (it uses the
longitudinal span of the center of each cell).

For more precise computations,  you can coerce the grid (or perhaps
the first cell of each row if the grid is large) to polygons and then
use "areaPolygons" in the development version (on R-Forge) of
"geosphere".
install.packages("geosphere", repos="http://R-Forge.R-project.org")

This function computes the area of spherical polygons.  How you would
implement this could very well depend on the size of your raster

rr = raster(r)
p = rasterToPolygons(rr)   # or an sp coerce method from the spgdf
rr[] = areaPolygon(p)

Here is a short illustration:
r = raster(nrow=30, ncol=4, ymn=-89)  # excluding South Pole; bugs there
p = rasterToPolygons(r)
library(geosphere)
r[] = areaPolygon(p)
plot(r)


Robert
On Wed, Apr 28, 2010 at 3:54 AM, Tim Erbrecht <timerbrecht at web.de> wrote:
#
P.S. The solution I just posted assumes the earth is a sphere, not an
ellipsoid (which would be more accurate).

        
On Wed, Apr 28, 2010 at 4:16 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote: