Header

Calculating the raster cell area of an unprojected raster layer

Author: Paulo van Breugel
Updated on: 2017-05-30

Issue

From Worldpop you can download raster layers with the number of people per raster cell. But what if you need the population density, e.g., the number of people per km2. Obviously, you need to divide the number of people by the surface area of the raster cells. However, the surface area of the raster cells of an unprojected (lat/lon) are not constant; they decrease with increasing latitude. So what you need is a raster layer with the surface areas of the cells.

Solution

In the latest (development) version of GRASS GIS, there is the area() function in r.mapcalc which does exactly that, it computes the area of the current cell in square meters. For earlier versions, including the current stable version of GRASS GIS, you'll need to compute the area yourself.

The formulate I am using below comes from an Q & A on the GIS StackExchange forum. There, you can read that the area of a cell spanning longitudes X0 toX1 (X1 > X0) and latitudes Y0 to Y1 (Y1 > Y0) equals;

(sin(Y1) - sin(Y0)) × (X1 - X0) × R²

where:

  • X0 and X1 are expressed in radians.
  • X1 – X0 is calculated modulo 2×pi (e.g., -179 – 181 = 2 degrees, not -362 degrees).
  • R is the authalic Earth radius

If you set the region based on your raster, you can then get the actual horizontal and vertical cell resolutions of your raster layer using g.region. Moreover, you can get the earth radius assumed for your spatial reference system using g.proj. Thus the only thing you need to set yourself is PI. Below, I want to compute the area of the raster cells of the raster layer 'humpop'.

g.region rast=humpop
PI=3.1415926536
export `g.region -g`
export `g.proj -j`
r.mapcalc "rcs = (sin(y()+$nsres/2) - sin(y()-$nsres/2)) * ($ewres * $PI/180) * float($a)^2"

Now, to get the population density (people / km2), you simply do:

r.mapcalc "humpop = humpop / (rcs/1000000)"

I am not sure QGIS has an equivalent for the y() and x() parameters in the raster calculator. But you can always use the GRASS toolbox. You probably can also use the raster calculators or their equivalents in other GIS packages. Or perhaps they have a dedicated function, like in the raster package in R (the area function)