How to do the loop when using GRASS module r.water.outlet?

Question

I have seen this post:Watershed Analysis with GRASS, and I have a .txt file which contains coordinates of many points, and I want to use r.water.outlet to get basin area of each point.

I type the following code in the post:

while read X Y; do
.... ;
done < cross_pts.txt

and in GRASS command line, it shows:

I already import the DEM file and successfully run r.watershed to get the drainage file in GRASS, but I do not know how to follow the guide in the post.

I am new to using command line and know nothing about bash script. How can I loop through all different coordinates? I am using GRASS 6.4.3 under windows8.

Answer 1

It looks like you're running GRASS using the default Windows shell. The loop code you posted is for a UNIX-like shell, and the Windows shell doesn't understand it. You have two options:

Use the Windows shell commands for all of your batch processing. This question/answer might help (although I've not done it myself): https://stackoverflow.com/questions/155932/how-do-you-loop-through-each-line-in-a-text-file-using-a-windows-batch-file

Or you could use a UNIX-like shell on Windows. I think that msys/mingw comes with the OSGeo4W installer so you probably already have it installed. In the command prompt you have there try typing "sh" (to start the new shell) and try your command again). If that doesn't work, try starting msys from the Start menu and starting GRASS from that shell.

Answer 2

here's my python script, and it runs fine,

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)

#csv file reading and importation
import csv
rows = list(open('tc_sta.csv'))
totalrows = len(rows) - 1

#loop through the csv(coordinates) file in r.water.outlet module
f = open('tc_sta.csv', 'r')
element = list(csv.reader(f))
i = 0
j = 0
while True:
    if i <= totalrows:
        g.run_command('r.water.outlet', drainage = 'dra', basin = 'b' + str(i + 1), east = element[i][j], north = element[i][j + 1])
        i = i + 1
    else:
        break

#print all files in grass database to check results of the module
print g.read_command('g.list', _type='rast')