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

Wed Nov 30 00:40:20 PST 2016

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";|

>From 
https://pvanb.wordpress.com/2012/12/17/calculating-the-raster-cell-area-of-an-unprojected-raster-layer/. 
In the comment section of that post there is a link to an addon on 
github that computes the (fractional) total area of all grid cells in a 
raster map. I haven't tried it out, but perhaps this has some relevant code.

>
> 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 <http://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 <http://r.in.xyz> method=sum
>
> The method options of r.in.xyz <http://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 <http://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
>
>
>
> _______________________________________________
> 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/bb0076d4/attachment.html>