[GRASS-dev] how to reproject a raster map with absolute numbers without losing data

Wed Nov 30 08:19:21 PST 2016

On Wed, Nov 30, 2016 at 2:12 PM, Markus Metz <markus.metz.giswork at gmail.com>
wrote:

>
>
> On Wed, Nov 30, 2016 at 1:22 PM, Moritz Lennert <
> mlennert at club.worldonline.be> wrote:
> >
> > On 30/11/16 09:40, Paulo van Breugel wrote:
> >>
> >>
> >>
> >> On 29-11-16 22:33, Markus Metz wrote:
> >>>
> >>>
> >>>
> >>> On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert
> >>> <mlennert at club.worldonline.be <mailto:mlennert at club.worldonline.be>>
> >>> wrote:
> >>> >
> >>> >
> >>> >
> >>> > Le 29 novembre 2016 15:10:11 GMT+01:00, Markus Metz
> >>> <markus.metz.giswork at gmail.com <mailto:markus.metz.giswork at gmail.com>>
> >>> a écrit :
> >>> > >On Mon, Nov 28, 2016 at 5:51 PM, Moritz Lennert <
> >>> > >mlennert at club.worldonline.be <mailto:mlennert at club.worldonline.be>>
> >>>
> >>> wrote:
> >>> > >>
> >>> > >> Hi,
> >>> > >>
> >>> > >> In discussions with students the following problem has arisen, and
> >>> > >I'm
> >>> > >not sure how to respond to this:
> >>> > >>
> >>> > >> When reprojecting raster data, one has to chose a method of
> >>> > >interpolation. However, when the data you have represents absolute
> >>> > >numbers
> >>> > >(such as population grids), you do not wish to interpolate. What is
> >>> > >needed
> >>> > >is a technique to redistribute all units (the counts) into a new
> >>> > >partition
> >>> > >of space. This is implement in r.resamp.stats, but is not possible
> with
> >>> > >r.proj.
> >>> > >>
> >>> > >> As an example, I took the 2010 Worldpop population grid [1] of
> >>> > >Madagascar. The data is described as:
> >>> > >>
> >>> > >> SPATIAL RESOLUTION: 0.000833333 decimal degrees (approx 100m at
> the
> >>> > >equator)
> >>> > >> PROJECTION: Geographic, WGS84
> >>> > >>
> >>> > >> I import the data into a epsg:4326 location and calculate the sum
> of
> >>> > >the
> >>> > >pixel values, i.e. the total population:
> >>> > >>
> >>> > >> r.univar worldpop_madagascar_2010
> >>> > >>
> >>> > >> total null and non-null cells: 143500959
> >>> > >> total null cells: 70052410
> >>> > >>
> >>> > >> Of the non-null cells:
> >>> > >> ----------------------
> >>> > >> n: 73448549
> >>> > >> minimum: 0
> >>> > >> maximum: 1163.43
> >>> > >> range: 1163.43
> >>> > >> mean: 0.282021
> >>> > >> mean of absolute values: 0.282021
> >>> > >> standard deviation: 4.3961
> >>> > >> variance: 19.3257
> >>> > >> variation coefficient: 1558.79 %
> >>> > >> sum: 20714000.1804286
> >>> > >>
> >>> > >> I then reproject the data to UTM 38S (epsg:32738), by first
> >>> > >estimating
> >>> > >the correct region settings with
> >>> > >>
> >>> > >> r.proj -g location=LL_WGS84 mapset=mlennert
> >>> > >input=worldpop_madagascar_2010
> >>> > >>
> >>> > >> which gives me
> >>> > >>
> >>> > >> n=8678821.71318899 s=7156445.80116155 w=302657.5841309
> >>> > >e=1098038.55818598
> >>> > >rows=16387 cols=8757
> >>> > >>
> >>> > >> I then set the region accordingly in order to stay close to the
> >>> > >original
> >>> > >grid and reproject
> >>> > >>
> >>> > >> g.region n=8678821.71318899 s=7156445.80116155 w=302657.5841309
> >>> > >e=1098038.55818598 rows=16387 cols=8757
> >>> > >> r.proj location=LL_WGS84 mapset=mlennert
> >>> > >input=worldpop_madagascar_2010
> >>> > >>
> >>> > >> When I then run r.univar on the resulting map, I get:
> >>> > >>
> >>> > >> r.univar worldpop_madagascar_2010 100%
> >>> > >> total null and non-null cells: 143500959
> >>> > >> total null cells: 73263298
> >>> > >>
> >>> > >> Of the non-null cells:
> >>> > >> ----------------------
> >>> > >> n: 70237661
> >>> > >> minimum: 0
> >>> > >> maximum: 1163.43
> >>> > >> range: 1163.43
> >>> > >> mean: 0.281982
> >>> > >> mean of absolute values: 0.281982
> >>> > >> standard deviation: 4.37828
> >>> > >> variance: 19.1693
> >>> > >> variation coefficient: 1552.68 %
> >>> > >> sum: 19805754.8529859
> >>> > >>
> >>> > >> Madagascar has just lost 908,246 inhabitants !
> >>> > >
> >>> > >You would need to convert absolute numbers to densities before
> >>> > >reprojection. Note that in latlong, cells have different sizes when
> >>> > >expressed in square meters.
> >>> >
> >>> > Yes. Do we have any existing module to calculate area by pixel ?
> >>>
> >>> Not that I know of. Only for vector areas.
> >>>
> >>> > An area() function in r.mapcalc would be nice...
> >>
> >>
> >> Would indeed be a nice addition, but it isn't too difficult to compute,
> >> using:
> >>
> >> |g.region rast=afripop10 at PERMANENT|
> >> |PI=3.1415926536 <0415%20926%20536>|
> >> |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";|
> >>
> >
> >
> > Yes, obviously, but:
> >
> > - Not many people will remember such formulas by heart, and so having a
> function implementing it would be a nice convenience.
> >
> > - This function does the calculations on the sphere. For most
> applications this is precise enough, but there might be some cases where
> taking into account the actual ellipsoid might be better. I don't know how
> far the geodesic distance functions in GRASS go concerning this.
>

True, good enough for Afripop given the large error margin in that data.
However, would be great if something better / more accurate could be
implemented in r.mapcalc.

>
> GRASS has G_area_of_cell_at_row() in lib/gis/area.c. A new area() function
> for r.mapcalc would only need to initialize calculations with
> G_begin_cell_area_calculations() which considers the actual ellipsoid,
> then call G_area_of_cell_at_row().
>

I saw G_area_of_cell_at_row() was used in this addon on github (
https://github.com/joergsteinkamp/r.area), but that is in C and I was
wondering if / how this can be used / accessed in a python script.

> Markus M
>
> >
> > In the case of Afripop/Worldpop the error of the estimation is probably
> much bigger than the error of approximating the earth to a sphere, so going
> beyond would probably be overkill.
> >
> > Moritz
>
>
> _______________________________________________
> grass-dev mailing list
> grass-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20161130/cc9eef40/attachment-0001.html>