Personal tools
Home » Members » scw » Basins Generation Process
Document Actions

Basins Generation Process

by Shaun Walbridge last modified 06-07-2011 14:46

Methods and process of developing a validated hydrologic basins dataset.

Basins Generation Process


Basin Overview

Prefix Cell Count (M) Cell Size (m) area (M km2) fills (at 100m) # of Basins
Africa af 4.65 922.8569813
Asia as 8.85 779.4213321
Australia au 0.98 908.9955774
Europe eu 9.64 702.5651599
North America na 4.44
779.4213321 27.0
Pacific pa 0.04 908.9955774 0.4
South America sa 2.19
Asia Islandsasi     28790

Clumping VMAP data

General VMAP information

Layer contents
Layer source
hyc=6 intermittent streams
f_code=BH000, hyc=8 lakes
f_code=BH000, hyc=6 intermittent lakes  'inwatera@hydro(*)_area'
f_code=BH090 floodland

Clumping method

Step One: ArcMap

  • SELECT vmap_inwatera WHERE f_code = BH000 & hyc = 8
    Export result to lakes.shp

    SELECT vmap_watrcrsl WHERE hyc = 8
    Export to streams.shp
  • Set extent with Spatial Analyst to be the same as the basins layer.
  • Buffer <1/2 cellsize> lakes.shp lakes_buffered.shp
    FeatureToRaster lakes_buffered.shp vmap_lakes
    FeatureToRaster streams.shp vmap_streams
    SingleOutputMapAlgebra "con(isnull(vmap_lakes), con(vmap_streams > 0, 1), 1)" \
    vmap_merge vmap_streams;vmap_lakes

Step Two: GRASS in=/cygdrive/f/dev/asia/as_setest/vmap_merge out=vmap_se_merge location=vmap_se
r.clump in=se_vmap_merge out=se_vmap_clump
r.out.arc in=se_vmap_clump out=/cygdrive/f/dev/wgs_globe/asia/as_setest/se_vmap_clump.asc

Step Three: ArcMap

ASCIIToRaster se_vmap_clump.asc vmap_clump INTEGER

Process basins

Delineating true boundries by visual inspection

Input datasets

Layer name
DEM, elevation inflated 1000m
se_norm DEM, true elevation values
se_fill DEM, processed with 100m fill, inflated 1000m
DEM, processed with 100m fill, true elevation values
se_facc flow accumulation
flow direction
fill depth: fill elevation - true elevation

  • Reproject GTN-basins into your locally centered Lambert projection
  • Clip GTN basins to continental mask
  • Generate a link to the local Landsat scene:
    1. Connect to Cabrillo
    2. [user@cabrillo] ln -s /mnt/raid/landsat/ ~/landsat
    3. Mount as network drive
    4. Use ECW driver to use data directly from mapped drive.
  • Create new ArcMap document with the following layers:
    1. vmap_clump: color by unique value
    2. gtn_basins: use provided style
    3. basins_elim_x3: use provided style
    4. se_fnorm: use provided style
    5. ... 
  • Add new column to basins layer `inspect', type: short Int.

Inspection key:

Key value
1 valid major basin
2 VMAP shows two or more basins
3 valid internal basin
4 GTN shows two or more basins
5 valid internal lake drainage

Manually forcing flow directions

Manual inspection of the basins and flow accumulation, in conjuction with Landsat imagery and the VMAP data, will show specific sections where flow was misrouted and thus the problem areas that cause erroneous basins. In order to adjust basins, it is necessary to adjust the DEM itself.

The modifications are stored in 4 seperate shapefiles:

  • burns (line; lower the DEM to aricifically create a river channel )
  • breaks (line; rasie the DEM to artificially create a "dam" or ridgeline )
  • drains (point; nullify the DEM at the given point in order to force water to an internal drain or to remove a false land bridge)
  • fills (polygon; create DEM values where there were originally nulls to avoid pour points occuring too far inland )

Re-running Flow

first create a new directory in the basins tree (something like "build_try1"). This is where the new re-run model will be output.

"Bumping" the basins (burns, breaks, drains, fills)

Using the BasinBump script, merge your modifications with the DEM. You'll need an input DEM, and your folder of modifications.  From the command line, run script changing the output prefix name according to the continent and output directory you want.

 python <prefix> <elevation> <modifications dir> <output dir>

Generating the basins

In ArcInfo GRID, run the basin_generate AML substituting your output prefix below:

basingenerate {outputname}_dem 800 12 {outputname}

Postprocessing the Basins

After the final basins layer has been generated, there are still a few more steps to finalize the product:

Eliminate Minor Basins

  • Define Projection: se_basin_nib.shp
  • Open se_basin_nib.shp Attribute table, Create new column `area', type double
  • Calculate Values and check the Advanced Box:
    Dim dblArea as double
    Dim pArea as IArea
    Set pArea = [shape]
    dblArea = pArea.area / 1000000
    In the result box at the bottom, type dblArea
  • Select By attributes: WHERE "area" < 12
  • Eliminate_management af_basin_n C:\Workspace\srtm_basins\af_basins\build_try0\af_basin_elim12.shp LENGTH
  • Rinse, lather, repeat: do last two steps twice.  Final file should be se_basin_elim_x3.shp.
  • Alternatively, you can keep repeating the eliminate until there is no change in the number of basins (with area < 12) between iterations.

