import datetime as dt
start = dt.datetime.now()
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
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'
from IPython.display import display, Math, Latex, Javascript
from IPython.core.display import Image
import grass.script as grass
import pandas as pd
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
GISDBASE = '/home/epifanio/dev/src/wython/grass7data'
!g.gisenv set="GISDBASE={GISDBASE}"
!g.gisenv set="LOCATION_NAME=nc_spm_08_grass7"
!g.gisenv set="MAPSET=PERMANENT"
!g.mapset mapset=PERMANENT location=nc_spm_08_grass7
WARNING: <PERMANENT> is already the current mapset
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.
New hydrologic modules have been developed recentely making GRASS a powerful toolset for hydrological analysis even for large datasets.
Inputs needed to run r.basin are :
Powered by r.stream.* hydrological tools, it performs:
Based on those maps, it calculates :
!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
Calculate flow accumulation map (MFD) - r.watershed
Extract the stream network - r.stream.extract
!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
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).
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.
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.
prefix='outx'
threshold=20
easting=636654.791181
northing=218824.126649
!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!
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.
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.
pd.read_csv('%s_elevation_parameters.csv' % prefix, skiprows=1)
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 |
makeImage(mapname='%s_elevation_accumulation' % prefix, maptitle='Accumulation')
Image(filename="%s_elevation_accumulation.png" % prefix)
makeImage(mapname='%s_elevation_aspect' % prefix, maptitle='Aspect')
Image(filename="%s_elevation_aspect.png" % prefix)
makeImage(mapname='%s_elevation_slope' % prefix, maptitle='Slope')
Image(filename="%s_elevation_slope.png" % prefix)
makeImage(mapname='%s_elevation_dist2out' % prefix, maptitle='Dist2out')
Image(filename="%s_elevation_dist2out.png" % prefix)
makeImage(mapname='%s_elevation_drainage' % prefix, maptitle='Drainage')
Image(filename="%s_elevation_drainage.png" % prefix)
makeImage(mapname='%s_elevation_hack' % prefix, maptitle='Hack')
Image(filename="%s_elevation_hack.png" % prefix)
makeImage(mapname='%s_elevation_hillslope_distance' % prefix, maptitle='Hillslopes')
Image(filename="%s_elevation_hillslope_distance.png" % prefix)
makeImage(mapname='%s_elevation_horton' % prefix, maptitle='Horton')
Image(filename="%s_elevation_horton.png" % prefix)
makeImage(mapname='%s_elevation_shreve' % prefix, maptitle='Shreve')
Forcing a smooth legend: too many categories for current window height
Image(filename="%s_elevation_shreve.png" % prefix)
makeImage(mapname='%s_elevation_strahler' % prefix, maptitle='Strahler')
Image(filename="%s_elevation_strahler.png" % prefix)
lcolor='blue'
pcolor='red'
acolor='red'
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)
mapname='catchment'
elev='elevation'
basin='%s_elevation_basin' % prefix
mainchannel='%s_elevation_mainchannel' % prefix
outlet='%s_elevation_outlet' % prefix
network="%s_elevation_network" % prefix
!rm -rf {mapname}.png
!d.mon start=cairo output={mapname}.png
!g.region rast={elev} n={n} s={s} w={w} e={e} -a #p
!d.rast map={elev} --q
!d.legend at={legend_at[0]},{legend_at[1]},{legend_at[2]},{legend_at[3]} map={elev}
!d.vect map={basin} color={acolor} type=boundary width=3 --q
!d.vect map={mainchannel} color={lcolor} width=3 --q
!d.vect map={network} color={lcolor} width=1 --q
!d.vect map={mainchannel} color={lcolor} width=3 --q
!d.vect map={outlet} color={pcolor} size=18 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
Image(filename="catchment.png")
The hypsographic curve provides the distribution of the areas at different altitudes. Each point on the hypsographic curve has on the y-axis the altitude and on the x-axis the percentage of basin surface placed above that altitude.
Image(filename='%s_elevation_Ipsographic.png' % prefix)
The hypsometric curve has the same shape but is dimensionless. Quantiles of the hypsometric curve are displayed:
Image(filename='%s_elevation_Ipsometric.png' % prefix)
The width function W(x) gives the area of the cells in the basin at a certain flow distance x from the outlet (distance-area function). Note that the distance is not intended in the euclidean sense, but it's calculated considering the hydrological path of the water.
Image(filename='%s_elevation_width_function.png' % prefix)
end = dt.datetime.now()
time = end-start
print 'processing performed in : ', time
processing performed in : 0:01:11.785753
display(Javascript("IPython.notebook.save_notebook()"))
The advantages offered by the implementation in Python and GRASS are manifold.
Python is a powerful scripting language with huge potential for researchers due to its relative simplicity, high flexibility and thanks to a broad availability of scientific libraries.
GRASS, and as a consequence, r.basin, is platform independent, so that it is available for GNU/Linux, MS Windows, Mac, etc.
Furthermore, the module is constantly maintained and improved according to users' feedback with the precious help of expert developers.
The code is available for review under the official GRASS add-ons repository, allowing hydrologists and researchers to knowingly use, inspect, modify, reuse, and even incorporate it in other projects, such as web services.
#lines = !jist -p AGU-2013-H52E02-MDS.ipynb
#print lines[0].replace("https://gist.github.com", "http://nbviewer.ipython.org")
#print lines[0].replace("https://gist.github.com", "https://gist.github.com/anonymous")
#!ipython nbconvert AGU-2013-H52E02-MDS.ipynb --to slides
#!ipython nbconvert --to html --template full AGU-2013-H52E02-MDS.ipynb
ipython nbconvert AGU-2013-H52E02-MDS.ipynb --to slides
ipython nbconvert --to html --template full AGU-2013-H52E02-MDS.ipynb
!jist -p AGU-2013-H52E02-MDS_v3.ipynb
https://gist.github.com/8426314