<div dir="ltr"><div><br><br>On Thu, Dec 1, 2016 at 1:16 PM, Paulo van Breugel <<a href="mailto:p.vanbreugel@gmail.com">p.vanbreugel@gmail.com</a>> wrote:<br>><br>><br>><br>> On 01-12-16 12:01, Moritz Lennert wrote:<br>>><br>>> On 30/11/16 09:40, Paulo van Breugel wrote:<br>>>><br>>>><br>>>>> On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert<br>>><br>>><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=3.1415926536|<br>>>> |export `g.region -g`|<br>>>> |export `g.proj -j`|<br>>><br>>><br>>> This doesn't work for me as the output is<br>>><br>>> g.proj -j<br>>><br>>> +proj=longlat<br>>> +no_defs<br>>> +a=6378137<br>>> +rf=298.257223563<br>>> +towgs84=0.000,0.000,0.000<br>>><br>>> and +a is not acceptable as variable name.<br>><br>><br>> You are right, that should  be g.proj -g I think<br>><br>>><br>>>> |r.mapcalc "rcs = (sin(y()+$nsres/2) - sin(y()-$nsres/2)) * ($ewres *<br>>>> $PI/180) * float($a)^2";|<br>>>><br>>><br>>> I tested this approach as well. I.e.<br>>><br>>> area = above formula<br>>> pop_density = pop / area<br>>><br>>> reproject pop_density<br>>><br>>> pop = pop_density * (ewres*nsres)<br>>><br>>> It is faster than the <a href="http://r.in.xyz">r.in.xyz</a> approach. But it does not seem to be as precise.<br>>><br>>> I still lost about 100.000 inhabitants, and more when I smooth more (the "loss" increases from nearest neighbor to bilinear to bicubic).<br><br></div>Note that some of the loss happens when a non-zero cell is bordering a NULL cell, unless you use one of the *_f interpolation methods.<br><div>><br>><br>> Just guessing, but an issue I had in the past with layers like Afripop is that it may contain cells with extremely high values (a result of the inter/extrapolation methods used to create those layers I guess, not sure it is an issue for e.g., Wordpop). Depending on the resampling method, this could also add to the differences in totals.<br>><br>>><br>>> In order to avoid precision issues, I tried with both density by m2 and density by km2, but results are similar.<br><br></div><div>Why not per hectare? The target cell size is close to one hectare, it could be set to exactly one hectare.<br></div><div>><br>> Does it make a difference if you use a much higher target resolution?<br><br></div><div>In theory, a higher target resolution should result in a higher total sum, a lower target resolution in a lower total sum.<br></div><div>>><br>>><br>>> I don't know which part of the error comes from the use of spheroid approximation and which part comes from the resampling at reprojection.<br>><br>><br>> It would be interesting to see how results are if using the same routine as above (convert to densities, reproject, convert back to totals), but using the more precise G_area_of_cell_at_row() function to compute the cell areas.<br><br></div><div>Try the new area() function of r.mapcalc ;-)<br><br></div><div>in the source location<br></div><div>record absolute count as src_count<br></div><div>convert absolute count to density per hectare with<br></div><div>r.mapcalc "pop_dens_ha = 10000.0 * pop_count / area()"<br><br></div><div>in the target location<br></div><div>set the region according to r.proj -g<br></div><div>use reasonable grid geometry with e.g. g.region -pa res=100<br></div><div>reproject pop_dens_ha with r.proj method=bilinear_f<br></div><div>now density is equal to the absolute count because density is per hectare and cell size is 10000 square meter<br></div><div><div>record absolute count as tgt_count<br></div></div><div><br></div><div>fix any deviations between src_count and tgt_count with<br></div><div>r.mapcalc "pop_count_calibrated = pop_dens_ha * $src_count / $tgt_count"<br></div><div><br></div><div>this calbration could also be done for administrative units separately<br></div><div><br></div><div>Markus M<br></div><div>><br>><br>>><br>>> But I do think that if it is important to maintain exact totals, then the r.to.vect - v.proj - v.out.ascii - <a href="http://r.in.xyz">r.in.xyz</a> solution seems to be more reliable.<br>>><br>>> Moritz<br>>><br>>><br>><br></div></div>