[GRASS-user] how to g.region

GRASS GRASS user
1

Hi all,

i'm using GRASS 6.2 for developing a Web Processing Service.
One of the operations offered by the service takes two raster maps as inputs.
The raster map resulting from this operation would have different spatial extents based on the input maps: the result would have the extent equal to the extent of one of the two input maps, or equal to the union of the extents of the two input maps or equal to the intersection of the extents of the two raster maps.
In order to develop the service i need to manage these four cases by changing the GRASS region.
Supposing that the input rasters are "raster1" and "raster2", the possible configurations are:

   1. region = raster1 --> g.region rast=raster1
   2. region = raster2 --> g.region rast=raster2
   3. region = the union of raster1 and raster2 --> g.region rast=raster1,raster2
   4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify the last point except to manually set the nord,sud,ovest and east values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region to set the current region as the intersection of two rasters?

Best regards

Raffaello Brondi

2

On 12.09.2007 20:01, Raffaello Brondi wrote:

  4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify
the last point except to manually set the nord,sud,ovest and east
values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region
to set the current region as the intersection of two rasters?

Hello,

I assume you have some kind of programing language which executes grass
commands. Can you parse the output? If so then you could do it like this:

g.region rast=raster1
g.region -p (parse the output, into rast1north, rast1east, rast1west,
rast1south, etc)
g.region rast=raster2
g.region -p (parse the output, into rast2north, rast2east, rast2west,
rast2south, etc)
g.region north=(min of rast1north, rast2northe) etc for other parameters.

Hope this helps,
--Wolf

--

<:3 )---- Wolf Bergenheim ----( 8:>

3

Raffaello Brondi wrote:

i'm using GRASS 6.2 for developing a Web Processing Service.

Have you seen Jachym's PyWPS?
  http://pywps.wald.intevation.org

One of the operations offered by the service takes two raster maps as
inputs.
The raster map resulting from this operation would have different
spatial extents based on the input maps: the result would have the
extent equal to the extent of one of the two input maps, or equal to
the union of the extents of the two input maps or equal to the
intersection of the extents of the two raster maps.
In order to develop the service i need to manage these four cases by
changing the GRASS region.
Supposing that the input rasters are "raster1" and "raster2", the
possible configurations are:

   1. region = raster1 --> g.region rast=raster1
   2. region = raster2 --> g.region rast=raster2
   3. region = the union of raster1 and raster2 --> g.region
rast=raster1,raster2
   4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify
the last point except to manually set the nord,sud,ovest and east
values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region
to set the current region as the intersection of two rasters?

This is not very efficient, and more like 'g.region zoom=', but you could
do:

g.region rast=map1,map2
r.mapcalc 'intersect= if(!isnull(map1) && !isnull(map2), 1, null())'
g.region zoom=intersect
g.remove intersect

But probably Wolf's method is what you want.

Hamish

4

Wolf Bergenheim wrote:

On 12.09.2007 20:01, Raffaello Brondi wrote:
  

  4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify
the last point except to manually set the nord,sud,ovest and east
values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region
to set the current region as the intersection of two rasters?
    
Hello,

I assume you have some kind of programing language which executes grass
commands. Can you parse the output? If so then you could do it like this:

g.region rast=raster1
g.region -p (parse the output, into rast1north, rast1east, rast1west,
rast1south, etc)
g.region rast=raster2
g.region -p (parse the output, into rast2north, rast2east, rast2west,
rast2south, etc)
g.region north=(min of rast1north, rast2northe) etc for other parameters.

Hope this helps,
--Wolf

Hello Wolf

Yes, the service i'm developing is written in Java.
Your solution is a good one but the problem is the overhead it introduces.
Since i can't (or maybe i'm not able to :frowning: ) interact directly with GRASS using Java, i would need to make 4 files reading/writing.

thanx

Raffaello Brondi

5

Hi Hamish.

Raffaello Brondi wrote:
  

i'm using GRASS 6.2 for developing a Web Processing Service.
    
Have you seen Jachym's PyWPS?
  http://pywps.wald.intevation.org
  

yes i saw it..but unfortunately i have to develop the service using Java.

  

One of the operations offered by the service takes two raster maps as inputs.
The raster map resulting from this operation would have different spatial extents based on the input maps: the result would have the extent equal to the extent of one of the two input maps, or equal to
the union of the extents of the two input maps or equal to the
intersection of the extents of the two raster maps.
In order to develop the service i need to manage these four cases by changing the GRASS region.
Supposing that the input rasters are "raster1" and "raster2", the possible configurations are:

   1. region = raster1 --> g.region rast=raster1
   2. region = raster2 --> g.region rast=raster2
   3. region = the union of raster1 and raster2 --> g.region rast=raster1,raster2
   4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify the last point except to manually set the nord,sud,ovest and east values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region to set the current region as the intersection of two rasters?
    
This is not very efficient, and more like 'g.region zoom=', but you could
do:

g.region rast=map1,map2
r.mapcalc 'intersect= if(!isnull(map1) && !isnull(map2), 1, null())'
g.region zoom=intersect
g.remove intersect

But probably Wolf's method is what you want.

Hamis

This is the solution i need :).
thanx a lot!

Raffaello Brondi

6

Raffaello Brondi wrote:

>> i'm using GRASS 6.2 for developing a Web Processing Service.
>
> Have you seen Jachym's PyWPS?
> http://pywps.wald.intevation.org
>
yes i saw it..but unfortunately i have to develop the service using
Java.

Have you seen jGRASS? http://www.jgrass.org/

Raffaello:

