------------------------------------------------------------------ Automatic delineation of drainage basins from contour elevation data using skeleton construction techniques ------------------------------------------------------------------ July 29, 2009 Stefano Orlandini Dipartimento di Ingegneria Meccanica e Civile Università degli Studi di Modena e Reggio Emilia Via Vignolese 905, I-41125 Modena, Italy Work phone: +39 059 2056105 Fax: +39 05902056126 E-mail: stefano.orlandini@unimore.it URL: http://www.idrologia.unimore.it/orlandini This "readme" file contains the following sections: 1. The Contour Line Analysis Method 2. Distribution and Copyright 3. The Archive cla.tar.gz 4. How to Use the Fortran Codes That Implement the Contour Line Analysis 5. Examples 6. References 1. The CLA Method -------------------- The method is described in Moretti and Orlandini (2008). The .pdf file of this paper (file 2007wr006309.pdf) and the Fortran codes that implement the method are provided in the archive cla.tar.gz described below. 2. Distribution and Copyright ----------------------------- Copyright (C) 2003 Stefano Orlandini This program is free software. You can redistribute it and/or modify it under the terms of the GNU General Public License version 2, 1991 as published by the Free Software Foundation. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY. Without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the full GNU General Public License is available at the address: http://www.gnu.org/copyleft/gpl.html, or from: The Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. If you wish to use or incorporate this program (or parts of it) into other software that does not meet the GNU General Public License conditions contact the author to request permission. Stefano Orlandini Dipartimento di Ingegneria Meccanica e Civile Università degli Studi di Modena e Reggio Emilia Via Vignolese 905, I-41125 Modena, Italy Work phone: +39 059 2056105 Fax: +39 05902056126 E-mail: stefano.orlandini@unimore.it URL: http://www.idrologia.unimore.it/orlandini 3. The Archive cla.tar.gz ---------------------------- All the material mentioned in this "readme" file is contained in the tar format archive compressed with gzip "cla.tar.gz". The gzip and tar programs can be obtained as free software on the Web for any operating system. To extract the "cla" directory from "cla.tar.gz" use the command lines: gzip -d cla.tar.gz tar xf cla.tar Alternatively, with Windows systems, one can use the Winzip program. Fortran codes can be found in "cla". Data related to a real case application can be found in "cla/examples". 4. How to Use the Fortran Codes That Implement the Contour Line Analysis ------------------------------------------------------------------------ The material provided in this package includes: A) Five main programs written in Fortran 90 and related modules/subroutines. B) The input files for the automatic delineation of drainage basin (cla/examples/basin_delineation) and the construction of flow net (cla/examples/flow_net). C) Examples of the output files that one may obtain from the Fortran programs can be visualized in a GIS environment (e.g., ArcView 3.2/8.x/9.x). The steps required to run the examples are finally summarized. Fortran Programs ---------------- Five main programs written in Fortran 90 are provided: 1) shp2dat_node2 (file shp2dat_node2.f90) 2) crust_skel_vrnode (file crust_skel_vrnode.f90) 3) orgskel_shp (file orgskel_shp.f90) 4) bbt_plus (file bbt_plus.f90) 5) cb_dem_eqraf_new(file cb_dem_eqraf_new.f90) These programs compile successfully with the following Fortran 90 compilers: Lahey/Fujitsu Fortran 95 Compiler Release 5.70c The related command lines are provided in the file "cla/compf90". Modules and subroutines associated to each program are those used in the file "cla/compf90". Note: The program ‘‘Triangle’’ by Jonathan R. Shewchuk (University of California, Berkeley, United States, http://www.cs.cmu.edu/~quake/triangle.html ) was used in this study for the computation of Delaunay triangulations and Voronoi diagrams. Four programs are common to the contour line analysis for both the delination of drainage basins and the construction of flow nets: 1) shp2dat_node2.f90 2) Triangle by Jonathan R. Shewchuk 3) crust_skel_vrnode.f90 4) orgskel_shp.f90 The drainage basin delineation is then performed by: 5) bbt_plus.f90 Instead, the flow net constraction is performed by: 5) cb_dem_eqraf_new.f90 A short descritpion of each program and its input and output files are reported below: shp2dat_node2.f90 ----------------- This program read the contour lines saved in shapefile format and store them in ascii files. Input files: contour.shp input.dat quote.dat contour.shp: shapefile polyline. Each record of this file must have a singlepart feature. To convert multipart to singlepart feature class is possible to use the ESRI tool "Multipart To Singlepart" placed in ArcToolbox/Data Managment Tools/Features. For more details about the single and multipart feature classes see the ArcGIS Desktop Help. input.dat: N_poly,M,passo (Number of polyline, Max dimension of arrays, Contour interval) quote.dat: obtained opening in MS Excel the contour.dbf file and saving only the field "quote" (elevation) in txt file format. Output files : cariso.dat iso.dat pdvet.dat qdvet.dat ptraf.pts nodes.node iso.dat: each record is a contour line described with the same set of points like in contour.shp enriched with other sample points as described in section 2.5, Moretti and Orlandini (2008). cariso.dat: each record is composed by features of the corresponding contour line: elevation; number of points; =0 contour is closed, =1 contour is open. The number of record is the id of contour line used in ip.dat and par_flow.dat nodes.node: The same set of points stored in iso.dat but saved with the format needed by the program "Triangle". pdvet.dat, qdvet.dat and ptraf.pts are auxiliary files not required for the analysis described here. Triangle.exe ------------ This program is free to download at http://www.cs.cmu.edu/~quake/triangle.html It generates Delaunay triangulations and Voronoi diagrams. For this analysis is suggested the following running options: ./triangle.exe nodes.node -ev Input files: nodes.node Output files: nodes.1.edge nodes.1.ele nodes.1.v.edge nodes.1.v.node More details about structure of input and output files are reported in the on-line help. crust_skel_vrnode.f90 --------------------- The Delaunay triangulations and Voronoi diagrams are elaborated and the skeleton is defined. Input files: nodes.1.edge nodes.1.ele nodes.1.v.edge nodes.1.v.node Output files: skeletonx.dat skeletony.dat branchx.dat branchy.dat qskel.dat qbran.dat dim_cs.dat briso.dat iddelbran.dat iso_qb_min.dat skeletonx.dat: x coordinates of each couple of vertices defining a single segment belonging to the skeleton stems. skeletony.dat: y coordinates of each couple of vertices defining a single segment belonging to the skeleton stems. qskel.dat: elevation of each segment belonging to the skeleton stems. branchx.dat: x coordinates of each couple of vertices defining a single segment belonging to the skeleton branch. branchy.dat: y coordinates of each couple of vertices defining a single segment belonging to the skeleton branch. qbran.dat: elevation of each segment belonging to the skeleton branch. dim_cs.dat: old variable always zero, # of segments belonging skeleton stems, # of segments belonging skeleton branches, # of junction between branches and contour lines. briso.dat: for each branch, x and y coordinates of the two vertices defining the last branch segment joined to the contour line and the elevation. iddelbran.dat: for each branch in this file are stored the identification number of the two Delauny triangles that it belongs. iso_qb_min.dat: min elevation of baranches that are joined to a contour line. orgskel_shp.f90 --------------- The skeleton stems and skelton branches are structured. Input files: skeletonx.dat skeletony.dat branchx.dat branchy.dat qskel.dat qbran.dat dim_cs.dat cariso.dat iso.dat input.dat briso.dat iddelbran.dat iso_qb_min.dat Output files: carskel.dat skel.dat carbran.dat bran.dat feat.dat carcs.dat cs.dat bran.shp skel.shp feat.dat: old variable always zero, # of records of carskel.dat and skel.dat, # of records of carbran.dat and bran.dat skel.dat: skeleton stems. carskel.dat: features of skeleton stems. For each element are stored: quota: elevation. n_nodi: number of nodes belonging the polyline defining each skeleton stem. n_bran: number of branches releted with each skeleton stem. bran.dat: skeleton branches. carbran.dat: features of skeleton branches. For each element are stored: quota: elevation. n_nodi: number of nodes belonging the polyline defining each skeleton branch. ord_ram: order of branching id_mainbr: identification number of main branch, if =0 is a main branch directly connected with the skeleton stem. el_mainbr: number of element of the main branch, where the current branch is connected. cima: if =1 the branch is a peak, if =0 isn't. sella: if >0 the branch is a saddle, if =0 isn't. up_down: if =1 the branch connect a skeleton stem with the upward conotur line, if =2 the branch connect a skeleton stem with the downward conotur line. cs.dat: In this file are stored the characteristic branch (section 2.2 Moretti and Orlandini 2008). carcs.dat: features of characteristic branches. Note: ---- At the end of the sequence of program reported above the Contour Lines Analysis is completed and contour lines are enriched with the skeleton. With this set of data is possible to carry out both the drainage basin delineation (with bbt_plus.f90) and the flow net construction (with cb_dem_eqraf_new.f90 ) as described below: bbt_plus.f90 ------------ Drainage basin delineation. Input files: ip.dat cariso.dat iso.dat carskel.dat skel.dat carbran.dat bran.dat feat.dat ip.dat: In this file are stored coordinates and values of the variables that define the two starting points. In particular: 1. left point x coordinate 2. left point y coordinate 3. right point x coordinate 4. right point y coordinate 5. downslope left contour line id 6. downslope right contour line id 7. = 0 if the left pt is between the skeleton stem and the upslope contour line = -1 if the left pt is between the skeleton stem and the downslope contour line 8. = 0 if the right pt is between the skeleton stem and the upslope contour line = -1 if the right pt is between the skeleton stem and the downslope contour line 9. contour interval Output files: co_baxy.dat coba_plus.shp co_baxy.dat: x and y coordinates of each vertex belonging to the drainage divide from the two assigned points defined in ip.dat coba_plus.shp: shape file polygon of the drainage basin defined by the drainage divide. cb_dem_eqraf_new.f90 -------------------- Flow net construction. Input files: par_flow.dat cariso.dat iso.dat carskel.dat skel.dat carbran.dat bran.dat feat.dat carcs.dat cs.dat iso_qb_min.dat par_flow.dat: lambda,coeff_m,coeff_p,passo_iso,quota_min,quota_max,id_iso_start lambda= spacing between flow lines. coeff_m= factor that multiplies lambda. If the distance between start points are less then coeff_m*lambda the flow line is terminated. coeff_p= factor that multiplies lambda. If the distance between start points are greater then coeff_p*lambda a new flow line is inserted. passo_iso= contour interval. quota_min= lower contour line elevation. quota_max= higher contour line elevation. id_iso_strat= identification number of the lower contour line; the id is the number of the record in file cariso.dat Output files: flow_line.shp flow_line.shp: shape file polyline where the flow net is stored. WARNING ------- All the programs described above are in progress: major improvements for merging and clarifying procedures, renaming variables and translating comments in English are planned but not yet done. For example: If the starting contour lines shape file is not well clipped the procedure described in the paragraph 28 of section 2.5 in Moretti and Orlandini (2008) and implemented in shp2dat_node2.f90 may fails and the delineation of drainage basin may be incorrect. With the current version of the program for the flow net construction (cb_dem_eqraf_new.f90) it is not possible to reproduce Figure 6a reported in Moretti and Orlandini (2008) because we have developed two different routines: a first one for the flow net construction around a simple saddle like in Figure 6a and a second one for more complex morphology analysis like that reported in Figure 6b of Moretti and Orlandini (2008). The first one is not implemented in the new program cb_dem_eqraf_new.f90. 5. Examples ----------- Two examples are provided. cla\examples\basin_delineation in this folder are stored input and output files that beginning from the contour line, are produced with the procedure described above and summerize below: 1) shp2dat_node2.exe 2) Triangle.exe 3) crust_skel_vrnode.exe 4) orgskel_shp.exe 5) bbt_plus.exe The file contour_multipart.shp is done only for showing the difference between multipart and singlepart. The subfolder ip contains two ip*.dat that shows different cross sections. a) in ip1.dat the two starting points are placed between the skeleton stem and the upslope contour line. b) in ip2.dat the two starting points are placed between the skeleton stem and the downslope contour line. The bbt_plus have been runned with both the ip*.dat input files and the output file coba_plus.shp are renamed coba_plus_ip*.shp. cla\examples\flow_net in this folder are stored input and output files that beginning from the contour line, are produced with the procedure described above and summerize below: 1) shp2dat_node2.exe 2) Triangle.exe 3) crust_skel_vrnode.exe 4) orgskel_shp.exe 5) cb_dem_eqraf_new.exe 6. References ------------- Moretti, G., and S. Orlandini (2008), Automatic delineation of drainage basins from contour elevation data using skeleton construction techniques, Water Resour. Res., 44, W05403, doi:10.1029/2007WR006309.