region boundaries & G_get_map_row

Bill Brown brown at gis.uiuc.edu
Wed Jun 28 08:00:00 EDT 1995


>
>
>if my region is: n=23 s=1 e=16 w=1 nsres=1 ewres=1, 
>should there be 21 or 22 rows in a raster file?
>
>I'm having problems with s.sample related to G_get_map_row:
>
>  row = (int) G_northing_to_row (north, &window);
>  if (G_get_map_row (fd, maprow, row) < 0)
>
>  WARNING: [elev in PERMANENT] - read request for row 22 is outside region
>  DIAG:   Region is: n=23 s=1 e=16 w=1
>          Data point is north=1 east=4

There should be 22 rows, 
STARTING WITH ROW 0 - so row #22 would be the 23rd row

The only way you should be getting row=22 is if
(window->north - north) / window->ns_res >= 22.0
which DOES happen when north = window->south - so it seems
that since points ON north are IN the region, then points ON
south are not.  Otherwise, for two contiguous regions a point
on south of the north region would be in both regions.

>
>The data point is at the southern boundary of the region. By my 
>definition, the data point (north,east) *is* contained within
>the region. Is my definition wrong?

That's wrong.
(north, west) *is* in the region - the other corners belong to
other regions.

>
>In other words:
>In the program, I check to see if a site is contained within
>the region using:
>    if (sites[i]->east >= window.west && 
>        sites[i]->east <= window.east &&
>        sites[i]->north <= window.north && 
>        sites[i]->north >= window.south)
>Should this be changed from "<=" and ">=" to "<" and ">", respectively?
>Or should I just change the comparisons for west and south?

I think this would be right:
if (sites[i]->east >= window.west && 
         sites[i]->east < window.east &&
         sites[i]->north <= window.north && 
         sites[i]->north > window.south)

- Bill







More information about the grass-dev mailing list