An Open-Source Approach for Catchment's Physiographic Characterization


Massimo Di Stefano - Rensselaer Polytechnic Institute (RPI), NY, USA

Margherita Di Leo - European Commission, Joint Research Centre, Italy

In [1]:
import datetime as dt
start = dt.datetime.now()

Enable access to GRASS GIS executables and configure GRASS Environment

In [2]:
import sys
sys.path.append('/home/epifanio/dev/src/wython/grass-7.0.svn/etc/python')
sys.path.append('/home/epifanio/dev/src/wython/bin')
sys.path.append('/home/epifanio/dev/src/wython/grass-7.0.svn/scripts')
sys.path.append('/home/epifanio/dev/src/wython/grass-7.0.svn/bin')
import grass
In [3]:
import os
os.environ['GISBASE'] = '/home/epifanio/dev/src/wython/grass-7.0.svn'
os.environ['GIS_LOCK'] = '$$'
os.environ['GISRC'] = '/home/epifanio/.grass7/rc'
os.environ['GISDBASE'] = '/home/epifanio/dev/src/wython/grass7data'
os.environ['GRASS_TRANSPARENT'] = 'TRUE'
os.environ['GRASS_TRUECOLOR'] = 'TRUE'
os.environ['GRASS_PNG_COMPRESSION'] = '9'
os.environ['GRASS_PNG_AUTO_WRITE'] = 'TRUE'
os.environ['PATH'] = '/home/epifanio/.gem/ruby/1.8/bin:/home/epifanio/dev/src/wython/grass-7.0.svn/bin:/home/epifanio/dev/src/wython/grass-7.0.svn/scripts:/usr/bin:/bin'
os.environ['LD_LIBRARY_PATH']='/home/epifanio/dev/src/wython/grass-7.0.svn/lib'

Import Python Modules

In [4]:
from IPython.display import display, Math, Latex, Javascript
from IPython.core.display import Image
import grass.script as grass
import pandas as pd

Define a function as 'd.mon'wrapper to render GRASS layers as PNG

In [5]:
def makeImage(mapname='', maptype='raster', lcolor='grey', vcolor='grey', vsize=1, icon='basic/diamond', maptitle='', n=223000, s=218000, w=632000, e=638000, gridsize=1000, title_at=(45,95), legend_at=(22,65,8,10)):
    !rm -rf {mapname}.png
    !d.mon start=cairo output={mapname}.png
    if maptype=='raster':
        !g.region rast={mapname} n={n} s={s} w={w} e={e} -a #p
        !d.rast map={mapname} --q
        !d.legend at={legend_at[0]},{legend_at[1]},{legend_at[2]},{legend_at[3]} map={mapname} 
    if maptype=='vector':
        !g.region vect={mapname} n={n} s={s} w={w} e={e} -a #p
        !d.vect map={mapname} color={vcolor} size={vsize} icon={icon} --q
    !d.grid size={gridsize} color=blue
    !d.text at={title_at[0]},{title_at[1]} text={maptitle} color=red
    !d.mon stop=cairo

Set the GRASS Working Environment

In [6]:
GISDBASE = '/home/epifanio/dev/src/wython/grass7data'
In [7]:
!g.gisenv set="GISDBASE={GISDBASE}"
!g.gisenv set="LOCATION_NAME=nc_spm_08_grass7"
!g.gisenv set="MAPSET=PERMANENT"
In [8]:
!g.mapset mapset=PERMANENT location=nc_spm_08_grass7
WARNING: <PERMANENT> is already the current mapset


Introduction


  • Extraction of physiographic descriptors of catchments is a fundamental element of hydrological studies.

  • They have been extracted manually for a long time, but currently Geographic Information System (GIS) tools ease such task and provide better results, offering a powerful instrument for hydrologists.

  • Also, extracting these descriptors can be a time-consuming task, which in turn needs to be repeated for every new analysis. For this type of challenge extending GIS capabilities through the use of scripting language can be extremely helpful.

Combining the flexibility of the Python programming language with the reliability of GRASS GIS, here we present a tool called r.basin, which automatically performs physiographic characterization of catchments.

  • GRASS GIS (Geographic Resource Analysis Support System) is a Free and Open Source tool for geospatial data management and analysis, image processing, graphics and maps production, spatial modeling and visualization. New hydrologic modules have been developed recentely making GRASS a powerful toolset for hydrological analysis even for large datasets.

  • Python is a scripting language widely used in computational science characterized by very clean syntax, rich modularization features, good support for numerical computing, and rapidly growing popularity.

Inputs needed to run r.basin are :

  • Digital Elevation Model
  • Coordinates of the outlet

