GRASS-Wiki - User contributions [en]

Python/pygrass

2015-12-19T17:07:59Z

⚠️Pmav99: /* Interface to listing maps (g.list/g.mlist) */

''PyGRASS'' is a library that extends the GRASS GIS 7 capabilities to allow users to access the low-level GRASS API.

== Documentation ==

Detailed PyGRASS documentation: http://grass.osgeo.org/grass71/manuals/libpython/

== How to test library ==

''pygrass'' has 'doctest' included in its code. You can run doctest in the [http://grass.osgeo.org/sampledata/north_carolina/nc_basic_spm_grass7.tar.gz North Carolina basic location], using the mapset '''user1'''.

To test a single module (file) you have to change into the source code directory of pygrass and run the following code (this is an example of functions.py):

<source lang="bash">
python -m doctest functions.py
</source>

== References ==

* Zambelli, P., Gebbert, S., Ciolli, M., 2013. ''Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS)''. ISPRS International Journal of Geo-Information 2, 201–219. ([http://dx.doi.org/10.3390/ijgi2010201 DOI] | [http://www.mdpi.com/2220-9964/2/1/201/pdf PDF])

* [http://grass.osgeo.org/grass71/manuals/libpython/ pyGrass documentation] (updated weekly from SVN)

== Sample PyGRASS scripts ==

=== Training material ===

Two workshops were presented during the [http://geomorfolab.arch.unige.it/genova2013/ XIV Italian meeting of GRASS] at Genova.
The workshops use '''ipython notebook''' to show the python and pygrass API, all the material are available on github ([https://github.com/zarch/workshop-python python], [https://github.com/zarch/workshop-pygrass pygrass]). All the execute examples are reported here:

* python training (Genova)
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/00_intro.ipynb Introduction]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/01_data_types.ipynb Data Types]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/02_syntax.ipynb Syntax]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/03_function.ipynb Function]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/04_class.ipynb Class]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/05_other.ipynb Other]

* python training (EURAC)
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/00_intro.ipynb Introduction]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/01_data_types.ipynb Data Types]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/02_syntax.ipynb Syntax]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/03_version_control_system.ipynb Version Control System (hg)]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/04_objects.ipynb Objects]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/05_python_debugger.ipynb Debugger]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/06_functions.ipynb Functions]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/07_numpy.ipynb Numpy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/08_matplotlib.ipynb Matplotlib]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/09_scipy.ipynb Scipy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/10_pandas.ipynb Pandas]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/11_sympy.ipynb Sympy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/12_user_interfaces.ipynb GUI]

* pygrass training
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/00_Modules.ipynb Modules]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/01_GIS_objects.ipynb GIS]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/02_Vector.ipynb Vector]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/03_Raster.ipynb Raster]

=== Interface to listing maps (g.list/g.mlist) ===

Here we call a GRASS module that writes to stdout and do not call a Python function that returns a Python object, therefore we can save stdout and then parse it with:

<source lang="python">
from grass.pygrass.modules.shortcuts import general as g
import subprocess as sub
gl = g.mlist # GRASS 6
gl = g.list # GRASS 7
gl(type='rast', pattern='elev-*', stdout=sub.PIPE)
gl.outputs.stdout.split()
['elev-000-000', 'elev-000-001', 'elev-000-002', 'elev-001-000',
'elev-001-001', 'elev-001-002', 'elev-002-000', 'elev-002-001',
'elev-002-002']
</source>

or you can use the "glist" method of the Mapset class:

<source lang="python">
from grass.pygrass.gis import Mapset
m = Mapset()
m.glist('rast', pattern='elev-*')
['elev-002-000', 'elev-000-000', 'elev-000-002', 'elev-000-001',
'elev-001-001', 'elev-002-002', 'elev-002-001', 'elev-001-000',
'elev-001-002']
</source>

... that uses the "G_list" function through ctypes and therefore is faster.

