g.region option to expand the computational region of about "some" pixels?

[nikos]

Let's accept that we have a set of point (vector) maps and we need to update one of their attributes by querying a raster map. We would use v.what.rast of course to get the job done.

Now, in order to save time and system resources, we would like to match the active region over the point map since we care about this area of the raster map (from which v.what.rast will grab the values and feed some column in the point map).

We could execute for example: "g.region res=ResOfRaster vect=PointMap -pa"

Matching the region over a point map is the problem! The points at the borders (happens quite often to me when the resolution of the raster map is low and so does the resolution of the region successively) are not taken into account.

To overcome the problem I had to find a way to "grow the region" so the region matching is "v.what.rast safe" :-). OK, in my case it was easy as the points where actually centroids extracted from a vector grid map. So matching the region to the original vector grid map did it.

But what if it's not the case to have a vector, enough larger than the point map, to perform a "v.what.rast safe" region matching?

(a) anybody else facing this problem?

Herewith I would like to suggest the implementation of a new option for g.region, a "v.what.rast-safe" option :-p

(b) is it difficult to build in a flag (let's say -x which will stand for eXpand) to "grow" the (2D) region in both, vertical and horizontal manner of about 1(or more?) pixel(s) (where pixel size here is meant to be the resolution of the 2D raster map to be queried)?

I imagine something that would use r.mapcalc (or r.mask), r.grow and g.region itself?

Thanks for reading, Nikos

[hamish]

did you try the g.region align=rastermap option?

fwiw, -a grows outwards.

Matching the region over a point map is the problem! The points at the borders (happens quite often to me when the resolution of the raster map is low and so does the resolution of the region successively) are not taken into account.

could you explain that more? I don't understand it.

can you provide an example using Spearfish's bugsites + elevation.dem or something from the NC standard dataset?

fwiw the region growing I've done has always been for cosmetic reasons. in those cases usually make the res=res*2 (or 10) with the -a flag, then run g.region again with the original region value.

Hamish

[hamish]

another trick is

eval `g.region -g`
res=$nsres
g.region n=n+$res s=s-$res w=w-$res e=e+$res

[nikos]

Replying to hamish:

did you try the g.region align=rastermap option?

fwiw, -a grows outwards.

Matching the region over a point map is the problem! The points at the borders (happens quite often to me when the resolution of the raster map is low and so does the resolution of the region successively) are not taken into account.

could you explain that more? I don't understand it.

1. Imagine a point-grid (i.e. points or centroids that would form a square if connected)
   
2. match the region over this point-grid (g.region -a vect=point-grid # the points at the edges are just to be seen when drawing the vector point-grid with d.vect for example).
   
3. The points at the edges are ignored by v.what.rast vector=point-grid rast=LowResRaster

can you provide an example using Spearfish's bugsites + elevation.dem or something from the NC standard dataset?

Right. I'll try to re-produce this in Spearfish.

fwiw the region growing I've done has always been for cosmetic reasons. in those cases usually make the res=res*2 (or 10) with the -a flag, then run g.region again with the original region value.

Would this (res=res*X) affect v.what.rast?

(Please, wait till I provide a Spearfish example. No reason to comment more on this. Also, the "trick" is an optimal solution I think. Have to try this out.)

[nikos]

can you provide an example using Spearfish's bugsites + elevation.dem or something from the NC standard dataset?

Right. I'll try to re-produce this in Spearfish.

# launch grass65 in spearfish60
nik@vertical:~$ grass65 /geo/grassdb/spearfish60/user1/

copy-paste-execute everything below!

# define a region
g.region n=4921000 s=4915000 w=591000 e=598500 res=500 -pa

# create a vector grid
v.mkgrid map="tempGrid" grid=12,15 position="region" angle=0 breaks=3 --o

# extract centroids from vector grid
v.extract input="tempGrid" output="tempGrid_centroids" type="centroid" layer=1 new=-1 --o

# draw maps
d.mon x0 && d.vect tempGrid && d.vect tempGrid_centroids col=blue

# convert centroids to points and draw
v.type in=tempGrid_centroids out=tempGrid_points type=centroid,point
d.vect tempGrid_points col=red

# add column to be updated from raster
v.db.addcol tempGrid_points column="landcover integer"

# match region with "-a" indeed, it "grows" outwards (or it just stays like it was before!?).
g.region vect=tempGrid_points -pa

# check visually: nothing changed
d.redraw

# now try _without_ "-a"
g.region vect=tempGrid_points -p

# obviously the region has been altered (as expected): check (also) the numbers!
d.redraw

# now try again with "-a"
g.region vect=tempGrid_points -pa

# there is no "returning back" or "growing"! (only when using res=500 again, keep reading)
d.redraw

# update points from landcover.30m raster and check table
v.what.rast tempGrid_points rast=landcover.30m column=landcover
db.select tempGrid_points
# works fine!

# assume that the raster was just about the same size of what the (grass-)region is:
r.mapcalc landcover_cropped = landcover.30m

# add a second column and update again
v.db.addcol map=tempGrid_points col="landcover2 integer"
v.what.rast vect=tempGrid_points rast=landcover_cropped col=landcover2

# check updated table
db.select tempGrid_points
# works fine too!


# now draw raster and points
d.erase
d.rast landcover_cropped
d.vect tempGrid_points col=red


# here is where the problem begins:
  # first check current region size
g.region -p
  # for some reason I was using the "res" parameter with g.region while looping over several point vector maps to get them updated from rasters
  # so, let's use res=500
g.region vect=tempGrid_points -pa res=500

# note the changes in the region:
  # the number of rows, cols. The region growed outwards (as expected)!
  # north is now larger
  # south has become smaller (!?)
  # west has become smaller (!?)
  # east is now larger 

# now repeat v.what.rast
v.what.rast vect=tempGrid_points rast=landcover_cropped col=landcover2

# check updated column 
db.select tempGrid_points
d.vect tempGrid_points col=blue where="landcover2 > 0"

# In this example I see points at east, south not being updated (actually updated with NULL)!
# Missing values only at (some) the edges! huh??

Nikos

[hamish]

Replying to hamish:

did you try the g.region align=rastermap option?

[hamish]

Hamish:

fwiw the region growing I've done has always been for cosmetic reasons. in those cases usually make the res=res*2 (or 10) with the -a flag, then run g.region again with the original region value.

that should read "original resolution value".

and unlike nsew=+-, that is multiplication by hand, not by the module.

Hamish

[nikos]

Replying to hamish:

Replying to hamish:

did you try the g.region align=rastermap option?

This works fine. Nevertheless, what I want to do is a bit more complicated:

• I have X raster and Y point vector maps.
  
• Each point vector map has one column for each raster, that is at least X columns (+cat +the typical rows, cols columns).
  
• I want to update the X columns in each point map with the values from the (X) raster maps (one column per raster map of course).
  

So what I do in my script is (in pseudocode):

for point_map in group_of_point_maps:
    g.region vect=point_map

        for Raster in group_of_rasters:
            v.what.rast vect=point_map rast=Raster column=Raster

This is where the problem arises since I can't use the align=ToRasterMap option. To make the long story short it is required:

1. to loop over a series of point vector maps whose columns are to be updated
   
2. to (sub-) loop over a series of raster maps from which the vector columns will be updated (for each point vector map).

How to solve this? Maybe it's a common problem but I don't know how to deal with it than just use a basic naming convention for all maps and play later with string manipulation utilities (such as cut, tr, sed, etc in bash or the likes in python). So, I was thinking that a "better" g.region option would make things easier when "align=ToRasterMap" is not an option.

Currently I don't have a problem because, as I described in the very beginning, I have the original vector grids from which the point maps derived and it's easy to use them (almost same name) in the "g.region" command in order to get a larger region.

Hope this is now clear. Thanks, Nikos

[hamish]

does the g.region.point addons script help?

​http://grass.osgeo.org/wiki/AddOns#g.region.point

Hamish

[hamish]

Replying to hamish:

does the g.region.point addons script help?

​http://grass.osgeo.org/wiki/AddOns#g.region.point

ping

[Nikos Alexandris]

Replying to hamish:

Replying to hamish:

does the g.region.point addons script help?

​http://grass.osgeo.org/wiki/AddOns#g.region.point

First guess is no. How would g.region.point help in setting the region so as to cover/match the extent of a point vector map and be a bit larger so as to include its points in the v.what.rast process? g.region.point sets the region around a pair of coordinates. Is the idea behind to loop over all points in a map? I guess not -- this would be a killer for loop.

[hamish]

Replying to nikosa:

How would g.region.point help in setting the region so as to cover/match the extent of a point vector map and be a bit larger so as to include its points in the v.what.rast process?

your initial idea of "g.region -a vect=mmmm res=xxxx" solves it. Set the resolution to 1000 there (or whatever) and the -a flag will always grow outwards in all 4 directions. So you run the above, then run g.region a second time with your finer resolution. If your second finer-resolution fits evenly in the initial coarser one, you now have a grown region with bounds rounded to nice whole numbers, a bit bigger than your points map so catching them all!

Hamish

[Nikos Alexandris]

Replying to hamish:

your initial idea of "g.region -a vect=mmmm res=xxxx" solves it. Set the resolution to 1000 there (or whatever) and the -a flag will always grow outwards in all 4 directions. So you run the above, then run g.region a second time with your finer resolution. If your second finer-resolution fits evenly in the initial coarser one, you now have a grown region with bounds rounded to nice whole numbers, a bit bigger than your points map so catching them all!

Maybe it does so. And, maybe I have placed, in the example above, the

# assume that the raster was just about the same size of what the (grass-)region is:
r.mapcalc landcover_cropped = landcover.30m

in the "wrong" place. This was part of my "pareto" implementation using MODIS data (low-res, to be assessed) and Landsat (higher-res, as a reference). I have to go through this again. If I remember well, I was forced to have "this" very step in-between of the grid creation/centroid extraction operations.

[martinl]

Is there anything we can do or just close this ticket?

[martinl]

Milestone renamed

[martinl]

Can we close this ticket as wontfix?

[neteler]

Ticket retargeted after milestone closed

[lucadelu]

I tried to implement the "expand the computational region of about "some" pixels" request, is my changes too simple?

[mmetz]

Replying to lucadelu:

I tried to implement the "expand the computational region of about "some" pixels" request, is my changes too simple?

The number of pixels is just a number, not a coordinate, therefore the new option should be of type TYPE_DOUBLE and scanned with sscanf(), not G_scan_easting()/G_scan_northing(). Should top and bottom also be adjusted?

[lucadelu]

Replying to mmetz:

The number of pixels is just a number, not a coordinate, therefore the new option should be of type TYPE_DOUBLE and scanned with sscanf(), not G_scan_easting()/G_scan_northing(). Should top and bottom also be adjusted?

Thanks Markus, I modify the patch with your suggestion, just one question, why TYPE_DOUBLE and not INTEGER (I used INT)? because the number of pixels should be integer, I don't think is good idea to permit to add 2.5 pixels ;-)

However I have no answer about top/bottom, I never used them. I try to run

g.region t=t+8

or

g.region t=800

but nothing changed in the output of g.region -p, so no idea...

[martinl]

Replying to lucadelu:

but nothing changed in the output of g.region -p, so no idea...

try

g.region -p3

[mmetz]

Replying to lucadelu:

Replying to mmetz:

The number of pixels is just a number, not a coordinate, therefore the new option should be of type TYPE_DOUBLE and scanned with sscanf(), not G_scan_easting()/G_scan_northing(). Should top and bottom also be adjusted?

Thanks Markus, I modify the patch with your suggestion, just one question, why TYPE_DOUBLE and not INTEGER (I used INT)? because the number of pixels should be integer, I don't think is good idea to permit to add 2.5 pixels ;-)

You are right, the number of pixels should be integer.

[lucadelu]

Replying to martinl:

Replying to lucadelu:

but nothing changed in the output of g.region -p, so no idea...

try

g.region -p3

Thanks Martin, top and bottom are used only for 3D region? If yes I don't think that pixels option should modify them....

[mlennert]

Replying to lucadelu:

Replying to martinl:

Replying to lucadelu:

but nothing changed in the output of g.region -p, so no idea...

try

g.region -p3

Thanks Martin, top and bottom are used only for 3D region? If yes I don't think that pixels option should modify them....

Why not ? If you have the same use case as the one described by Nikos, but with a 3D point map, wouldn't you need the same functionality also in top and bottom direction ?

[lucadelu]

Replying to mlennert:

Why not ? If you have the same use case as the one described by Nikos, but with a 3D point map, wouldn't you need the same functionality also in top and bottom direction ?

For me it is the same, but probably I would like to have a flag to modify the 3D region to use with pixels option

my2cents

[mmetz]

Replying to lucadelu:

Replying to mlennert:

Why not ? If you have the same use case as the one described by Nikos, but with a 3D point map, wouldn't you need the same functionality also in top and bottom direction ?

For me it is the same, but probably I would like to have a flag to modify the 3D region to use with pixels option

Either have one option that modifies extents of all 3 dimensions, or have separate options for each dimension, e.g. xpixels, ypixels, zpixels. I favour the first solution, one option for all 3 dimensions.

[lucadelu]

Replying to mmetz:

Either have one option that modifies extents of all 3 dimensions, or have separate options for each dimension, e.g. xpixels, ypixels, zpixels. I favour the first solution, one option for all 3 dimensions.

Ok, implemented in r72274. If it works correctly please close the ticket

[Nikos Alexandris]

I am out of the loop here. g.region pixels=1 goes from

rows:       1193
cols:       981

to

rows:       1195
cols:       983

This means expansion, by the given number of pixels=1 in all directions. It works for me, thank you.

[lucadelu]

Closing. Please reopen if needed

[wenzeslaus]

Should this be g.region pixels= or cells= or something else?

Pixel is not really used in GRASS overall, see search in the HTML doc folder:

$ cd dist.../docs/html
$ grep -Irn cell | wc
   1360   16896  162701
$ grep -Irn cells | wc
    654    8520   79660
$ grep -Irn pixel | wc
    200    3035   32650
$ grep -Irn pixels | wc
    108    1169    9922
$ grep -Irn cell | grep -E "^i\." | wc
     35     400    3083
$ grep -Irn pixel | grep -E "^i\." | wc
    103    1115   11107
# (the singular search contains plural too)

Cell(s) seems to be much more common than pixel(s). And half of the pixel matches seems to be in imagery modules (which is not the case for cell). g.region manual page (G7:g.region) does not contain the word pixel and neither does the raster intro page (G7:rasterintro). (In case of the intro it should, but that's another story.)

Based on that it seems that cells should be used instead of pixels but cells is in the g.region output telling us the total number of cell, so perhaps that would be confusing and not enough descriptive. (Is it total cells/pixels or something your are adding or removing?) Looking at the keywords here, expand, extend, enlarge, increase, add, make_larger, broaden, plus, and grow seems possible. My suggestion is:

enlarge=value
    Number of cells to enlarge the bounding box

This does not account for the opposite case when a negative number is provided (from the code I assume it works) but that can be implicit or handled by a detailed description (using label and description).

[mmetz]

Replying to wenzeslaus:

Should this be g.region pixels= or cells= or something else?

Pixel is not really used in GRASS overall, see search in the HTML doc folder:

$ cd dist.../docs/html
$ grep -Irn cell | wc
   1360   16896  162701
$ grep -Irn cells | wc
    654    8520   79660
$ grep -Irn pixel | wc
    200    3035   32650
$ grep -Irn pixels | wc
    108    1169    9922
$ grep -Irn cell | grep -E "^i\." | wc
     35     400    3083
$ grep -Irn pixel | grep -E "^i\." | wc
    103    1115   11107
# (the singular search contains plural too)

Cell(s) seems to be much more common than pixel(s). And half of the pixel matches seems to be in imagery modules (which is not the case for cell). g.region manual page (G7:g.region) does not contain the word pixel and neither does the raster intro page (G7:rasterintro). (In case of the intro it should, but that's another story.)

Based on that it seems that cells should be used instead of pixels but cells is in the g.region output telling us the total number of cell, so perhaps that would be confusing and not enough descriptive. (Is it total cells/pixels or something your are adding or removing?) Looking at the keywords here, expand, extend, enlarge, increase, add, make_larger, broaden, plus, and grow seems possible. My suggestion is:

enlarge=value
    Number of cells to enlarge the bounding box

This does not account for the opposite case when a negative number is provided (from the code I assume it works) but that can be implicit or handled by a detailed description (using label and description).

How about grow, shorter than enlarge, also in line with r.grow, r.grow.distance?

grow=value
label: Number of cells to add to each side of the current region extents
description: A negative number shrinks the current region extents

[wenzeslaus]

In 72296:

g.region: pixels not really used, grow says what it does (wording by mmetz, see #789)

[wenzeslaus]

I did the change in r72296. Here is a diff of only the naming change:

main.c

@@ -337 +337 @@
     parm.pixels = G_define_option();
-    parm.pixels->key = "pixels";
+    parm.pixels->key = "grow";
     parm.pixels->key_desc = "value";
     parm.pixels->required = NO;
     parm.pixels->multiple = NO;
     parm.pixels->type = TYPE_INTEGER;
-    parm.pixels->description =  
-    _("Number of pixels to increase the bounding box");
+    parm.pixels->label =
+        _("Number of cells to add to each side of the current region extent");
+    parm.pixels->description =
+        _("A negative number shrinks the current region extent");
     parm.pixels->guisection = _("Bounds");

The is an issue with shrinking now reported in #3509.