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

[Moritz Lennert]

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.

> |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).

In order to avoid precision issues, I tried with both density by m2 and 
density by km2, but results are similar.

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.

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