Note: If you choose the first solution it is better if you use directly the
function in [http://grass.osgeo.org/programming7/namespacepython_1_1script_1_1core.html python.script.core.mlist_grouped() etc.].

=== Interface to copying maps (g.copy) ===
<source lang="python">
from grass.pygrass.modules.shortcuts import general as g

vectinmap = 'zipcodes_wake'
vectoutmap = 'zip_elev_zonal'
g.copy(vect="%s,%s" % (vectinmap,vectoutmap) )
</source>

=== Sample script for opening, query and closing of a vector map ===

<source lang="python">
#!/usr/bin/env python

# Example for pyGRASS usage - vector API

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.vector import VectorTopo

g.message("Assessing vector topology...")

vectmap = 'zipcodes_wake'
zipcodes = VectorTopo(vectmap)

# Open the map with topology:
zipcodes.open()

# query number of topological features
areas = zipcodes.number_of("areas")
islands = zipcodes.number_of("islands")
print 'Map: <' + vectmap + '> with %d areas and %d islands' % (areas, islands)

# http://www.ing.unitn.it/~zambelli/projects/pygrass/attributes.html
# (note that above documentation is slightly outdated)
dblink = zipcodes.dblinks[0]
print 'DB name:'
print dblink.database
table = dblink.table()
print 'Column names:'
print table.columns.names()
print 'Column types:'
print table.columns.types()

zipcodes.close()
</source>

=== Sample script for Landsat 7 processing ===

<source lang="python">
#!/usr/bin/env python
#Date: 7th February, 2013
#Public domain, GRASS GIS

#Run this script in the terminal as:
#python python-grass.py
#For debugging, run as:
#python -i python-grass.py
#PURPOSE
#This script processes LANDSAT 7 ETM+ images
#1 - unzip *.gz files
#2 - import files in GRASS GIS Location of your choice (r.in.gdal)
#3 - DN to Top of Atmosphere reflectance (i.landsat.toar)
#4 - TOA reflectance to Surface reflectance (i.atcorr)
#5 - NDVI (i.vi), Albedo (i.albedo), Emissivity (i.emissivity)

#USER HAS TO SET THOSE
#QUIET REPORTS
QIET=True
#OVERWRITE EXISTING FILES
OVR=False
#Define Landsat 7 sensor for i.landsat.toar
LSENSOR="tm7"
#Setup the path to the Landsat 7 Directories
rsdatapath="/home/yann/RSDATA/Myanmar/L7"
#Setup your GRASS GIS working directory
gisdb="/home/yann/GRASSDATA"
location="L7_Myanmar"
print location
mapset="PERMANENT"
#DEM input to atmospheric correction
inDEM=rsdatapath+"/dem.tif"
#set L7 Metadata wildcards
wldc_mtl="*_MTL.txt"
#wldc_met="*.met"
#Visibility distance [Km]
vis=18

#Set python path to enable finding of grass.script lib

#END OF USER CHANGES
###DO NOT CHANGE ANYTHING AFTER THIS LINE !
###
#Load necessary libraries
import os, glob, time, re
from grass import script as g
from grass.script import setup as gsetup
gisbase=os.environ['GISBASE']
gsetup.init(gisbase,gisdb,location,mapset)
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules.shortcuts import display as d
#Needed for floor()
from math import *
#env = os.environ.copy()
#env['GRASS_MESSAGE_FORMAT'] = 'gui'
#Function to get a list of L7 Directories in the rsdatapath
def fn(path):
for top, dirs, files in os.walk(path):
return [os.path.join(top, dir) for dir in dirs]

#START PROCESS
### PART 0: PRE-PROCESSING STUFF ###
#import DEM for atmospheric correction
#r.in.gdal(input=inDEM,output="dem",overwrite=OVR)
#r.mapcalc(expression="dem=25",overwrite=OVR)
#create a visibility map
r.mapcalc(expression="vis=18",overwrite=OVR)
#Find the central location of the Landsat file from metadata
metadata=[]
fileList=[]
L7Dirs=fn(rsdatapath)
for L7Dir in L7Dirs:
#Ungzip all of your Landsat7 images in all your directories
# print "Ungzip Landsat files in\t",L7Dir
# p=os.system("gzip -d -q "+L7Dir+"/*.gz")
#Using pthreads on multi-core machines
#p=os.system("pigz -d "+L7Dir+"/*.gz")
#Wait ten seconds for gzip to create the tif images
# time.sleep(10)
print "Import in GRASS GIS"
for L7f in glob.glob(os.path.join(L7Dir,"*.[tT][iI][fF]")):
f1=L7f.replace(L7Dir+"/","")
f2=f1.replace(".TIF","")
f3=f2.replace("_B10",".1")
f4=f3.replace("_B20",".2")
f5=f4.replace("_B30",".3")
f6=f5.replace("_B40",".4")
f7=f6.replace("_B50",".5")
f8=f7.replace("_B61",".61")
f9=f8.replace("_B62",".62")
f10=f9.replace("_B70",".7")
f11=f10.replace("_B80",".8")
f12=f11.replace("L72","L71")
L7r=f12.replace("_VCID_","")
print "\t> ",L7r
r.in_gdal(input=L7f,output=L7r,flags="e",overwrite=OVR)
fileList.append(L7r)

#reproject the DEM World map for the new PERMANENT location extents
r.proj(input="dem",location="Myanmar",memory=10000,resolution=90.0,overwrite=OVR)
#Get list of metadata files
for metaf in glob.glob(os.path.join(L7Dir,wldc_mtl)):
metadata.append(metaf)
print "Metadata in:\n",metadata[0]
with open(metadata[0],"r") as f:
data=f.read()

f.close()
dt=data.split("\n")
for idx in range(len(dt)):
found=dt[idx].find("CORNER_UL_LON_PRODUCT = ")
if found > 0:
ulx=float(dt[idx].replace("CORNER_UL_LON_PRODUCT = ",""))
found=dt[idx].find("CORNER_UL_LAT_PRODUCT = ")
if found > 0:
uly=float(dt[idx].replace("CORNER_UL_LAT_PRODUCT = ",""))
found=dt[idx].find("CORNER_LR_LON_PRODUCT = ")
if found > 0:
lrx=float(dt[idx].replace("CORNER_LR_LON_PRODUCT = ",""))
found=dt[idx].find("CORNER_LR_LAT_PRODUCT = ")
if found > 0:
lry=float(dt[idx].replace("CORNER_LR_LAT_PRODUCT = ",""))
#Two cases, because a string starting by 0 is not converted well
found=dt[idx].find("SCENE_CENTER_TIME = ")
if found > 0:
tim=dt[idx].replace("SCENE_CENTER_TIME = ","")
print "timestamp=",str(tim)
gmth=float(tim.split(":")[0])
gmtm=float(tim.split(":")[1])
gmtdec=float(tim.split(":")[2].replace("Z",""))
print gmth,gmtm,gmtdec
found=dt[idx].find("DATE_ACQUIRED = ")
if found > 0:
dat=dt[idx].replace("DATE_ACQUIRED = ","")
date=dat.split("-")
found=dt[idx].find("SUN_ELEVATION = ")
if found > 0:
sunza=90-float(dt[idx].replace("SUN_ELEVATION = ",""))

#L7 DN->Rad->Ref
print "Convert DN to Rad to TOARef"
f1=sorted(fileList)[0]
f2=f1.replace(L7Dir+"/","")
pref=f2[:-1]
outpref=pref[:-2]+".toar."
print "Prefix IN:\t",pref
print "Prefix OUT:\t", outpref
i.landsat_toar(input_prefix=pref,output_prefix=outpref,metfile=metadata[0],sensor=LSENSOR,quiet=QIET,overwrite=OVR)
#Atmospheric Correction
print "Atmospheric Correction"
# Basic script for i.atcorr for L 7 ETM+
#Geometrical conditions (L7ETM+)
geom=7
#Sensor height (satellite is -1000)
sens_height=-1000
#Here we suppose you have altitude (DEM) and Visibility (VIS) maps ready
#---------------------------------------------
#Visibility dummy value (overwritten by VIS raster input)
vis=15
#Altitude dummy value (in Km should be negative in this param file)
#(overwritten by DEM raster input)
alt=-1.200
#datetime of satellite overpass (month, day, GMT decimal hour)
mdh=str(int(date[1]))+" "+str(int(date[2]))+" "+str("%.2f" % (gmth+gmtm/60.0+gmtdec/3600.0))
print "MM DD hh.ddd:\t",mdh
# Central Lat/Long
Long=ulx+(lrx-ulx)/2.0
Lat=lry+(uly-lry)/2.0
print "Center:\t(",Long, ",", Lat,")"
#Atmospheric mode
atm_mode=1 #Tropical
#Aerosol model
aerosol_mode=2 #Sri Lanka is a small island (maritime)
#satellite band number (L5TM [25,26,27,28,29,30], L7ETM+ [61,62,63,64,65,66,67])
satbandno=61 #Band 1 of L7ETM+ is first to undergo atmospheric correction
#make time stamp for use in time series analysis
dat1=str(date[2])+"-"+str(date[1])+"-"+str(date[0])
dat=dat1.replace(" ","")
timestamp=time.strftime("%d %b %Y",time.strptime(dat,"%d-%m-%Y"))
print "Timestamp:\t",timestamp
for idx in [1,2,3,4,5,7]:
b=pref[:-2]+".toar."+str(idx)
b_out=b.replace(".toar.",".surf.")
param=[]
param.append(str(geom)+" - geometrical conditions=Landsat 7 ETM+\n")
param.append(str(mdh)+" "+str("%.2f" % Long)+" "+str("%.2f" % Lat)+" - month day hh.ddd longitude latitude (hh.ddd is in decimal hours GMT)\n")
param.append(str(atm_mode)+" - atmospheric mode=tropical\n")
param.append(str(aerosol_mode)+" - aerosols model=maritime\n")
param.append(str(vis)+" - visibility [km] (aerosol model concentration), not used as there is raster input\n")
param.append(str(alt)+" - mean target elevation above sea level [km] (here 600m asl), not used as there is raster input\n")
param.append(str(sens_height)+" - sensor height (here, sensor on board a satellite)\n")
param.append(str(satbandno)+" - i th band of ETM+ Landsat 7 (atcorr internal no)\n")
f=open(os.path.join(L7Dir,"param_L7.txt"),"w")
f.writelines(param)
f.close()
prm=os.path.join(L7Dir,"param_L7.txt")
print "\t> ",b
print "\t> ",b_out
i.atcorr(input=b, elevation="dem", visibility="vis", parameters=prm, output=b_out, flags="ra", range=[0,1],quiet=QIET,overwrite=OVR)
r.timestamp(map=b_out,date=timestamp,finish_=False)
satbandno = satbandno + 1

#Allocate surface reflectance names
b1=pref[:-2]+".surf.1"
b2=pref[:-2]+".surf.2"
b3=pref[:-2]+".surf.3"
b4=pref[:-2]+".surf.4"
b5=pref[:-2]+".surf.5"
b61=pref[:-2]+".toar.61"
b62=pref[:-2]+".toar.62"
b7=pref[:-2]+".surf.7"
b8=pref[:-2]+".surf.8"

### PART 1: BASIC STUFF ###
b_in=pref[:-2]+".toar."
b_clouds=pref[:-2]+".toar.acca"
print "Clouds:\t\t>",b_clouds
i.landsat_acca(input_prefix=b_in,output=b_clouds,overwrite=OVR)
#png_clouds=L7Dir+"/"+pref[:-2]+".clouds.png"
#d.mon(start='png',output=png_clouds,width=800,height=800)
print "MASK:\t\tON"
#Should always be rewritten!
r.mask(raster=b_clouds,flags="i",overwrite=True)

b_ndvi=pref[:-2]+".surf.ndvi"
print "NDVI:\t",b_ndvi
i.vi(red=b3,nir=b4,output=b_ndvi,viname="ndvi",quiet=QIET,overwrite=OVR,finish_=False)

b_in=[b1,b2,b3,b4,b5,b7]
b_albedo=pref[:-2]+".surf.albedo"
print "Albedo:\t",b_albedo
i.albedo(input=b_in,output=b_albedo,flags="lc",quiet=QIET,overwrite=OVR,finish_=False)

b_emissivity=pref[:-2]+".surf.emissivity"
print "Emissivity:\t",b_emissivity
i.emissivity(input=b_ndvi, output=b_emissivity,quiet=QIET,overwrite=OVR,finish_=False)
</source>

=== Multiprocessing example: parallelized SHAPE file import ===

<source lang="python">
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Multiprocessing example: parallelized SHAPE file import

Created on Fri, Oct 18, 2013, posted to grass-user@lists.osgeo.org

@author: Pietro Zambelli, freely inspired by: http://stackoverflow.com/a/16071616
http://lists.osgeo.org/pipermail//grass-user/2013-October/069130.html
"""

# Directory containing SHAPE files to import
DIR = '/data/shp'

###
from multiprocessing import Queue, Process, cpu_count
from os.path import split
from subprocess import Popen

from grass.pygrass.functions import findfiles

def spawn(func):
def fun(q_in, q_out):
while True:
path, cmdstr = q_in.get()
if path is None:
break
q_out.put(func(path, cmdstr))
return fun

def mltp_importer(dirpath, match, cmdstr, func, nprocs=cpu_count()):
q_in = Queue(1)
q_out = Queue()
procs = [Process(target=spawn(func), args=(q_in, q_out))
for _ in range(nprocs)]
for proc in procs:
proc.daemon = True
proc.start()

# set the parameters
sent = [q_in.put((path, cmdstr)) for path in findfiles(dirpath, match)]
# set the end of the cycle
[q_in.put((None, None)) for proc in procs]
[proc.join() for proc in procs]
return [q_out.get() for _ in range(len(sent))]

def importer(path, cmdstr):
name = split(path)[-1][:-4]
print name
popen = Popen(cmdstr.format(path=path, name=name), shell=True)
popen.wait()
return path, name, False if popen.returncode else True

CMD = 'v.in.ogr dsn={path} layer={name} output={name}'

processed = mltp_importer(DIR, '*.shp', CMD, importer)
# check for errors
errors = [p for p in processed if not p[2]]
if errors:
# do something (print list of failed SHP files at end)
pass
</source>

{{Python}}

Creating watersheds

2015-12-07T21:13:28Z

⚠️Pmav99: /* Combining watersheds into one patched vector */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

And here's how it works for GRASS 7:
<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do
i=$(( ${i} + 1 ))
r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i
r.null tmp_$i setnull=0
echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=-
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
# r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 6
# r.to.vect in=wshed_patch out=wshed_patch feature=area # Grass 6
r.patch in=`g.list raster pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 7
r.to.vect in=wshed_patch out=wshed_patch type=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
# g.remove rast=`g.mlist pattern=tmp* sep=,` # Grass 6
g.remove -f type=raster name=`g.list raster pattern=tmp* sep=, ` # Grass 7
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
# v.db.addcol map=wshed_final col="area_sqkm double" # Grass 6
v.db.addcolumn map=wshed_final col="area_sqkm double" # Grass 7

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T21:11:14Z

⚠️Pmav99: /* Combining watersheds into one patched vector */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

And here's how it works for GRASS 7:
<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do
i=$(( ${i} + 1 ))
r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i
r.null tmp_$i setnull=0
echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=-
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
# r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 6
# r.to.vect in=wshed_patch out=wshed_patch feature=area # Grass 6
r.patch in=`g.list raster pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 7
r.to.vect in=wshed_patch out=wshed_patch type=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
# g.remove rast=`g.mlist pattern=tmp* sep=,` # Grass 6
g.remove -f type=raster name=`g.list raster pattern=tmp* sep=, ` # Grass 7
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T21:09:14Z

⚠️Pmav99: /* Combining watersheds into one patched vector */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

And here's how it works for GRASS 7:
<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do
i=$(( ${i} + 1 ))
r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i
r.null tmp_$i setnull=0
echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=-
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
# r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 6
# r.to.vect in=wshed_patch out=wshed_patch feature=area # Grass 6
r.patch in=`g.list raster pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 7
r.to.vect in=wshed_patch out=wshed_patch type=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
# g.remove rast=`g.mlist pattern=tmp* sep=,` # Grass 6
g.remove rast=`g.list raster pattern=tmp* sep=, ` # Grass 7
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T21:08:10Z

⚠️Pmav99: /* Combining watersheds into one patched vector */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

And here's how it works for GRASS 7:
<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do
i=$(( ${i} + 1 ))
r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i
r.null tmp_$i setnull=0
echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=-
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
# r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 6
# r.to.vect in=wshed_patch out=wshed_patch feature=area # Grass 6
r.patch in=`g.list raster pattern=tmp_reclass* separator=,` out=wshed_patch # Grass 7
r.to.vect in=wshed_patch out=wshed_patch type=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T21:01:44Z

⚠️Pmav99: /* Looping thru drainage outlet points */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

And here's how it works for GRASS 7:
<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do
i=$(( ${i} + 1 ))
r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i
r.null tmp_$i setnull=0
echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=-
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:55:29Z

⚠️Pmav99: /* Determine drainage points */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:52:19Z

⚠️Pmav99: /* Determine drainage points */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point separator=space # Grass 7
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:43:44Z

⚠️Pmav99: /* Display first results */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:43:29Z

⚠️Pmav99: /* Display first results */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area # Grass 6
r.to.vect in=catch out=catchments type=area # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# If you are using GRASS 6, then you should use the following line
# r.to.vect in=str_thin out=streams feature=line # Grass 6
r.to.vect in=str_thin out=streams type=line # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# r.shaded.relief map=dem shade=dem_shade zmult=1.5 # Grass 6
r.relief input=dem output=dem_shade zscale=1.5 # Grass 7

# Now display layers
# d.mon x0 # Grass 6
d.mon start=wx0 # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:39:23Z

⚠️Pmav99: /* Display first results */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
# If you are using GRASS 6, then you should use the following line
# r.to.vect in=catch out=catchments feature=area
r.to.vect in=catch out=catchments type=area

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# If you are using GRASS 6, then you should use the following line
# r.to.vect in=str_thin out=streams feature=line
r.to.vect in=str_thin out=streams type=line
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# If you are using GRASS 6, then you should use the following line
# r.shaded.relief map=dem shade=dem_shade zmult=1.5
r.relief input=dem output=dem_shade zscale=1.5

# Now display layers
d.mon x0
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]

Creating watersheds

2015-12-07T20:36:56Z

⚠️Pmav99: /* Display first results */

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

== Creating watersheds with specific drainage outlets==
As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

=== Preparation ===

<source lang="bash">
#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>
</source>

So far, pretty straightforward. There's abundant information on {{cmd|r.watershed}}. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

=== Display first results ===
When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

<source lang="bash">
#Convert the steams and catchments to vectors
r.to.vect in=catch out=catchments feature=area

# the stream raster usually requires thinning
r.thin in=str out=str_thin
r.to.vect in=str_thin out=streams feature=line
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D"
# If you are using GRASS 6, then uncomment the following line
#r.shaded.relief map=dem shade=dem_shade zmult=1.5
r.relief input=dem output=dem_shade zscale=1.5

# Now display layers
d.mon x0
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2
</source>

=== Determine drainage points ===
Now we need to find all the points where streams cross roads. The {{cmd|v.overlay}} module does not deal with point vectors (hint: {{cmd|v.select}} does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run {{cmd|v.clean}}, and we get all those intersection points in a new vector.

<source lang="bash">
# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space
</source>

=== Looping thru drainage outlet points ===
Once we have the crossing points in a file, we simply run {{cmd|r.water.outlet}} in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have '''null values''' outside of the watersheds, and each watershed must use a '''different value''' in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the {{cmd|r.reclass}} function to make a reclassed raster with different values for each watershed. Here's how it works:

<source lang="bash">
i=0 #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 )) \
r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i \
r.null tmp_$i setnull=0 \
echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"
</source>

=== Combining watersheds into one patched vector ===
Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

<source lang="bash">
r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch
r.to.vect in=wshed_patch out=wshed_patch feature=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150
</source>

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

<source lang="bash">
g.remove rast=`g.mlist pattern=tmp* sep=,`
</source>

Most likely we'll want to calculate the area for each watershed.

<source lang="bash">
# Add a table with a column for area in sq.km.
v.db.addcol map=wshed_final col="area_sqkm double"

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k
</source>

And finally, we can view the catchments, and their area values (you may use the [[wxGUI]]):

<source lang="bash">
d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange
</source>

[[Category: FAQ]]
[[Category: Documentation]]
[[Category: Hydrology]]