Opened 6 years ago
Closed 6 years ago
#3680 closed enhancement (fixed)
[PATCH] r.object.geometry: add output of mean coordinates of objects - difference to r.mapcalc/r.univar output
Reported by: | mlennert | Owned by: | |
---|---|---|---|
Priority: | normal | Milestone: | |
Component: | Addons | Version: | unspecified |
Keywords: | r.object.geometry cell coordinates | Cc: | |
CPU: | Unspecified | Platform: | Unspecified |
Description
For some applications it can be useful to obtain the mean coordinates of raster objects. I would like to add this functionality to r.object.geometry. For this I use the attached patch.
However, when I use this to calculate mean coordinates of objects, I get one unit of resolution difference between these results and the results of the following r.mapcalc/r.univar call:
g.region rast=MyMap r.mapcalc "x = x()" r.univar -t x zones=MyMap
Here's some output for comparison:
r.mapcalc/r.univar r.object.geometry zone|mean cat|mean_x 1|270347.103514056 1|270347.353514 2|270355.499712833 2|270355.749713 3|270243.128030303 3|270243.378030 4|270270.053380443 4|270270.303380 5|270280.90239726 5|270281.152397 6|270277.778301887 6|270278.028302
The result of r.object.geometry is always shifted one resolution (0.25) to the East.
However, y values are not:
r.mapcalc/r.univar r.object.geometry zone|mean cat|mean_x 1|154056.164959839 1|154056.164960 2|153996.903084224 2|153996.903084 3|154041.199242424 3|154041.199242 4|154016.742456073 4|154016.742456 5|154017.476598174 5|154017.476598 6|154017.057193396 6|154017.057193
Is the error in my patch, or maybe in the r.mapcalc x() function ?
Any objections to applying this patch ? It changes the output of the module, but it add the additional values at the end, so unless some scripts rely on fd being the last value, this should not cause too much trouble.
Moritz
Attachments (1)
Change History (5)
by , 6 years ago
Attachment: | r_object_geometry_coords.diff added |
---|
follow-up: 2 comment:1 by , 6 years ago
comment:2 by , 6 years ago
Replying to marisn:
Comments in the code of r.mapcalc (https://trac.osgeo.org/grass/browser/grass/trunk/raster/r.mapcalc/xcoor.c#L27):
x() easting at center of column y() northing at center of rowAnd the code itself:
x = Rast_col_to_easting(0.5, ¤t_region2); y = Rast_row_to_northing(current_row + 0.5, ¤t_region2);
Yes, because in the programmer's manual, it says:
/*! 250 * \brief Column to easting. 251 * 252 * Converts a <i>col</i> relative to a <i>window</i> to an easting. 253 * 254 * <b>Note:</b> <i>col</i> is a <i>double</i>: 255 * - col+0.0 will return the easting for the western edge of the column. 256 * - col+0.5 will return the easting for the center of the column. 257 * - col+1.0 will return the easting for the eastern edge of the column. 258 * 259 * \param col column number 260 * \param[in] window pointer to Cell_head 261 * 262 * \return east coordinate 263 */ 264 double Rast_col_to_easting(double col, const struct Cell_head *window) 265 { 266 return window->west + col * window->ew_res; 267 }
And in the code of r.mapcalc, you can see
27 x = Rast_col_to_easting(0.5, ¤t_region2); 28 29 for (i = 0; i < columns; i++) { 30 res[i] = x; 31 x += current_region2.ew_res; 32 }
i.e. the first x is calculated for col = 0.5 and all others for respectively one ew_res to the east, so always col + 0.5.
That is why I'm surprised to not get the same answer.
Maybe it is due more to r.univar than to r.mapcalc ?
follow-up: 4 comment:3 by , 6 years ago
Replying to mlennert:
For some applications it can be useful to obtain the mean coordinates of raster objects. I would like to add this functionality to r.object.geometry. For this I use the attached patch.
However, when I use this to calculate mean coordinates of objects, I get one unit of resolution difference between these results and the results of the following r.mapcalc/r.univar call:
g.region rast=MyMap r.mapcalc "x = x()" r.univar -t x zones=MyMapHere's some output for comparison:
r.mapcalc/r.univar r.object.geometry zone|mean cat|mean_x 1|270347.103514056 1|270347.353514 2|270355.499712833 2|270355.749713 3|270243.128030303 3|270243.378030 4|270270.053380443 4|270270.303380 5|270280.90239726 5|270281.152397 6|270277.778301887 6|270278.028302The result of r.object.geometry is always shifted one resolution (0.25) to the East.
the reason is that r.object.geometry pads the current region with one column to the west and east, and the module iterates over columns with
for (col = 1; col <= ncols; col++) { ...
that also means that 0.5 needs to be subtracted:
obj_geos[cur - min].x += Rast_col_to_easting(col - 0.5, &cellhd);
to get the real center for each column in the current region.
comment:4 by , 6 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
Replying to mmetz:
Replying to mlennert:
The result of r.object.geometry is always shifted one resolution (0.25) to the East.
the reason is that r.object.geometry pads the current region with one column to the west
Right, yes, I should have remembered that (or identified it in the code). Sorry.
Thanks !
I committed the change in r73523.
Comments in the code of r.mapcalc (https://trac.osgeo.org/grass/browser/grass/trunk/raster/r.mapcalc/xcoor.c#L27):
And the code itself: