<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body smarttemplateinserted="true" bgcolor="#FFFFFF" text="#000000">
<p><br>
</p>
<br>
<div class="moz-cite-prefix">On 29-11-16 22:33, Markus Metz wrote:<br>
</div>
<blockquote
cite="mid:CAG+h=FH2+HGzPjha1ud7A01NRr=0sdR9zBZ210Yeu=1kqiksiQ@mail.gmail.com"
type="cite">
<div dir="ltr">
<div><br>
<br>
On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert <<a
moz-do-not-send="true"
href="mailto:mlennert@club.worldonline.be">mlennert@club.worldonline.be</a>>
wrote:<br>
><br>
><br>
><br>
> Le 29 novembre 2016 15:10:11 GMT+01:00, Markus Metz <<a
moz-do-not-send="true"
href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>>
a écrit :<br>
> >On Mon, Nov 28, 2016 at 5:51 PM, Moritz Lennert <<br>
> ><a moz-do-not-send="true"
href="mailto:mlennert@club.worldonline.be">mlennert@club.worldonline.be</a>>
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_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_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>
</div>
Not that I know of. Only for vector areas.<br>
<div><br>
> An area() function in r.mapcalc would be nice...<br>
</div>
</div>
</blockquote>
<br>
Would indeed be a nice addition, but it isn't too difficult to
compute, using:<br>
<br>
<div class="container">
<div class="line number1 index0 alt2"><code class="plain plain">g.region
rast=afripop10@PERMANENT</code></div>
<div class="line number2 index1 alt1"><code class="plain plain">PI=3.1415926536</code></div>
<div class="line number3 index2 alt2"><code class="plain plain">export
`g.region -g`</code></div>
<div class="line number4 index3 alt1"><code class="plain plain">export
`g.proj -j`</code></div>
<div class="line number5 index4 alt2"><code class="plain plain">r.mapcalc
"rcs = (sin(y()+$nsres/2) - sin(y()-$nsres/2)) * ($ewres *
$PI/180) * float($a)^2";</code></div>
</div>
<br>
From
<a class="moz-txt-link-freetext" href="https://pvanb.wordpress.com/2012/12/17/calculating-the-raster-cell-area-of-an-unprojected-raster-layer/">https://pvanb.wordpress.com/2012/12/17/calculating-the-raster-cell-area-of-an-unprojected-raster-layer/</a>.
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.<br>
<br>
<blockquote
cite="mid:CAG+h=FH2+HGzPjha1ud7A01NRr=0sdR9zBZ210Yeu=1kqiksiQ@mail.gmail.com"
type="cite">
<div dir="ltr">
<div><br>
Indeed.<br>
</div>
<div><br>
[...]<br>
> >><br>
> >> - Would it be possible (desirable?) to implement
some equivalent of<br>
> >r.resamp.stats, ideally with the -w flag for
weighting the sum.<br>
> ><br>
> >IMHO, the -w flag does not make sense because cell
borders are no<br>
> >longer<br>
> >rectangles after reprojection, they are rather
something like<br>
> >trapezoids.<br>
><br>
><br>
> And why does that imply that -w doesn't make sense ?<br>
<br>
</div>
<div>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 <a moz-do-not-send="true" href="http://r.in.xyz">r.in.xyz</a>
could work:<br>
<br>
</div>
<div>1. convert raster cells to vector points with raster cell
value as vector z value.<br>
</div>
<div>2. reproject these vector points to the target location<br>
</div>
<div>3. export the vector points with v.out.ascii format=point<br>
</div>
<div>4. import the exported points with <a
moz-do-not-send="true" href="http://r.in.xyz">r.in.xyz</a>
method=sum<br>
<br>
</div>
<div>The method options of <a moz-do-not-send="true"
href="http://r.in.xyz">r.in.xyz</a> could be added to r.proj
for spatial aggregation, but that would mean that most of the
code of <a moz-do-not-send="true" href="http://r.in.xyz">r.in.xyz</a>
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.<br>
</div>
<div><br>
</div>
<div>Markus M<br>
</div>
<div><br>
</div>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<br>
<pre wrap="">_______________________________________________
grass-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="http://lists.osgeo.org/mailman/listinfo/grass-dev">http://lists.osgeo.org/mailman/listinfo/grass-dev</a></pre>
</blockquote>
<br>
</body>
</html>