Powered by r.stream.* hydrological tools, it performs:

  • Flow calculation
  • Basin boundaries
  • Drainage network
  • Flow direction and accumulation
  • Distance to outlet
  • Hill slope length maps.

Based on those maps, it calculates :

  • Hydrologically meaningful shape factors
  • Topological diameter
  • Drainage density
  • Horton's ratios
  • Concentration time
  • Statistics on main channel and elevation
  • Centroid's coordinates
  • Rectangle containing the basin, etc.
  • Hypsographic and hypsometric curve
  • Width Function

Preparation


As a first step, we set the computational region to match the elevation raster map:

In [9]:
!g.region rast=elevation@PERMANENT -ap
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      10
ewres:      10
rows:       1350
cols:       1500
cells:      2025000

Stream network extraction


In [10]:
!r.watershed -a elevation=elevation@PERMANENT accumulation=accum --o --q
!r.stream.extract elevation=elevation@PERMANENT accumulation=accum threshold=20 stream_rast=stream_network stream_vect=streams --o --q
WARNING: Writing out only positive flow accumulation values.
WARNING: Cells with a likely underestimate for flow accumulation can no
         longer be identified.
WARNING: Vector map <streams> already exists and will be overwritten


r.basin Usage


Parameters :

  • Prefix : string given by the user to distinguish the maps produced by every run of the program

  • Easting Northing1 : a pair of coordinates for the outlet belonging to the calculated stream network

  • Threshold2 : number of cells that determine where the river begins (same as in r.watershed and r.stream.extract).


  1. For the basin's delineation, a pair of coordinates is required. Usually coordinates belonging to the natural river network don't exactly match with the calculated stream network. What we should do is to calculate the stream network first, and then to find the coordinates on the calculated stream network closest to the coordinates belonging to the natural stream network.

  2. Threshold value is usually determined by trial and error. r.basin has a flag -a, which uses an automatic threshold value of 1km2 (suitable for Italy). For more detailed information on how to determine the threshold, see also r.threshold and r.stream.preview.

Run IT!

In [11]:
prefix='outx'
threshold=20
easting=636654.791181 
northing=218824.126649
In [12]:
!r.basin map=elevation@PERMANENT prefix={prefix} easting={easting} northing={northing} threshold={threshold}
WARNING: Writing out only positive flow accumulation values.
WARNING: Cells with a likely underestimate for flow accumulation can no
         longer be identified.
Running r.-stream.extract
running r.stream.basin
Delineation of basin done
WARNING: Vector map <outx_elevation_network> already exists and will be
         overwritten
WARNING: Vector map <outx_elevation_basin> already exists and will be
         overwritten
Reading areas...
 100%
Reading areas...
 100%
Creating outx_elevation_hack
WARNING: Vector map <outx_elevation_outlet> already exists and will be
         overwritten
##################################
Tot. cells 78593.0
===========================
Ipsometric | quantiles
===========================
145 | 0.025
143 | 0.05
141 | 0.1
137 | 0.25
129 | 0.5
119 | 0.75
122 | 0.7
110 | 0.9
100 | 0.975
##################################
Tot. cells 78593.0
Tot. area 7859300.0
Max distance 6769.330609
===========================
Whidth Function | quantiles
===========================
802 | 0.05
1508 | 0.15
2168 | 0.3
2573 | 0.4
2960 | 0.5
3632 | 0.6
4186 | 0.7
5131 | 0.85
6089 | 0.95
##################################
r.statistics2 done
r.info done
r.mapcalc done
r.volume done
g.region done
Rectangle containing basin done
Directing vector done
Prevalent orientation done
Compactness coefficient done
Circularity ratio done
doing r.thin
WARNING: Vector map <outx_elevation_mainchannel> already exists and will be
         overwritten
doing v.what
doing v.to.points
##################################
Output of r.stream.stats: 

Summary:
Max order | Tot.N.str. | Tot.str.len. | Tot.area. | Dr.dens. | Str.freq. 
  (num)   |    (num)   |     (km)     |   (km2)   | (km/km2) | (num/km2) 
        5 |        612 |      80.8272 |    7.8593 |  10.2843 | 77.8695 

Stream ratios based on regresion coefficient:
 Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. 
  5.0940 |  2.1837 |   5.6571 |  1.4418 |  1.5420

Avaraged stream ratios with standard deviations:
 Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. 
  5.5417 |  2.6271 |   5.1916 |  1.4625 |  1.6693
  3.5678 |  1.9342 |   4.3638 |  0.1448 |  0.6477

