VIEWS: 8 PAGES: 6 POSTED ON: 7/8/2010
Prof. Dr. D. Schröder Page 1 GIS in Hydrology and Water ResourceManagement - ENWAT Tutorial: ArcGIS Tools for Hydrology – Hydrological Analysis of a Watershed In this tutorial, the tools for hydrologic analysis of ArcGIS for surface water will be introduced. These tools are part of the extension Spatial Analyst, which is an ArcGIS extension for raster analysis. Thus most of the analysis will be based on raster maps. The study area is a small sub basin of the Rems watershed, the “Strümpfelbach” sub basin. The main tasks of this tutorial are: Delineation of the Strümpfelbach sub basin Determination of the sub basin characteristics, and Calculation of a flow time map and a time area diagram for the outlet point As input data the following raster and vector layers are available: DEM with 25 m resolution georeferenced topographic raster map digitized streams from topographic map Corine Land Cover vector map Workflow: Remember: all the rasters we are working with in our analysis should have the same resolution and the same extent. Here all the analysis is based on the DEM. So it’s a good idea to use its resolution and extent for all the analysis. As we will create a number of new maps, it’s a good idea to use a strict file naming convention scheme! Step 1: preparation of the DEM For preparation our DEM for hydrological analysis, we should burn-in the existing stream network and fill the existing sinks. Burn in existing channels In flat areas, it is difficult for the accumulation algorithm to determine the flow channel path. To improve the result a “burn-in” of the existing stream network can be used, i.e. subtracting from the DEM a meaningful value for the pixel covered by the streams. So we need a raster, where all the stream pixels have a value of let’s say 2 m, where as all the other pixels have the valued 0 (Attention: not noData!) First we have to “rasterize” the stream ploylines, e.g. convert them to a raster: Conversion Tools To Raster Polyline to Raster Click the environment button to set the extent of the new raster map: Environment General settings Extent set same as layer DEM For the cell size, you should set the exact value of the DEM (25.170305836741), best copy it, e.g. from the cell size setting in the raster analysis settings! After applying the settings, a new raster layer will be created. The raster cells will have cell values according the field you have used as input (by default the FID field), all cells not covered by a stream will just have noData as value. Prof. Dr. D. Schröder Page 2 GIS in Hydrology and Water ResourceManagement - ENWAT Now we have to reclassify the raster. All the “stream” cells should have the value 2, all the others just 0, so that we can subtract this layer from the DEM layer. Spatial Analyst Tools Reclass Reclassify Now we can burn in our streams using map Algebra: Spatial Analyst Tools Map Algebra Single Output Map Algebra The syntax is straight forward, just subtract the two layers! Fill sinks In the next step we have to correct for errors in the DEM which will have effects on the hydrological analysis. Spatial Analyst Tools Hydrology Fill To check the effect of the fill of sinks, calculate the difference of the original DEM and with the filled sinks! Now we are ready for the hydrological analysis! Prof. Dr. D. Schröder Page 3 GIS in Hydrology and Water ResourceManagement - ENWAT Step 2: Hydrological analysis Flow direction map and flow accumulation map Now we will calculate the flow direction map and the flow accumulation. With that we can extract the (synthetic) stream network and delineate the watershed. Spatial Analyst Tools Hydrology Flow Direction (input: filled_DEM) Spatial Analyst Tools Hydrology Flow Accumulation (input: flow_direction) Extraction of the stream network There is no special tool for the extraction of the network. Instead we have to use map algebra that is creating a new raster layer containing only those pixels, which have an accumulation value bigger than a fixed threshold. For setting the threshold we should look at our topographic map, checking where here the streams start. On the other hand, we have to consider the resolution of our DEM of about 25 m, too. For the extracted streams later gully flow will be assumed, for all the other pixels sheet flow. But our Strümpfelbach is only a small brook, with a top with of about 5 m at its junction to its headwater. Thus we have to use some compromise for fixing a meaningful threshold. For Map Algebra, we have to use the conditional function con, with the general syntax: Con(<condition>, <true_expression>, <false_expression>) Here e.g. con(flowacc_flow2 > 1500,1,0) with the threshold set to 1500 accumulated pixels corresponding to a drainage area of about 94ha. As the true expression is 1, all pixel values with an accumulation bigger than 1500 will be assigned 1, all others 0. With the stream network we now are able to calculate: - different segments for the network (Spatial Analyst Tools Hydrology Stream Link) - stream order according to Strahler or Shreve (Spatial Analyst Tools Hydrology Stream Order) - convert the (raster) stream network to vector polylines (Spatial Analyst Tools Hydrology Stream to Feature) Watershed Delineation Now we are ready to delineate the watershed. In the Spatial Analyst, there are two tools: basin and watershed. The first one will create all sub basins for the whole study area, where as with the second you can define watershed upstream to defined outlet points (pour points). So we will use watershed, as we want to delineate the watershed upstream to the outlet point, i.e. the junction of the Strümpfelbach to its headwater. The pour point hast to be located exactly on the pixel of the highest accumulation to determine the correct watershed. For that, the tool Snap Pour Point has to be used. Here the steps in detail: - use ArcCatalog to create a new shape file (geometry type Point) - digitize the pout point Prof. Dr. D. Schröder Page 4 GIS in Hydrology and Water ResourceManagement - ENWAT - use the snap to Pour point tool to drag the point exactly on the pixel with highest accumulation at the outlet. - delineate the watershed using the watershed tool - convert the raster to a (vector) polygon (Conversions Tool Raster to Polygon) Step 3: Calculate watershed/stream network characteristics After extracting the streams and delineating the watershed, we can calculate some general parameters of the watershed: - Channel length: (see also tutorial 1 on vector data!) add a new field to the table of your vector streams and use Calculate Geometry (only in ArcGIS 9.2! look in the help of the field calculator for older versions how to calculate the length). Select the main channel segments and use Statistics to get the sum of the segments (about 5 km) - Channel slope: use the identify tool to get the height values for the outlet point and the source of your main stream (2.7%) - Drainage area: calculate Geometry for the polygon or number of pixels x area of a pixel (949ha polygon, pixel 14981 x 625= 936ha) - Drainage density: Select all segments of your streams inside the polygon using Selection by location. Divide the length by the drainage area (0.6 17 1/km) - Hydrological length: Extract by mask from the (filled) DEM the pixels inside the watershed, recalculate flow direction map, use the tool Flow length, the maximum value is the hydrological length (6.3km). - Watershed slope: From the pixels inside the watershed take the difference of the highest and the lowest value and divide it by the hydrological length (480-247)/6300 = 3.6% - Length of the centroid: use Data Management Tools Features Feature to Point to calculate the centroid. Use the measurement tool to measure the approximate distance to the perpendicular toe point on the stream. (There are tools to make it better, but remember, a GIS is not a CAD, so it is rather cumbersome to de it exact!) (2.9km) - Shape Factor Ll = (LLca)0.3 2.39 km2/3 - Circularity ratio Fc = P/(4πA)0.5 (P=14392m; Fc=1.32) - Circularity ration Rc = A/Ao (U=2π r r = 2274 A0 = pi r2 A0 = 162 Rc =5.85) - Ruggedness Number = basin relief times drainage density.(m/km) 143 Step 4: Flow velocities Based on the flow velocities, the travel time from each pixel to the outlet can be calculated. Calculation of the velocity map is the most crucial step in the workflow. The flow velocity for each pixel depends on the slope and the roughness of the surface. To calculate the velocities, for instance the Manning-Strickler equation can be used. Here we have carefully to distinguish between sheet and gully flow, i.e. different parameters for the roughness as well as for the hydraulic radius should be used. The roughness itself is determined by the land use. In the literature table for the Manning value n can be found (Kst is the inverse of the Manning value n), for instance Haestad Methods: Floodplain modelling using HEC-RAS, 2003. Sheet flow Surface description Kst in m1/3s -1 Smooth surface (concrete, asphalt, gravel, or 90 bare soil) Fallow 20 Cultivated soils 10 Grass 5 Wood 17 Prof. Dr. D. Schröder Page 5 GIS in Hydrology and Water ResourceManagement - ENWAT Channel flow Surface description Kst in m1/3s -1 Minor streams, some grass and weeds 30 The Corine Land Cover layer can be used for assigning meaningful averaged roughness values by adding a new field (Kst) to the attribute table and edit it manually. After adding the values, we are ready to rasterize the map using the Kst values. Now we have to do some map algebra to calculate the flow velocities. In complex expression, It is a good idea to split up the calculation into different steps. Use map Algebra to calculate first the Kst x r 2/3 map, where r=0.03, so r2/3=0.0965 For modelling the channel flow, we have to keep in mind that the map resolution is 25m, so that only a part of a pixel is really covered by the channel. For the rest we will still have sheet flow. So again, some averaged value should be used as roughness and hydraulic radius. For the first order streams, let’s assume a top width of 1 m and a depth of 0.3, so that for a rectangular channel shape we have a wetted parameter of 1.6m and the hydraulic radius (A/Pw) becomes 0.1875 m. For second order streams assume a top width of 2 m and a depth of 0.5 m, i.e. a wetted parameter of 3m and a hydraulic radius 0.333m. As not the whole pixel will be covered by the channel, we use a “corrected” Kst of 20. Thus we have as kst r2/3 for our channels: first order 7 and second order 10. Now we can Reclassify the Strahler order map assigning the values above to the stream pixels, all other pixels should have the value 0. Now we have to combine the two maps with sheet and gully flow. One possibility is to first “burn out” the streams in the sheet flow map by using a binary stream map (1 for no stream, 0 for stream) and multiply this with the sheet flow map. Secondly just sum the two maps to get the final kst r2/3 map for combined sheet and gully flow. To get the final velocities map, we have to multiply this map with the square root slope map. Creating the slope map, you should take care of the units. Using percentage we have to divide the numbers first by 100, than we can multiply to get the final velocity map. Step 5 Calculating flow time and the simple IUH As the flow time is the quotient of path and velocity, we can use Flow length to calculate it. We have only to introduce the inverse of the velocity as weights. As a weight raster has to be an integer raster, we have first to convert it to integer. Dividing by 3600 we will get the travel time in h as an isochron map. The time area diagram we get from the histogram of this map. For an overview, use the histogram (accumulated) from the symbology dialog. If you need more details, export (copy and paste) the VAT of the raster layer to Excel. As the ordinates of the IUH are the differential changes of the drainage area with respect to time, the ordinate of the IUH are the gradients of the time area diagram and can be easily calculated. Prof. Dr. D. Schröder Page 6 GIS in Hydrology and Water ResourceManagement - ENWAT Flow time map Time area diagram (histogram from ArcGIS) 1400 1200 1000 pixel count 800 IUH in minutes (Excel based) 600 400 200 0 0 10 20 30 40 50 60 time (min)