[GRASS-dev] how to reproject a raster map with absolute numbers without losing data

[Paulo van Breugel]

On 01-12-16 12:01, Moritz Lennert wrote:
> On 30/11/16 09:40, Paulo van Breugel wrote:
>>
>>> On Tue, Nov 29, 2016 at 5:44 PM, Moritz Lennert
>
>>>
>>> > 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 at PERMANENT|
>> |PI=3.1415926536|
>> |export `g.region -g`|
>> |export `g.proj -j`|
>
> This doesn't work for me as the output is
>
> g.proj -j
>
> +proj=longlat
> +no_defs
> +a=6378137
> +rf=298.257223563
> +towgs84=0.000,0.000,0.000
>
> and +a is not acceptable as variable name.

You are right, that should  be g.proj -g I think

>
>> |r.mapcalc "rcs = (sin(y()+$nsres/2) - sin(y()-$nsres/2)) * ($ewres *
>> $PI/180) * float($a)^2";|
>>
>
> I tested this approach as well. I.e.
>
> area = above formula
> pop_density = pop / area
>
> reproject pop_density
>
> pop = pop_density * (ewres*nsres)
>
> It is faster than the r.in.xyz approach. But it does not seem to be as 
> precise.
>
> I still lost about 100.000 inhabitants, and more when I smooth more 
> (the "loss" increases from nearest neighbor to bilinear to bicubic).

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.

>
> In order to avoid precision issues, I tried with both density by m2 
> and density by km2, but results are similar.
Does it make a difference if you use a much higher target resolution?
>
> 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.

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.

>
> But I do think that if it is important to maintain exact totals, then 
> the r.to.vect - v.proj - v.out.ascii - r.in.xyz solution seems to be 
> more reliable.
>
> Moritz
>
>