Re: [GRASS-dev] error v.what.rast
On Fri, Oct 10, 2014 at 1:44 PM, Glynn Clements wrote: [...] Thanks for the detailed analysis, Glynn. > 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). To avoid that it gets lost here, I have submitted that in r62254 (trunk for now, please test). Markus ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Fri, Oct 10, 2014 at 1:44 PM, Glynn Clements wrote: > > 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. > Yes, the way grass handles points outside the region makes sense to me; it follows the general philosophy of grass to ignore what is outside the region. > > 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. > What would make sense to me is if for those points to update the table with a NULL value, as that is what it really is (considering that all values within the defined region are valid). > 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 > ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
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 ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Fri, Oct 10, 2014 at 9:00 AM, Paulo van Breugel
wrote:
>
>
> On Thu, Oct 9, 2014 at 11:49 PM, Markus Neteler wrote:
>
>> On Wed, Oct 8, 2014 at 9:44 AM, Paulo van Breugel
>> wrote:
>> > On Wed, Oct 8, 2014 at 9:21 AM, Markus Neteler
>> wrote:
>> >>
>> >> On Tue, Oct 7, 2014 at 11:20 PM, Paulo van Breugel
>> >> wrote:
>> >> > I am trying to upload raster values of the raster 'grid' to the
>> >> > attribute
>> >> > table of the gbif vector layer, as below:
>> >> >
>> >> > g.region rast=grid
>> >>
>> >> can you please post
>> >>
>> >> g.region rast=grid -p
>> >
>> >
>> > GRASS 7.1.svn (latlon):~ > g.region rast=grid -p
>> > projection: 3 (Latitude-Longitude)
>> > zone: 0
>> > datum: wgs84
>> > ellipsoid: a=6378137 es=0.006694379990141316
>> > north: 37:22:30N
>> > south: 39S
>> > west: 20:12:30W
>> > east: 57:40E
>> > nsres: 0:02:30
>> > ewres: 0:02:30
>> > rows: 1833
>> > cols: 1869
>> > cells: 3425877
>>
>> Can you please post with "-g" which is decimal?
>>
>
> n=37.375
> s=-34.875
> w=-17.54167
> e=51.41667
> nsres=0.0417
> ewres=0.0417
> rows=1734
> cols=1655
> cells=2869770
>
>
>>
>> >>
>> >>
>> >> > v.what.rast map=gbif raster=grid column=cell2arcm
>> >> > Reading features from vector map...
>> >> > 100%
>> >> > ERROR: Reading raster map request for row 1830 is outside
>> region
>> >>
>> >> ... and also the boundary box of gbif?
>> >
>> > Projection: Latitude-Longitude
>> >
>> > N: 37:19:58.08N
>> > S: 38:52:30S
>> > E: 57:38:24E
>> > W: 20:02:49.92W
>>
>> Also here, can you please post with "-g" which is decimal?
>>
>
> n=37.3328
> s=-38.875
> w=-20.0472
> e=57.64
> nsres=0.0416663750683434
> ewres=0.0416776824034335
> rows=1829
> cols=1864
> cells=3409256
>
>>
>> The error message is generated in
>> lib/raster/get_row.c
>>
>> static int compute_window_row(int fd, int row, int *cellRow)
>> {
>> struct fileinfo *fcb = &R__.fileinfo[fd];
>> double f;
>> int r;
>>
>> /* check for row in window */
>> if (row < 0 || row >= R__.rd_window.rows) {
>> G_fatal_error(_("Reading raster map <%s@%s> request for row %d
>> is outside region"),
>> fcb->name, fcb->mapset, row);
>> }
>>
>> 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? 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.
>
> WARNING: 2287133 points outside current region were skipped
> Update vector attributes...
> 100%
> v.what.rast complete. 57621 records updated.
>
>
Sorry, I wasn't completely awake yet. It is about the extend of the raster
layer to be sampled being smaller than the region... and yes, that was the
case. Reducing the region to fit within the bounds of the raster layer
solved the problem. Seems obvious now I know the error I made, sorry.
Thanks a lot for your help, much appreciated!
>
>
>>
>> Markus
>>
>
>
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Thu, Oct 9, 2014 at 11:49 PM, Markus Neteler wrote:
> On Wed, Oct 8, 2014 at 9:44 AM, Paulo van Breugel
> wrote:
> > On Wed, Oct 8, 2014 at 9:21 AM, Markus Neteler
> wrote:
> >>
> >> On Tue, Oct 7, 2014 at 11:20 PM, Paulo van Breugel
> >> wrote:
> >> > I am trying to upload raster values of the raster 'grid' to the
> >> > attribute
> >> > table of the gbif vector layer, as below:
> >> >
> >> > g.region rast=grid
> >>
> >> can you please post
> >>
> >> g.region rast=grid -p
> >
> >
> > GRASS 7.1.svn (latlon):~ > g.region rast=grid -p
> > projection: 3 (Latitude-Longitude)
> > zone: 0
> > datum: wgs84
> > ellipsoid: a=6378137 es=0.006694379990141316
> > north: 37:22:30N
> > south: 39S
> > west: 20:12:30W
> > east: 57:40E
> > nsres: 0:02:30
> > ewres: 0:02:30
> > rows: 1833
> > cols: 1869
> > cells: 3425877
>
> Can you please post with "-g" which is decimal?
>
n=37.375
s=-34.875
w=-17.54167
e=51.41667
nsres=0.0417
ewres=0.0417
rows=1734
cols=1655
cells=2869770
>
> >>
> >>
> >> > v.what.rast map=gbif raster=grid column=cell2arcm
> >> > Reading features from vector map...
> >> > 100%
> >> > ERROR: Reading raster map request for row 1830 is outside
> region
> >>
> >> ... and also the boundary box of gbif?
> >
> > Projection: Latitude-Longitude
> >
> > N: 37:19:58.08N
> > S: 38:52:30S
> > E: 57:38:24E
> > W: 20:02:49.92W
>
> Also here, can you please post with "-g" which is decimal?
>
n=37.3328
s=-38.875
w=-20.0472
e=57.64
nsres=0.0416663750683434
ewres=0.0416776824034335
rows=1829
cols=1864
cells=3409256
>
> The error message is generated in
> lib/raster/get_row.c
>
> static int compute_window_row(int fd, int row, int *cellRow)
> {
> struct fileinfo *fcb = &R__.fileinfo[fd];
> double f;
> int r;
>
> /* check for row in window */
> if (row < 0 || row >= R__.rd_window.rows) {
> G_fatal_error(_("Reading raster map <%s@%s> request for row %d
> is outside region"),
> fcb->name, fcb->mapset, row);
> }
>
> 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? 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.
WARNING: 2287133 points outside current region were skipped
Update vector attributes...
100%
v.what.rast complete. 57621 records updated.
>
> Markus
>
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Wed, Oct 8, 2014 at 9:44 AM, Paulo van Breugel
wrote:
> On Wed, Oct 8, 2014 at 9:21 AM, Markus Neteler wrote:
>>
>> On Tue, Oct 7, 2014 at 11:20 PM, Paulo van Breugel
>> wrote:
>> > I am trying to upload raster values of the raster 'grid' to the
>> > attribute
>> > table of the gbif vector layer, as below:
>> >
>> > g.region rast=grid
>>
>> can you please post
>>
>> g.region rast=grid -p
>
>
> GRASS 7.1.svn (latlon):~ > g.region rast=grid -p
> projection: 3 (Latitude-Longitude)
> zone: 0
> datum: wgs84
> ellipsoid: a=6378137 es=0.006694379990141316
> north: 37:22:30N
> south: 39S
> west: 20:12:30W
> east: 57:40E
> nsres: 0:02:30
> ewres: 0:02:30
> rows: 1833
> cols: 1869
> cells: 3425877
Can you please post with "-g" which is decimal?
>>
>>
>> > v.what.rast map=gbif raster=grid column=cell2arcm
>> > Reading features from vector map...
>> > 100%
>> > ERROR: Reading raster map request for row 1830 is outside region
>>
>> ... and also the boundary box of gbif?
>
> Projection: Latitude-Longitude
>
> N: 37:19:58.08N
> S: 38:52:30S
> E: 57:38:24E
> W: 20:02:49.92W
Also here, can you please post with "-g" which is decimal?
The error message is generated in
lib/raster/get_row.c
static int compute_window_row(int fd, int row, int *cellRow)
{
struct fileinfo *fcb = &R__.fileinfo[fd];
double f;
int r;
/* check for row in window */
if (row < 0 || row >= R__.rd_window.rows) {
G_fatal_error(_("Reading raster map <%s@%s> request for row %d
is outside region"),
fcb->name, fcb->mapset, row);
}
Perhaps a vector point just hits the extents of the computational region or map?
Markus
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Wed, Oct 8, 2014 at 9:21 AM, Markus Neteler wrote: > On Tue, Oct 7, 2014 at 11:20 PM, Paulo van Breugel > wrote: > > I am trying to upload raster values of the raster 'grid' to the attribute > > table of the gbif vector layer, as below: > > > > g.region rast=grid > > can you please post > > g.region rast=grid -p > GRASS 7.1.svn (latlon):~ > g.region rast=grid -p projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: a=6378137 es=0.006694379990141316 north: 37:22:30N south: 39S west: 20:12:30W east: 57:40E nsres: 0:02:30 ewres: 0:02:30 rows: 1833 cols: 1869 cells: 3425877 > > > v.what.rast map=gbif raster=grid column=cell2arcm > > Reading features from vector map... > > 100% > > ERROR: Reading raster map request for row 1830 is outside region > > ... and also the boundary box of gbif? > Projection: Latitude-Longitude N: 37:19:58.08N S: 38:52:30S E: 57:38:24E W: 20:02:49.92W > > > Any idea what the error message is about. I have done this done before > > without problems (same location/mapset and script), and I am really not > sure > > what I am doing different now. > > The computational region comes to mind. If that holds true, we might > improve that error message. > Markus > > I should add that it runs well on some resolutions (e.g., using 10 or 5 arc-minute resolution, not on others (like the 2.5 arc-minute above) > > Running GRASS 7.1 on Ubuntu > > > > ___ > > grass-dev mailing list > > grass-dev@lists.osgeo.org > > http://lists.osgeo.org/mailman/listinfo/grass-dev > ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] error v.what.rast
On Tue, Oct 7, 2014 at 11:20 PM, Paulo van Breugel wrote: > I am trying to upload raster values of the raster 'grid' to the attribute > table of the gbif vector layer, as below: > > g.region rast=grid can you please post g.region rast=grid -p > v.what.rast map=gbif raster=grid column=cell2arcm > Reading features from vector map... > 100% > ERROR: Reading raster map request for row 1830 is outside region ... and also the boundary box of gbif? > Any idea what the error message is about. I have done this done before > without problems (same location/mapset and script), and I am really not sure > what I am doing different now. The computational region comes to mind. If that holds true, we might improve that error message. Markus > Running GRASS 7.1 on Ubuntu > > ___ > grass-dev mailing list > grass-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/grass-dev ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
