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

[Paulo van Breugel]

On Thu, Dec 1, 2016 at 2:31 PM, Markus Metz <markus.metz.giswork at gmail.com>
wrote:

>
>
> On Thu, Dec 1, 2016 at 1:16 PM, Paulo van Breugel <p.vanbreugel at gmail.com>
> wrote:
> >
> >
> >
> > 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 <0415%20926%20536>|
> >>> |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).
>
> 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.
> >
> >
> > 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.
>
> Why not per hectare? The target cell size is close to one hectare, it
> could be set to exactly one hectare.
> >
> > Does it make a difference if you use a much higher target resolution?
>
> In theory, a higher target resolution should result in a higher total sum,
> a lower target resolution in a lower total sum.
> >>
> >>
> >> 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.
>
> Try the new area() function of r.mapcalc ;-)
>

Missed that in the meantime it was already implemented, awesome!


>
> in the source location
> record absolute count as src_count
> convert absolute count to density per hectare with
> r.mapcalc "pop_dens_ha = 10000.0 * pop_count / area()"
>
> in the target location
> set the region according to r.proj -g
> use reasonable grid geometry with e.g. g.region -pa res=100
> reproject pop_dens_ha with r.proj method=bilinear_f
> now density is equal to the absolute count because density is per hectare
> and cell size is 10000 square meter
> record absolute count as tgt_count
>
> fix any deviations between src_count and tgt_count with
> r.mapcalc "pop_count_calibrated = pop_dens_ha * $src_count / $tgt_count"
>
> this calbration could also be done for administrative units separately
>
> Markus M
> >
> >
> >>
> >> 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
> >>
> >>
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20161201/ad8e8b97/attachment.html>