<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 30, 2016 at 2:12 PM, Markus Metz <span dir="ltr"><<a href="mailto:markus.metz.giswork@gmail.com" target="_blank">markus.metz.giswork@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><br><br>On Wed, Nov 30, 2016 at 1:22 PM, Moritz Lennert <<a href="mailto:mlennert@club.worldonline.be" target="_blank">mlennert@club.worldonline.be</a>> wrote:<br>><br>> On 30/11/16 09:40, Paulo van Breugel wrote:<br>>><br>>><br>>><br>>> On 29-11-16 22:33, Markus Metz wrote:<br>>>><br>>>><br>>>><br>>>> On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert<br>>>> <<a href="mailto:mlennert@club.worldonline.be" target="_blank">mlennert@club.worldonline.be</a> <mailto:<a href="mailto:mlennert@club.worldonline.be" target="_blank">mlennert@club.<wbr>worldonline.be</a>>><br>>>> wrote:<br>>>> ><br>>>> ><br>>>> ><br>>>> > Le 29 novembre 2016 15:10:11 GMT+01:00, Markus Metz<br>>>> <<a href="mailto:markus.metz.giswork@gmail.com" target="_blank">markus.metz.giswork@gmail.com</a> <mailto:<a href="mailto:markus.metz.giswork@gmail.com" target="_blank">markus.metz.giswork@<wbr>gmail.com</a>>><br>>>> a écrit :<br>>>> > >On Mon, Nov 28, 2016 at 5:51 PM, Moritz Lennert <<br>>>> > ><a href="mailto:mlennert@club.worldonline.be" target="_blank">mlennert@club.worldonline.be</a> <mailto:<a href="mailto:mlennert@club.worldonline.be" target="_blank">mlennert@club.<wbr>worldonline.be</a>>><br>>>><br>>>> wrote:<br>>>> > >><br>>>> > >> Hi,<br>>>> > >><br>>>> > >> In discussions with students the following problem has arisen, and<br>>>> > >I'm<br>>>> > >not sure how to respond to this:<br>>>> > >><br>>>> > >> When reprojecting raster data, one has to chose a method of<br>>>> > >interpolation. However, when the data you have represents absolute<br>>>> > >numbers<br>>>> > >(such as population grids), you do not wish to interpolate. What is<br>>>> > >needed<br>>>> > >is a technique to redistribute all units (the counts) into a new<br>>>> > >partition<br>>>> > >of space. This is implement in r.resamp.stats, but is not possible with<br>>>> > >r.proj.<br>>>> > >><br>>>> > >> As an example, I took the 2010 Worldpop population grid [1] of<br>>>> > >Madagascar. The data is described as:<br>>>> > >><br>>>> > >> SPATIAL RESOLUTION: 0.000833333 decimal degrees (approx 100m at the<br>>>> > >equator)<br>>>> > >> PROJECTION: Geographic, WGS84<br>>>> > >><br>>>> > >> I import the data into a epsg:4326 location and calculate the sum of<br>>>> > >the<br>>>> > >pixel values, i.e. the total population:<br>>>> > >><br>>>> > >> r.univar worldpop_madagascar_2010<br>>>> > >><br>>>> > >> total null and non-null cells: 143500959<br>>>> > >> total null cells: 70052410<br>>>> > >><br>>>> > >> Of the non-null cells:<br>>>> > >> ----------------------<br>>>> > >> n: 73448549<br>>>> > >> minimum: 0<br>>>> > >> maximum: 1163.43<br>>>> > >> range: 1163.43<br>>>> > >> mean: 0.282021<br>>>> > >> mean of absolute values: 0.282021<br>>>> > >> standard deviation: 4.3961<br>>>> > >> variance: 19.3257<br>>>> > >> variation coefficient: 1558.79 %<br>>>> > >> sum: 20714000.1804286<br>>>> > >><br>>>> > >> I then reproject the data to UTM 38S (epsg:32738), by first<br>>>> > >estimating<br>>>> > >the correct region settings with<br>>>> > >><br>>>> > >> r.proj -g location=LL_WGS84 mapset=mlennert<br>>>> > >input=worldpop_madagascar_<wbr>2010<br>>>> > >><br>>>> > >> which gives me<br>>>> > >><br>>>> > >> n=8678821.71318899 s=7156445.80116155 w=302657.5841309<br>>>> > >e=1098038.55818598<br>>>> > >rows=16387 cols=8757<br>>>> > >><br>>>> > >> I then set the region accordingly in order to stay close to the<br>>>> > >original<br>>>> > >grid and reproject<br>>>> > >><br>>>> > >> g.region n=8678821.71318899 s=7156445.80116155 w=302657.5841309<br>>>> > >e=1098038.55818598 rows=16387 cols=8757<br>>>> > >> r.proj location=LL_WGS84 mapset=mlennert<br>>>> > >input=worldpop_madagascar_<wbr>2010<br>>>> > >><br>>>> > >> When I then run r.univar on the resulting map, I get:<br>>>> > >><br>>>> > >> r.univar worldpop_madagascar_2010 100%<br>>>> > >> total null and non-null cells: 143500959<br>>>> > >> total null cells: 73263298<br>>>> > >><br>>>> > >> Of the non-null cells:<br>>>> > >> ----------------------<br>>>> > >> n: 70237661<br>>>> > >> minimum: 0<br>>>> > >> maximum: 1163.43<br>>>> > >> range: 1163.43<br>>>> > >> mean: 0.281982<br>>>> > >> mean of absolute values: 0.281982<br>>>> > >> standard deviation: 4.37828<br>>>> > >> variance: 19.1693<br>>>> > >> variation coefficient: 1552.68 %<br>>>> > >> sum: 19805754.8529859<br>>>> > >><br>>>> > >> Madagascar has just lost 908,246 inhabitants !<br>>>> > ><br>>>> > >You would need to convert absolute numbers to densities before<br>>>> > >reprojection. Note that in latlong, cells have different sizes when<br>>>> > >expressed in square meters.<br>>>> ><br>>>> > Yes. Do we have any existing module to calculate area by pixel ?<br>>>><br>>>> Not that I know of. Only for vector areas.<br>>>><br>>>> > An area() function in r.mapcalc would be nice...<br>>><br>>><br>>> Would indeed be a nice addition, but it isn't too difficult to compute,<br>>> using:<br>>><br>>> |g.region rast=afripop10@PERMANENT|<br>>> |PI=<a href="tel:0415%20926%20536" value="+31415926536" target="_blank">3.1415926536</a>|<br>>> |export `g.region -g`|<br>>> |export `g.proj -j`|<br>>> |r.mapcalc "rcs = (sin(y()+$nsres/2) - sin(y()-$nsres/2)) * ($ewres *<br>>> $PI/180) * float($a)^2";|<br>>><br>><br>><br>> Yes, obviously, but:<br>><br>> - Not many people will remember such formulas by heart, and so having a function implementing it would be a nice convenience.<br>><br>> - 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.<br></div></div></div></blockquote><div><br></div><div>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.<br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><br></div>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_<wbr>calculations() which considers the actual ellipsoid, then call G_area_of_cell_at_row().<br></div></div></blockquote><div><br></div><div>I saw G_area_of_cell_at_row() was used in this addon on github (<a href="https://github.com/joergsteinkamp/r.area">https://github.com/joergsteinkamp/r.area</a>), but that is in C and I was wondering if / how this can be used / accessed in a python script. <br><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><br></div>Markus M<br><div><div><br>><br>> 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.<br>><br>> Moritz<br><br></div></div></div>
<br>______________________________<wbr>_________________<br>
grass-dev mailing list<br>
<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/grass-dev" rel="noreferrer" target="_blank">http://lists.osgeo.org/<wbr>mailman/listinfo/grass-dev</a><br></blockquote></div><br></div></div>