Friday, 4 September 2015

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


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:


enter image description here


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



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')

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...