What if you get a raster layer with number of people per raster cell, like for example the population layer from Afripop, and you want to convert it to a population density layer?

Well, 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.

I thought I had seen a function in GRASS GIS to do this, but that might have been a typical case of the wish being the father to the thought. But anyway, it isn’t terribly difficult to calculate it yourself using the map calculator.

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

Using this information, we can calculate a raster layer `rcs`

with the surface area per raster cell as (note that I am defining R in meters here):

PI=3.1415926536 RES=0.0008333 R=6371007.2 r.mapcalc "rcs = (sin(y()+$RES/2) - sin(y()-$RES/2)) * ($RES * $PI/180) * $R^2";

In the first three lines I define PI, the cell resolution and the earth radius respectively. I use these variables in the r.mapcalc expression (you can use those numbers directly in the expression if you like).

You can now see that for a raster with a resolution of 1/1200 × 1/1200 arc degrees, the area of a raster cell at the equator will be approximately 8585 m², but 8310 m² in northern Ethiopia (one should of course always realize that these are approximations).

Edit: just realized that 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 still need to set yourself is PI.

g.region rast=afripop10adj@PERMANENT 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;

r.mapcalc "afripopdens = afripop / (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)

Thank you for this information – it was very helpful! I have a small issue to resolve with the code however and am hoping you can provide some insight. In your Edit, the formula calculating the area relies on attributes coming from g.region and g.proj. The “$a” you use for the earth radius is supposed to come from g.proj however when I try it, I get the following attributes in the output: name, proj, no_defs, datum, ellps, towgs84, unit, units, and meters; no “a” … where or how did you get $a ??

It should have been g.proj -j (j instead of g flag). I corrected it in the text.

I once had the same problem, to calculate total annual Carbon and Nitrogen fixation for a paper (Fig.3 in http://dx.doi.org/10.1038/ngeo1486). Therefore I wrote my own function as a new Grass GIS command “r.area”. You can find it on GitHub (https://github.com/joergsteinkamp/r.area) and anybody is welcome to contribute and further improve it.

Thanks for sharing. To suggestions; why not adding your addon to the official GRASS GIS addon repository? And you might want to consider renaming the function as there is already a function r.area [https://grass.osgeo.org/grass70/manuals/addons/r.area.html].

Nice article by the way. I have to admit that I know very little about cryptogamic covers, but it provides in interesting read. I was particularly interested to learn about the relative importance of cryptogamic cover in the two tropical forests basins in terms of flux intensity of carbon net uptake (if I understand it correctly).

Hi, I wrote them to include it but sadly never got a reply. Maybe I should retry.

They (the active developers of GRASS GIS) are normally very responsive. Just ask write-access to the addons SVN repository via the GRASS-PSC mailing list (see https://trac.osgeo.org/grass/wiki/HowToContribute#WriteaccesstotheGRASS-Addons-SVNrepository). After you got the write-access, you can add the addon yourself (as said, you probably need to rename it first as there is already an addon under that name).