Automating the tagging of internal drains

Using the raster seamask as the definition of the shoreline:

  1. convert raster to features (polygon)
  2. convert the polygons into lines
  3. assign the srs
  4. buffer them by 2km
  5. Then select all basins intersect the shoreline buffer.
  6. Invert the selection
  7. Add an attribute column to the basins layer ('inspect'; short int)
  8. assign inspect = '3' to all selected basins.

Creating Pour Points

Run, i.e:

python sa_bas_n_elim_2.shp sa_facc f:\dev\wgs_globe\samerica\bump_v6 sa
This script will generate pour points for the maximum flow accumulation within each basin.  Unfortunately, some islands and other small coastal basins have multiple pours, since the facc is the same in multiple locations. Often, this is a side-effect of nibbled or eliminated locations. To modify the basins for these psuedo pours, the followed hackery is needed:
  1. Run HawthsTools: PointInPolygon on: basins (sa_bas_n_elim2) + drains (sa_processed) to field PNTPOLYCNT
  2. Run HawthsTools: Intersect Point Tool: drains + basins: choose ID, GRIDCODE, PNTPOLYCNT
  3. Export points with more than one drain to a new dataset:
    • Select by Attributes: "PNTPOLYCNT" > 1

    • Data: Export data: create new dataset sa_multi_points.shp
  4. Export points with a single drain to a seperate dataset:
    • Selection: Switch selection
    • Data: Export data: create new dataset sa_single_points.shp
  5. For the multiple drains layer, we now want to select an arbitrary pour, so we can maintain a 1:1 basin to pour ratio:
    • Open attribute table on sa_multi_points.shp
    • Right click ID column, select Summarize
    • Summarize on POINTID: minimum
    • Save as sa_multi_sum.dbf
  6. Join: sa_multi_points `POINTID' to sa_multi_sum `Min_POINTID'
  7. Select by Attributes: "sa_multi_sum.Min_POINTI" IS NOT NULL
  8. Double check that basin count matched pour point count. South America was short four pour points due to nibble mishaps.
  9. Export as sa_multi_fixed_points.shp
  10. In ArcCatalog copy sa_single_points.shp  and rename sa_merged_points.shp  Run Append (Management) from Data Management Toolbox.  Input Dataset: sa_multi_fixed_points.shp, Target Dataset: sa_merged_points.shp Schema Type: NO_TEST.  Check sa_merged_points.shp Attribute Table to make sure number of Records is equal to number of basins in sa_bas_n_elim2
  11. A few basins were falsely collapsed -- still short six after editings, rerun step #1 on merged points + basins.
  12. Finally, make sure all the basins hit the coastline: Select By Location: sa_merged_points.shp "intersect" sa_mask_ocean <buffer 500m>.

Basins Naming

After all post-processing is completed, we need to name the basins with unique identifiers.

  1. Open the completed basins layer, right click: Open Attribute Table
  2. Options: Add Field: Name the field basin_id, type text, length 8

  3. Calculate Values on basin_id: Choose advanced, then paste this snippet, replacing "sa" with your continent abbreviation:
    Dim id as string
    Dim zeros as string
    Dim cont as string
    Dim basinId as string

    cont = "sa"
    id = [ID]
    zeros = "00000"

    basinId = cont & "_" & Left(zeros, 5 - Len(id)) & id

Built with Plone