>> One of the operations offered by the service takes two raster maps
>> as inputs.
>> The raster map resulting from this operation would have different
>> spatial extents based on the input maps: the result would have the
>> extent equal to the extent of one of the two input maps, or equal to
>> the union of the extents of the two input maps or equal to the
>> intersection of the extents of the two raster maps.
>> In order to develop the service i need to manage these four cases by
>> changing the GRASS region.
>> Supposing that the input rasters are "raster1" and "raster2", the
>> possible configurations are:
>>
>> 1. region = raster1 --> g.region rast=raster1
>> 2. region = raster2 --> g.region rast=raster2
>> 3. region = the union of raster1 and raster2 --> g.region
>> rast=raster1,raster2
>> 4. region = the intersection of raster1 and raster2 --> ?
>>
>>
>> As far as i see, using the g.region command, there is no way to
>> specify the last point except to manually set the nord,sud,ovest and
>> east values, but unfortunately i can not use this way.
>> Is there any other command in GRASS that i can use instead of
>> g.region to set the current region as the intersection of two
>> rasters?
>

Hamish:

> This is not very efficient, and more like 'g.region zoom=', but you
> could do:
>
> g.region rast=map1,map2
> r.mapcalc 'intersect= if(!isnull(map1) && !isnull(map2), 1, null())'
> g.region zoom=intersect
> g.remove intersect

Raffaello:

This is the solution i need :).

Wolf:

g.region rast=raster1
g.region -p (parse the output, into rast1north, rast1east, rast1west,
rast1south, etc)
g.region rast=raster2
g.region -p (parse the output, into rast2north, rast2east, rast2west,
rast2south, etc)
g.region north=(min of rast1north, rast2northe) etc for other
parameters.

Raffaello:

Your solution is a good one but the problem is the overhead it
introduces. Since i can't (or maybe i'm not able to :frowning: ) interact
directly with GRASS using Java, i would need to make 4 files
reading/writing.

Running g.region twice is going to have way less overhead than running
r.mapcalc, especially if the maps are of a reasonable size.

GRASS now has SWIG bindings to help closely link GRASS's libgis etc.
into higher level languages than C/C++. Right now we have some
parts in place for perl & python, but a motivated contributor could
add java to that list. look in the swig/ directory in the 6.3 source code.
  http://www.swig.org/Doc1.3/Java.html

Hamish

7

On 9/16/07 1:29 AM, "Raffaello Brondi" <raffaello.brondi@pisa.intecs.it>
wrote:

Wolf Bergenheim wrote:

On 12.09.2007 20:01, Raffaello Brondi wrote:
  

  4. region = the intersection of raster1 and raster2 --> ?

As far as i see, using the g.region command, there is no way to specify
the last point except to manually set the nord,sud,ovest and east
values, but unfortunately i can not use this way.
Is there any other command in GRASS that i can use instead of g.region
to set the current region as the intersection of two rasters?
    
Hello,

I assume you have some kind of programing language which executes grass
commands. Can you parse the output? If so then you could do it like this:

g.region rast=raster1
g.region -p (parse the output, into rast1north, rast1east, rast1west,
rast1south, etc)
g.region rast=raster2
g.region -p (parse the output, into rast2north, rast2east, rast2west,
rast2south, etc)
g.region north=(min of rast1north, rast2northe) etc for other parameters.

Hope this helps,
--Wolf

Hello Wolf

Yes, the service i'm developing is written in Java.
Your solution is a good one but the problem is the overhead it introduces.
Since i can't (or maybe i'm not able to :frowning: ) interact directly with
GRASS using Java, i would need to make 4 files reading/writing.

Here is probably the fastest way to do this, using Python syntax (I think
I've got it right, but it IS Monday morning), which is similar to Java
(which I don't know).

# use -ugp to output region parameters in an easy to parse format without
# changing the region
map1 = subprocess.Popen(cmd=['g.region', '-ugp', 'rast=rast1'])
map2 = subprocess.Popen(cmd=['g.region', '-ugp', 'rast=rast2'])

# parse g.region output for each raster
map1 = map1.split('\n')
map2 = map1.split('\n')

for coord in map1:
    if 'n=' in coord:
        north1 = coord.split('=')[1]
    elif 's=' in coord:
        south1 = coord.split('=')[1]
    elif 'e=' in coord:
        east1 = coord.split('=')[1]
    elif 'w=' in coord:
        west1 = coord.split('=')[1]

for coord in map2:
    if 'n=' in coord:
        north2 = coord.split('=')[1]
    elif 's=' in coord:
        south2 = coord.split('=')[1]
    elif 'e=' in coord:
        east2 = coord.split('=')[1]
    elif 'w=' in coord:
        west2 = coord.split('=')[1]

# create a new region that is the intersection of rast1 and rast2
subprocess.Popen(cmd=['g.region', 'n=%s' % min(north1, north2), \
    's=%s' % max(south1, south2), 'e=%s' % min(east1, east2), \
    'w=%s' % max(west1, west2) ]

Michael
__________________________________________
Michael Barton, Professor of Anthropology
Director of Graduate Studies
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

8

Raffaello Brondi wrote:

> This is not very efficient, and more like 'g.region zoom=', but you could
> do:
>
> g.region rast=map1,map2
> r.mapcalc 'intersect= if(!isnull(map1) && !isnull(map2), 1, null())'
> g.region zoom=intersect
> g.remove intersect
>
>
> But probably Wolf's method is what you want.

This is the solution i need :).

No it isn't. You should be using Wolf's approach and computing the
intersection yourself.

However, rather than using "g.region rast=... ; g.region -p", it's
simpler to just use "r.info -sg map=...".

Creating a whole new map just so that you can find its bounds is
gratuitously inefficient.

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