Order | Avg.len |  Avg.ar  |  Avg.sl |  Avg.grad. | Avg.el.dif
 num  |   (km)  |  (km2)   |  (m/m)  |    (m/m)   |     (m)   
    1 |  0.0880 |   0.0100 |  0.0428 |     0.0331 |  3.1995
    2 |  0.2025 |   0.0501 |  0.0303 |     0.0250 |  5.3206
    3 |  0.5337 |   0.2538 |  0.0212 |     0.0177 |  8.5684
    4 |  2.7421 |   2.7104 |  0.0159 |     0.0136 | 25.1265
    5 |  1.1871 |   7.8593 |  0.0095 |     0.0051 |  6.1054

Order | Std.len |  Std.ar  |  Std.sl |  Std.grad. | Std.el.dif
 num  |   (km)  |  (km2)   |  (m/m)  |    (m/m)   |     (m)   
    1 |  0.0771 |   0.0082 |  0.0362 |     0.0275 |  3.5790
    2 |  0.1841 |   0.0398 |  0.0185 |     0.0145 |  5.6062
    3 |  0.6182 |   0.2493 |  0.0131 |     0.0120 |  8.5329
    4 |  2.8546 |   3.1286 |  0.0062 |     0.0085 | 15.5228
    5 |  0.0000 |   0.0000 |  0.0000 |     0.0000 |  0.0000

Order | N.streams | Tot.len (km) | Tot.area (km2)
    1 |       490 |      43.1013 |  4.8895
    2 |        98 |      19.8470 |  4.9085
    3 |        21 |      11.2077 |  5.3290
    4 |         2 |       5.4841 |  5.4207
    5 |         1 |       1.1871 |  7.8593

Order | Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. | d.dens. | str.freq.
    1 |  5.0000 |  2.3024 |   0.0000 |  1.4154 |  1.3236 |  8.8151 | 100.2147
    2 |  4.6667 |  2.6353 |   5.0194 |  1.4279 |  1.4086 |  4.0434 | 19.9654
    3 | 10.5000 |  5.1378 |   5.0664 |  1.3358 |  1.3065 |  2.1031 |  3.9407
    4 |  2.0000 |  0.4329 |  10.6807 |  1.6710 |  2.6385 |  1.0117 |  0.3690
    5 |  0.0000 |  0.0000 |   2.8997 |  0.0000 |  0.0000 |  0.1510 |  0.1272


##################################
Morphometric parameters of basin :
##################################
Easting Centroid of basin : 268000
Northing Centroid of Basin : 76477
Rectangle containing basin N-W : 632880 , 222890
Rectangle containing basin S-E : 637010 , 218720
Area of basin [km^2] : 7.8592625
Perimeter of basin [km] : 17.7583322395
Max Elevation [m s.l.m.] : 151.7396
Min Elevation [m s.l.m.]: 94.82206
Elevation Difference [m]: 56.91754
Mean Elevation [m s.l.m.]: 127.7042
Mean Slope : 3.50
Length of Directing Vector [km] : 395.182311757
Prevalent Orientation [degree from north, counterclockwise] :
0.368488942888
Compactness Coefficient : 5.61379074702
Circularity Ratio : 0.313175157621
Topological Diameter : 1.0
Elongation Ratio : 3163.34060883
Shape Factor : 7859.2625
Concentration Time (Giandotti, 1934) [hr] : 1.85821486566
Length of Mainchannel [km] : 0.001
Mean slope of mainchannel [percent] : 1.043132
Mean hillslope length [m] : 5932.153
Magnitudo : 415.0
Max order (Strahler) : 5
Number of streams : 612
Total Stream Length [km] : 80.8272
First order stream frequency : 52.8039367562
Drainage Density [km/km^2] : 10.2843237518
Bifurcation Ratio (Horton) : 5.0940
Length Ratio (Horton) : 2.1837
Area ratio (Horton) : 5.6571
Slope ratio (Horton): 1.4418
##################################

Done!

r.basin outputs :

  • Several morphometric parameters, which are printed in the terminal and also stored in a csv file called out_elevation_parameters.csv, in the working directory;

  • Charts;

  • Maps.

