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

[Markus Metz]

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

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().

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20161130/3c9dbf7e/attachment-0001.html>