[GRASS-dev] error v.what.rast

[Glynn Clements]

Paulo van Breugel wrote:

> > Perhaps a vector point just hits the extents of the computational region
> > or map?
> 
> You are right that the extends do not match completely. But I was under the
> impression that points outside the region simply got skipped?

I believe that *ought* to happen. The question is whether the fix
should go into v.what.rast or libraster.

> For example, if I reduce the region to a small area, and run
> v.what.rast, I get the warning message "WARNING: 2287133 points
> outside current region were skipped", after which the it happily
> continues to run updating the points that are within the region
> bounds, as per below.

The behaviour of Rast_get_row() etc is that requesting a row which is
outside of the map's bounds returns a row of nulls, but requesting a
row which is outside of the current region is a fatal error.

The former is just how raster maps are conceived. Conceptually, a
raster map has infinite extent, but all non-null values are contained
within a finite region. A raster map's stored bounds provide a loose
bound on the region containing the non-null values, as well as
avoiding requiring infinite disk space for each map.

The latter indicates a programming error. Any module which accesses a
raster map first specifies (implicitly or explicitly) a finite grid of
sample locations. Subsequent Rast_get_row() etc calls return a
one-dimensional array of values obtained by sampling the map at the
locations for one row of the grid. The error occurs when the program
requests data for a row which doesn't exist in the sampling grid which
it specified.

It would be fairly straightfoward to work around the issue within
libraster, by making compute_window_row() just return 0 instead of
raising a fatal error. But this could mask more significant issues,
resulting in modules silently returning wrong results.

Personally, I feel that v.what.rast should be fixed. Specifically, the
code

	row = Rast_northing_to_row(Points->y[0], &window);
	col = Rast_easting_to_col(Points->x[0], &window);

should be followed by:

	if (col < 0 || col >= window.cols || row < 0 || row >= window.rows)
	    continue;

The Vect_point_in_box() call just above that code doesn't appear to be
sufficient in the case where the point lies either on the boundary of
the region or just inside.

This is partly because Vect_point_in_box() uses >= and <= rather than
> and < (i.e. a point on the boundary is considered inside), and
partly because Rast_northing_to_row() etc are affected by rounding
error (window.north - window.south is not necessarily equal to
window.ns_res * window.rows, similarly for the horizontal direction).

-- 
Glynn Clements <glynn at gclements.plus.com>