Basin morphometric parameters :

  • Bounding box: coordinates of the vertices of the rectangle containing the basin.
  • Center of gravity: coordinates of the pixel nearest to the center of gravity of the geometric figure resulting from the projection of the basin on the horizontal plane.
  • Area: area of a single cell multiplied by the number of cells belonging to the basin.
  • Perimeter: length of the contour of the figure resulting from the projection of the basin on the horizontal plane.
  • Elevation descriptors: max, min, range and mean elevation values
  • Mean slope: slope average
  • Length of the directing vector: length of the vector linking the outlet to basin barycenter.
  • Prevalent orientation: orientation of the directing vector.
  • Length of main channel: length of the longest succession of segments that connect a source to the outlet of the basin.
  • Mean slope of main channel: it is calculated as follows: \[S_{mean} = \frac{1}{N} \sum \frac{\Delta z_{i}}{L_{i}} \cdot 100\] where N is the topological diameter, i.e. the number of links in which the main channel can be divided on the basis of the junctions.
  • Circularity ratio : is the ratio between the area of the basin and the area of the circle having the same perimeter of the basin.
  • Elongation ratio : is the ratio between the diameter of the circle having the same area of the basin and the length of the main channel.
  • Compactness coefficient : is the ratio between the perimeter of the basin and the diameter of the circle having the same area of the basin.
  • Shape factor : is the ratio between the area of the basin and the square of the length of the main channel.
  • Concentration time (Giandotti, 1934) : \[t_c = \frac{4 \cdot \sqrt[]{A} + 1,5 L }{0.8 \cdot \sqrt[]{{H}}}\] Where A is the area, L the length of the main channel and H the difference between the highest and the lowest elevation of the basin.
  • Mean hillslope length : is the mean of the distances calculated along the flow direction of each point non belonging to the river network from the point in which flows into the network.
  • Magnitudo : is the number of the branches of order 1 following the Strahler hierarchy.
  • Max order : is the order of the basin, following the Strahler hierarchy.
  • Number of streams : is the number of the branches of the river network.
  • Total stream length : is the sum of the length of every branches.
  • First order stream frequency : is the ratio between the magnitudo and the area of the basin.
  • Drainage density : is the ratio between the total length of the river network and the area.
  • Horton ratios (Horton, 1945; Strahler, 1957).
In [13]:
 pd.read_csv('%s_elevation_parameters.csv'  % prefix, skiprows=1)
Out[13]:
Easting Centroid of basin 268000
Northing Centroid of basin 76477
Rectangle containing basin N-W ('632880', '222890')
Rectangle containing basin S-E ('637010', '218720')
Area of basin [km^2] 7.8592625
Perimeter of basin [km] 17.7583322395147
Max Elevation [m s.l.m.] 151.7396
Min Elevation [m s.l.m.] 94.82206
Elevation Difference [m] 56.91754
Mean Elevation 127.7042
Mean Slope 3.50
Length of Directing Vector [km] 395.18231175741295
Prevalent Orientation [degree from north, counterclockwise] 0.3684889428877813
Compactness Coefficient 5.613790747022824
Circularity Ratio 0.31317515762091674
Topological Diameter 1.0
Elongation Ratio 3163.3406088270253
Shape Factor 7859.2625
Concentration Time (Giandotti, 1934) [hr] 1.858214865664734
Length of Mainchannel [km] 0.001
Mean slope of mainchannel [percent] 1.0431320016345824
Mean hillslope length [m] 5932.153
Magnitudo 415.0
Max order (Strahler) 5
Number of streams 612
Total Stream Length [km] 80.8272
First order stream frequency 52.803936756152375
Drainage Density [km/km^2] 10.284323751751517
Bifurcation Ratio (Horton) 5.0940
Length Ratio (Horton) 2.1837
Area ratio (Horton) 5.6571
Slope ratio (Horton) 1.4418

Raster maps :

- Flow accumulation

In [14]:
makeImage(mapname='%s_elevation_accumulation' % prefix, maptitle='Accumulation')
In [15]:
Image(filename="%s_elevation_accumulation.png" % prefix)
Out[15]:

- Aspect

In [16]:
makeImage(mapname='%s_elevation_aspect' % prefix, maptitle='Aspect')
In [17]:
Image(filename="%s_elevation_aspect.png" % prefix)
Out[17]:

- Slope

In [18]:
makeImage(mapname='%s_elevation_slope' % prefix, maptitle='Slope')
In [19]:
Image(filename="%s_elevation_slope.png" % prefix)
Out[19]:

- Distance to outlet (width function)

In [20]:
makeImage(mapname='%s_elevation_dist2out' % prefix, maptitle='Dist2out')
In [21]:
Image(filename="%s_elevation_dist2out.png" % prefix)
Out[21]:

- Drainage

In [22]:
makeImage(mapname='%s_elevation_drainage' % prefix, maptitle='Drainage')
In [23]:
Image(filename="%s_elevation_drainage.png" % prefix)
Out[23]:

- Main (Hack) Stream Order

In [24]:
makeImage(mapname='%s_elevation_hack' % prefix, maptitle='Hack')
In [25]:
Image(filename="%s_elevation_hack.png" % prefix)
Out[25]:

- Length of hillslopes

In [26]:
makeImage(mapname='%s_elevation_hillslope_distance' % prefix, maptitle='Hillslopes')
In [27]:
Image(filename="%s_elevation_hillslope_distance.png" % prefix)
Out[27]: