1

I wrote the following code to use r.water.outlet module in GRASS GIS:

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

I stored the coordinates in tc_sta.csv, and input them to r.water.outlet to calculate drainage area of each point.

Although the output raster files were fine, I had to export them as ESRI shapefile and use QGIS to calculate area of polygons, and finally store the area value of each point in a csv file containing coordinates and its related drainage area.

I felt that the procedure was not efficient, and is there a way I could calculate and store the output area immediately after I run r.water.outlet?

1

You can use r.report on the outlet map right away after the r.water.outlet step. It calculates the area size in a user defined unit (m^2, ha, etc.)

  • It's a nice advice! I am wondering if I use python to save the area value r.report outputs in a CSV file? – Heinz Oct 21 '16 at 9:23
  • I tried r.stats and it also worked! Thank you for the advice, could you help me with this post, thanks? gis.stackexchange.com/questions/215075/… – Heinz Oct 22 '16 at 11:28
  • Here is existing GRASS GIS Python code for parsing the r.report output: trac.osgeo.org/grass/browser/grass-addons/grass7/raster/r.green/… which you may use and adapt. – markusN Oct 23 '16 at 11:20
  • thank you, but how can I use this? I downloaded utils.py and put it in C:\OSGeo4W64\apps\grass\grass-7.0.5\scripts, but I still can't run perc_of_overlay in GRASS. And I still can not store just the area value. – Heinz Oct 23 '16 at 12:58
1

I routinely use this script:

#!/bin/bash
if [[ $# -lt 3  ]]; then
        echo "Parameters missing"
        exit 1
else
        BASIN=$1
        X=$2
        Y=$3
fi
g.region --q rast=<flow_dir>
r.stream.basins --o direction=<flow_dir> coordinates="${X},${Y}" basin=${BASIN}

r.to.vect --o -s $BASIN output=${BASIN} type=area
v.db.addcolumn $BASIN column="area_sqkm double"
v.db.update $BASIN column=label value=${BASIN}
v.to.db $BASIN option=area unit=kilometers column=area_sqkm
[[ -f ${BASIN}.shp ]] && rm -f ${BASIN}.*
v.out.ogr -s -e $BASIN output=${BASIN}.shp

Maybe it will be of some help.

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.