[GRASS-dev] [GRASS GIS] #789: g.region option to expand the computational region of about "some" pixels?

GDAL · October 14, 2009, 10:46pm

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@lists.osgeo.org
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------
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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 3:12am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:1>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 3:14am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by hamish):

another trick is

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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:2>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 8:18am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by nikos):

Replying to [comment:1 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.

a. Imagine a point-grid (i.e. points or centroids that would form a square
if connected)[[BR]]
b. 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).[[BR]]
c. 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?[[BR]]

(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.)

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789#comment:3>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 10:40am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by 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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789#comment:4>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 10:57am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by hamish):

Replying to [comment:1 hamish]:
> did you try the g.region align=rastermap option?

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:5>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 11:00am

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:6>
GRASS GIS <http://grass.osgeo.org>

GDAL · October 15, 2009, 12:07pm

#789: g.region option to expand the computational region of about "some" pixels?
--------------------------+-------------------------------------------------
  Reporter: nikos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: normal | Milestone: 7.0.0
Component: default | Version: unspecified
Resolution: | Keywords: g.region, expand computational region
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by nikos):

Replying to [comment:5 hamish]:
> Replying to [comment:1 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.[[BR]]
  * Each point vector map has one column for each raster, that is at least
X columns (+cat +the typical rows, cols columns).[[BR]]
  * 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).[[BR]]

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:[[BR]]
a. to loop over a series of point vector maps whose columns are to be
updated[[BR]]
b. 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.[[BR]]

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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789#comment:7>
GRASS GIS <http://grass.osgeo.org>

GDAL · March 9, 2011, 7:25am

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------

Comment(by hamish):

does the g.region.point addons script help?
http://grass.osgeo.org/wiki/AddOns#g.region.point

Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:8>
GRASS GIS <http://grass.osgeo.org>

GDAL · November 12, 2013, 9:53am

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------

Comment(by hamish):

Replying to [comment:8 hamish]:
> does the g.region.point addons script help?
> http://grass.osgeo.org/wiki/AddOns#g.region.point

ping

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789#comment:9>
GRASS GIS <http://grass.osgeo.org>

GDAL · November 14, 2013, 8:13pm

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------

Comment(by nikosa):

Replying to [comment:9 hamish]:
> Replying to [comment:8 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/789#comment:10>
GRASS GIS <http://grass.osgeo.org>

GDAL · November 14, 2013, 8:53pm

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------

Comment(by hamish):

Replying to [comment:10 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:11>
GRASS GIS <http://grass.osgeo.org>

GDAL · November 15, 2013, 9:03am

#789: g.region option to expand the computational region of about "some" pixels?
---------------------------------------------------+------------------------
Reporter: nikos | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: unspecified
Keywords: g.region, expand computational region | Platform: Unspecified
Cpu: Unspecified |
---------------------------------------------------+------------------------

Comment(by nikosa):

Replying to [comment:11 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.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/789#comment:12>
GRASS GIS <http://grass.osgeo.org>