Hi
I am trying the generate watersheds for several outlets by looping
over r.water.outlet in a python script (attached). The basin procedure
is that for each outlet point (read using v.info -t), my script runs
r.water.outlet, r.to.vect and v.patch, to combine the watersheds into
a simgle vector map. I am connecting the watersheds to the outlet
points by setting the category value in the output from r.to.vect,
i.e. before patching. The trouble is that some of the watersheds
overlap, e.g. where an outlet lies inside another watershed. v.patch
doesn't seem to handle this too well, and the result is that I have a
vector map with a few areas, a few boundaries, a few centroids and a
lot of errors.
Is there a better way of approaching this? And how do I allow overlapping areas?
Regards
David
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
############################################################################
#
# MODULE: r.water.outlet.points
# AUTHOR(S): David Townshend
# PURPOSE: Runs r.water.outlet for a list of points
# COPYRIGHT: (C) 2004,2008,2009 by the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
#############################################################################
#%Module
#% description: Runs r.water.outlet for a list of points
#%End
#%option
#% key: drainage
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of input raster map
#% required: yes
#%end
#%option
#% key: basin
#% type: string
#% gisprompt: new,vector,vector
#% description: Name of vector map to contain results
#% required: yes
#%end
#%option
#% key: points
#% type: string
#% gisprompt: old,vector,vector
#% description: Points at which to calculate catchments
#% required: yes
#%end
import sys
from grass.script import core as gg
from grass.script import vector as gv
from grass.script import raster as gr
def getColumnType(map, column):
info = gg.parse_command('v.info', flags='c', map=map, quiet=True).keys()
dinfo = {}
for i in info:
t, n = i.split('|')
dinfo[n] = t
return dinfo[column]
def percent(i, n, s):
# TODO: Figure out why percent doesn't work
gg.message("%d %d %d" % (i, n, s))
def main():
gg.message('Initializing')
# Get options
r_drainage = options['drainage']
v_basin = options['basin']
v_points = options['points']
total_points = int(gg.parse_command('v.info', flags='t', map=v_points, quiet=True)['points'])
# Create new vector
gg.parse_command('v.edit', tool='create', map=v_basin, quiet=True)
# Get coordinates
coords = gg.parse_command('v.to.db', flags='p', map=v_points, option='coor', type='point', quiet=True).keys()
count = 0
gg.message('Calculating')
for line in coords:
percent(count, total_points, 1)
# Create temporary files
r_basin = 'tmp_r_basin'
v_tmp_basin = 'tmp_v_basin'
v_tmp_basin2 = 'tmp_v_basin2'
# Get point info
coord = line.split('|')
cat, x, y = int(coord[0]), float(coord[1]), float(coord[2])
# Calculate watershed
gg.parse_command('r.water.outlet', drainage=r_drainage, basin=r_basin, easting=x, northing=y)
# Convert to vector
gg.parse_command('r.to.vect', flags='v', input=r_basin, output=v_tmp_basin, feature='area', quiet=True)
# Update vectors
gg.parse_command('v.edit', map=v_tmp_basin, cats=0, tool='delete', quiet=True)
gg.parse_command('v.category', input=v_tmp_basin, output=v_tmp_basin2, option='sum', cat=cat-1, quiet=True)
gg.parse_command('v.patch', flags='a', input=v_tmp_basin2, output=v_basin, overwrite=True, quiet=True)
# Clean
gg.parse_command('g.remove', rast=r_basin, vect=','.join([v_tmp_basin, v_tmp_basin2]), quiet=True)
count += 1
percent(count, total_points, 1)
if __name__ == '__main__':
options, flags = gg.parser()
sys.exit(main())
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user
