[GRASS-user] how to g.region

[Michael Barton]

On 9/16/07 1:29 AM, "Raffaello Brondi" <raffaello.brondi at 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 :( ) 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