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

[Markus Metz]

On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert <
mlennert at club.worldonline.be> wrote:
>
>
>
> Le 29 novembre 2016 15:10:11 GMT+01:00, Markus Metz <
markus.metz.giswork at gmail.com> a écrit :
> >On Mon, Nov 28, 2016 at 5:51 PM, Moritz Lennert <
> >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...

Indeed.

[...]
> >>
> >> - Would it be possible (desirable?) to implement some equivalent of
> >r.resamp.stats, ideally with the -w flag for weighting the sum.
> >
> >IMHO, the -w flag does not make sense because cell borders are no
> >longer
> >rectangles after reprojection, they are rather something like
> >trapezoids.
>
>
> And why does that imply that -w doesn't make sense ?

Because the -w flag in r.resamp.stats works with a rectangle overlapping
other rectangles. It only works with rectangles, i.e. different cell
extents in the same coordinate reference system (GRASS location). For
reprojection, something like r.in.xyz could work:

1. convert raster cells to vector points with raster cell value as vector z
value.
2. reproject these vector points to the target location
3. export the vector points with v.out.ascii format=point
4. import the exported points with r.in.xyz method=sum

The method options of r.in.xyz could be added to r.proj for spatial
aggregation, but that would mean that most of the code of r.in.xyz would
need to be copied to r.proj, then adapted. That would be an interesting
enhancement, also considering that gdalwarp offers as resampling methods
near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med,
Q1, Q3.

Markus M
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20161129/beedce28/attachment-0001.html>