[GRASS-user] raster area calculations in lat/long

Dear List,

I'm a new grass user and have been learning on the command line. I have a raster file of the world in which the cells are floating point values of the percentage of a crop type within that cell (ex., 0.01, 0.29, 0.45). It's in Lat/Long. I want to calculate the total area of crop in specific regions that are defined in a separate raster file. I want to multiply the value of the cell by the cell area and then sum the cell values in each region to get total acreage. My questions are: 1) Can I do this in r.stats, r.report or r.coin, or do I need to use r.mapcalc?, and 2) Because my files are in lat/long, the area covered by each cell changes based on its latitude. Can this be accounted for in trying to come with total area?

Thanks in advance for any suggestions to a grass newbie,

Jerry Griffith

Hi Jerry,

I'm a new grass user and have been learning on the command line.
I have a raster file of the world in which the cells are floating
point values of the percentage of a crop type within that cell (ex.,
0.01, 0.29, 0.45). It's in Lat/Long. I want to calculate the
total area of crop in specific regions that are defined in a separate
raster file. I want to multiply the value of the cell by the cell
area and then sum the cell values in each region to get total acreage.
My questions are: 1) Can I do this in r.stats, r.report or r.coin, or
do I need to use r.mapcalc?,

see also r.statistics, r.sum, r.univar, r.reclass(.area), ....
I'm not enough of an expert with using the cover map functions to answer this well. Regardless r.mapcalc should be flexible enough to do what you like as a worst case scenario.
sounds like a fairly standard & straight forward task though.

and 2) Because my files are in lat/long, the area covered by each
cell changes based on its latitude. Can this be accounted for
in trying to come with total area?

AFAICT from a quick look at the code, r.report and r.coin are just
friendly front-ends to r.stats. So the question reduces to: does
r.stats account for changing cell area with latitude? -- It does.

specifically r.stats calls these libgis functions:
G_begin_cell_area_calculations()
G_area_of_cell_at_row()

so it recalculates the cell area for each cell row for lat/lon.
for more detail read the code comments for those functions in:
http://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/area.c
(or the GRASS programmer's manual)

If you like, reproject the maps to an equal-area map projection and double check that the results are similar.

if you were really keen you could try to write your own module in C or with the swig/python interface for the calc. But the existing modules should give you all the tools you need from the command line.

Hamish

Hi,

Ok, since GRASS does account for differing cell size in lat/long raster files, what I am after is mapcalc is to do something like:

croparea = [cellsize] * corn_map (my corn_map has a floating point value, ex. 0.29 that gives the fraction of the cell that is in corn cultivation)

I have tried simply typing in for [cellsize] the functions listed below, but that doesn't work or my syntax is wrong.

G_begin_cell_area_calculations()
G_area_of_cell_at_row()

and I haven't found in r.stats or related commands that allows an intermediate multiplication step. I have tried using r.proj (and with -n) to reproject to an equal area projection, but my crop maps are global coverages and the other Locations in lambert equal area are for the US, so the result is my global crop maps get clipped to the US.

Thanks very much,

Jerry

Hamish wrote:

Hi Jerry,

I'm a new grass user and have been learning on the command line.
I have a raster file of the world in which the cells are floating
point values of the percentage of a crop type within that cell (ex.,
0.01, 0.29, 0.45). It's in Lat/Long. I want to calculate the
total area of crop in specific regions that are defined in a separate
raster file. I want to multiply the value of the cell by the cell
area and then sum the cell values in each region to get total acreage.
My questions are: 1) Can I do this in r.stats, r.report or r.coin, or
do I need to use r.mapcalc?,
    
see also r.statistics, r.sum, r.univar, r.reclass(.area), ....
I'm not enough of an expert with using the cover map functions to answer this well. Regardless r.mapcalc should be flexible enough to do what you like as a worst case scenario.
sounds like a fairly standard & straight forward task though.

and 2) Because my files are in lat/long, the area covered by each

cell changes based on its latitude. Can this be accounted for
in trying to come with total area?
    
AFAICT from a quick look at the code, r.report and r.coin are just
friendly front-ends to r.stats. So the question reduces to: does
r.stats account for changing cell area with latitude? -- It does.

specifically r.stats calls these libgis functions:
G_begin_cell_area_calculations()
G_area_of_cell_at_row()

so it recalculates the cell area for each cell row for lat/lon.
for more detail read the code comments for those functions in:
http://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/area.c
(or the GRASS programmer's manual)

Hamish

Jerry Griffith wrote:

Ok, since GRASS does account for differing cell size in lat/long raster
files, what I am after is mapcalc is to do something like:

croparea = [cellsize] * corn_map (my corn_map has a floating point
value, ex. 0.29 that gives the fraction of the cell that is in corn
cultivation)

I have tried simply typing in for [cellsize] the functions listed below,
but that doesn't work or my syntax is wrong.

G_begin_cell_area_calculations()
G_area_of_cell_at_row()

These are C functions; they don't exist in r.mapcalc.

and I haven't found in r.stats or related commands that allows an
intermediate multiplication step. I have tried using r.proj (and with
-n) to reproject to an equal area projection, but my crop maps are
global coverages and the other Locations in lambert equal area are for
the US, so the result is my global crop maps get clipped to the US.

If you can't use a global equal-area projection, you could just use
the underlying equations within r.mapcalc. Depending upon the required
accuracy, multiplying by cos(y()) might be sufficient. Otherwise, see
lib/gis/area_ellipse.c for the equations used.

--
Glynn Clements <glynn@gclements.plus.com>

Jerry Griffith wrote:
> Ok, since GRASS does account for differing cell size
> in lat/long raster files, what I am after is mapcalc is
> to do something like:
>
> croparea = [cellsize] * corn_map (my corn_map has a
> floating point value, ex. 0.29 that gives the fraction of the
> cell that is in corn cultivation)
>
> I have tried simply typing in for [cellsize] the
> functions listed below,
> but that doesn't work or my syntax is wrong.
>
> G_begin_cell_area_calculations()
> G_area_of_cell_at_row()

Glynn:

These are C functions; they don't exist in r.mapcalc.

> and I haven't found in r.stats or related commands that allows an
> intermediate multiplication step. I have tried using r.proj (and with
> -n) to reproject to an equal area projection, but my crop maps are
> global coverages and the other Locations in lambert equal area are for
> the US, so the result is my global crop maps get clipped to the US.

If you can't use a global equal-area projection, you could just use
the underlying equations within r.mapcalc. Depending upon the required
accuracy, multiplying by cos(y()) might be sufficient.
Otherwise, see lib/gis/area_ellipse.c for the equations used.

perhaps this is a nice task for the SWIG/python interface?
see swig/python/examples/rasteraccess.py

Hamish