1

I type those code in pyscripter,

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'acc1', drainge = 'dra1', basin = 'bas1', stream = 'str1', half_basin = 'hbas1', visual = 'vis1')
print g.read_command('g.list', _type = 'rast')

PS. I set the output half.basin as half_basin because the following error will occur as I set it as former one,

Message File Name   Line    Position    
SyntaxError             
    keyword can't be an expression (module2.py, line 27)    C:\Users\Heinz\Documents\pyscript\module2.py    27

and here's the output,

enter image description here

there are no new raster file generated, and I don't know where I make mistake.


Update#1 This is my code:

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'accu1', drainge = 'dran1', basin = 'basi1', stream = 'stre1', half.basin = 'hbasi1', visual = 'visu1')

#check result
print g.read_command('g.list', _type = 'rast')

I use drainage and half.basin as my parameter, and here's the output:

enter image description here


Update#2 I tried the script below,

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'water_test_2'
gsetup.init(gisbase, gisdb, location, mapset)

#import raster data
#g.run_command('r.in.gdal', input = 'hc', out = 'hc')

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'accu', drainage = 'dran', basin = 'basi', stream = 'stre', half = 'hbas', visual = 'visu')

#check result
print g.read_command('g.list', _type = 'rast')

And the output looks fine,

*** Remote Interpreter Reinitialized  ***
>>> 
----------------------------------------------
raster files available in mapset <water_test_2>:
accu basi dran hbas stre visu

raster files available in mapset <PERMANENT>:
hc

----------------------------------------------

My conclusion is, we can use python to call r.watershed in GRASS 6.4.3 from outside, but one of the input name must be 'half'

Thanks everyone that helps me successfully get the result I want.

3

Documentation of r.watershed module for GRASS 6 says that there is a half.basin option (parameter):

half.basin   Output map: each half-basin is given a unique value

So, this means that in command line (Bash etc.) one can write:

r.watershed elevation=elev_lid792_1m drainage=elev_lid_drainage half.basin=elev_lid_half_basin threshold=10

GRASS provides a clear way how to call the GRASS module in Python similarly to Bash (e.g., using run_command function). Unfortunately, Python does not allow a dot (.) in function keyword parameter name, so using half.basin causes an (Python) error:

>>> gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', half.basin='elev_lid_half_basin', threshold=10)
  File "<stdin>", line 1
SyntaxError: keyword can't be an expression

Writing there underscore (_) or nothing instead of a dot is a good idea, however this replacement would have to be supported by the module or GRASS in general. And since this was not implemented, it results in error:

>>> gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', half_basin='elev_lid_half_basin', threshold=10)
Sorry, <half_basin> is not a valid parameter
...

Fortunately, what is possible and implemented in GRASS 6, is the possibility to shorten the name of the parameter, so you can write half instead of half.basin (as long as it is not causing ambiguity with some other parameter). So, in GRASS 6 you can write the following and it will work:

>>> gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', half='elev_lid_half_basin', threshold=10)

For GRASS GIS 7, the module options were revised and standardized and dot should not be a part of the module option (parameter) anymore (compatibility with Python is one of the reasons). As a consequence, r.watershed manual for GRASS 7 says:

half_basin   Name for output half basins raster map
               Each half-basin is given a unique value

So, in GRASS 7, in Python, you can really write half_basin and it will work:

 >>> gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', half_basin='elev_lid_half_basin', threshold=10)

Moreover, GRASS 7 allows advanced shortening of module options, so you can also use halfbasin or hbasin and it will work:

>>> gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', halfbasin='elev_lid_half_basin', threshold=10)

In case of r.watershed, you cannot shorten half_basin to basin since it would conflict with the basin option (basin Unique label for each watershed basin).

To get all outputs from r.watershed (in GRASS 6 or 7) you can use the following (used map names works in the GRASS NC sample dataset):

gcore.run_command('r.watershed', elevation='elev_lid792_1m', drainage='elev_lid_drainage', half='elev_lid_half_basin', accumulation='elev_lid_accumulation', basin='elev_lid_basin', stream='elev_lid_stream', visual='elev_lid_visual', length='elev_lid_length_slope', slope='elev_lid_slope_steepness', threshold=10)

In my opinion, the best practice is to not to use shortened options in scripts and use only the full option names. This means that (in both GRASS 6 and GRASS 7) you should use the shortened options only when writing commands in command line (Bash or Python) but not when writing scripts (Bash or Python).

However, there is one exception. When the shortened option is the way to achieve compatibility between (current stable release) GRASS GIS 6 and GRASS GIS 7 (currently beta release), I think that you can use it. So, in our example of r.watershed case you can use half (which will mean half.basin in GRASS 6 and half_basin in GRASS 7). But I think that this is the only when shortened options should be used in scripts. For sure, an additional comment with an explanation should be provided in this case.

  • thank you, so you mean these shortened parameters can't be used in python script? – Heinz Mar 31 '14 at 3:52
  • I mean that the best practice is not to use them, although it is possible to use them. I edited my answer to be more clear. Please, accept the answer if you are satisfied. – wenzeslaus Mar 31 '14 at 15:12
  • I use the shortened parameter in my python script but still can't generate 6 output files from r.watershed. GRASS can't just take the position of parameters, no matter what their name is, to execute the module. – Heinz Mar 31 '14 at 16:31
  • 1
    To your question about just providing parameter values and not the names, the answer is no, you must provide both names and values. The reason for this limitation is that it improves readability of the commands and it avoids confusion when there is a lot of options and you would accidentally change their order. – wenzeslaus Mar 31 '14 at 17:38
  • 1
    According to your other other emails, I would say that you should double check how you are creating a location and mapset and how you are switching to them. Explore your locations, mapsets and maps using some other method, e.g. GRASS GIS GUI. Also check outputs of g.list in different mapsets to understand which data you have there and what g.list is showing. If all of that fails, ask a new question (this no longer looks like r.watershed interface problem). – wenzeslaus Apr 1 '14 at 3:16
2

You are apparently not running GRASS in its own session (fine) but you are working with GRASS without starting it explicitly. There is a Python example for this.

Concerning your example: check the command line output for the error message (unfortunately not seen in the Python shell).

I tried and got these errors listed:

Sorry, <drainge> is not a valid parameter

--> "...drainge = 'dra1'..." should be "drainage"

Sorry, <half_basin> is not a valid parameter

--> "...half_basin = 'hbas1'..." should be "half.basin"

So you have simply two typos there.

  • I don't know what you mean, because I have already set drainage and halfbasin(half.basin will cause error) as my parameters. – Heinz Mar 29 '14 at 13:06
  • Sure but you wrote "drainge" and that is a typo. Can you post the error of "half.basin"? Then I can check with the developers. – markusN Mar 29 '14 at 16:28
  • I have edited my post. – Heinz Mar 30 '14 at 4:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.