Opened 5 years ago

Closed 5 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: grass-dev@…
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)

r_object_geometry_coords.diff (1.9 KB ) - added by mlennert 5 years ago.

Download all attachments as: .zip

Change History (5)

by mlennert, 5 years ago

comment:1 by marisn, 5 years ago

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 row

And the code itself:

x = Rast_col_to_easting(0.5, &current_region2);
y = Rast_row_to_northing(current_row + 0.5, &current_region2);

in reply to:  1 comment:2 by mlennert, 5 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 row

And the code itself:

x = Rast_col_to_easting(0.5, &current_region2);
y = Rast_row_to_northing(current_row + 0.5, &current_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, &current_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 ?

in reply to:  description ; comment:3 by mmetz, 5 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=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.

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.

in reply to:  3 comment:4 by mlennert, 5 years ago

Resolution: fixed
Status: newclosed

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.

Note: See TracTickets for help on using tickets.