On Thu, Dec 1, 2016 at 2:31 PM, Markus Metz <markus.metz.gisw...@gmail.com> wrote:
> > > On Thu, Dec 1, 2016 at 1:16 PM, Paulo van Breugel <p.vanbreu...@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@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 > >> > >> > > >
_______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
