On 02-12-16 11:19, Moritz Lennert wrote:
On 01/12/16 14:31, Markus Metz wrote:
On Thu, Dec 1, 2016 at 1:16 PM, Paulo van Breugel
<p.vanbreu...@gmail.com <mailto: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
No this does not provide a, but the name of the ellipsoid. I'm not
sure that you can currently directly parse g.proj output for a...
It does to me:
GRASS 7.3.svn (latlon):~ > g.proj -g
name=Lat/Lon
proj=ll
datum=wgs84
a=6378137
es=0.006694379990141316
no_defs=defined
unit=degree
units=degrees
meters=1.0
GRASS 7.3.svn (latlon):~ >
This works for me:
eval $(g.proj -j | awk '{if(match($1,"=")>0){sub("+", "", $1); print
$1}}')
|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 <http://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.
Why is that if we work in densities ?
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 ;-)
Thanks !!
FYI, for all of Madagascar, for pixels of resolution 0:00:02.99988
(i.e. +/- 100m at the equator), I get an RMSE of about 42m2 between
Paulo's formula and area() in r.mapcalc. The error varies between 31
and 52m2, with highest values in the North, i.e. the closer you get to
the equator.
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
Actually, using area(), and a 100m target resolution, I get closer to
the source total with bilinear than with bilinear_f, but in all cases
the target total is now _above_ the source total:
Source total: 20714000
Target total - Source total:
nn: 22473
bilinear: 2300
bilinear_f: 26873
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
Yes, obviously, by calibrating this way we will get the same total by
definition. But this means that we postulate that the "error" is
distributed equally across space...
Thanks to both of you for the very helpful discussion. I'm afraid few
people dealing with these data think much about this problem.
Moritz
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
