Basins Generation Process

Resources:

• Wikipedia river entries
• VMAP naming conventions
• Source code
  
  

Basin Overview

Name	Prefix	Cell Count (M)	Cell Size (m)	area (M km2)	fills (at 100m)	# of Basins
Africa	af	4.65	922.8569813	39.6	295826	17375
Asia	as	8.85	779.4213321	52.2	~600000	12327
Australia	au	0.98	908.9955774	8.1	46978	8226
Europe	eu	9.64	702.5651599	47.6	419772	23872
North America	na	4.44	779.4213321	27.0	299513	58082
Pacific	pa	0.04	908.9955774	0.4	3280	5350
South America	sa	2.19	908.9955774	18.1	179023	10883
Asia Islands	asi					28790

Clumping VMAP data

General VMAP information

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

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

r.in.gdal 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	Description
se_ult_clip	DEM, elevation inflated 1000m
se_norm	DEM, true elevation values
se_fill	DEM, processed with 100m fill, inflated 1000m
se_fnorm	DEM, processed with 100m fill, true elevation values
se_facc	flow accumulation
se_fdir	flow direction
as_fdepth	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	Description
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 basinBump.py script changing the output prefix name according to the continent and output directory you want.

 python basinBump.py <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 basinPours.py, i.e:

python basinPours.py 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