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