On Wed, Nov 30, 2016 at 2:12 PM, Markus Metz <markus.metz.gisw...@gmail.com> wrote:
> > > On Wed, Nov 30, 2016 at 1:22 PM, Moritz Lennert < > mlenn...@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 > >>> <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 <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@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
