On 29-11-16 22:33, Markus Metz wrote:


On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert <mlenn...@club.worldonline.be <mailto:mlenn...@club.worldonline.be>> wrote:
>
>
>
> Le 29 novembre 2016 15:10:11 GMT+01:00, Markus Metz <markus.metz.gisw...@gmail.com <mailto:markus.metz.gisw...@gmail.com>> a écrit :
> >On Mon, Nov 28, 2016 at 5:51 PM, Moritz Lennert <
> >mlenn...@club.worldonline.be <mailto:mlenn...@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@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@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

Reply via email to