/* NAME 
/*
/* STREAMNETWORK.AML
/*
/*----------------------------------------------------------------------------
/*
/* AUTHOR 
/* Jacek Blaszczynski, Physical Scientist, BLM NARSC, 04/99, Original Coding
/*
/*----------------------------------------------------------------------------
/*
/* DESCRIPTION:
/*
/* This AML develops a grid and a coverage stream network from digital
/* elevation data. It first removes all the internal drainages or sinks from the
/* terrain model. Then it routes the flow of water finding a flow direction 
/* for every cell of the elevation grid following the steepest paths between cells.  
/* Flows from cell to cell on the same flowpath result in a flow accumulation
/* map in which the value for each cell represents the accumulated flow along 
/* a particular path.  We can use flow accumulation map to delineate stream 
/* networks by setting a lower limit to our flow accumulation value since
/* that flow accumulation on hillslopes is lower than in the stream channels.  
/* Depending on the flow accumulation threshold we decide to use our network 
/* "grows" tributaries in the upland direction (toward the hill crests or 
/* interfluves).  The lower the flow accumulation threshold used, the 
/* more detailed is the resulting network.  In the final grid product each stream 
/* segment or reach is provided with a different value.  Based on this grid 
/* map a line coverage is also produced representing the network.  The program
/* then proceeds to do Strahler and Shreve stream network calculations
/* and attaches the values of stream order to the grid and line vector
/* stream data.
/* 
/* WARNING:  If the extent of the project area that is being processed
/* is large or flow thresh low so that more than 100,000 of links get generated
/* the program might not be able to proceed.  A streamlined version of
/* STREAMNETWORK.AML which avoids this problem is available from the author.
/* The program REACHSHEDMODEL contains this streamlined version and can
/* be used to generate stream network data as well.
/*
/* NOTE:
/* The "z" suffix on product maps indicates that this operation was 
/* performed on a filled DEM.
/*

/*----------------------------------------------------------------------------
/* 
/* VARIABLES:
/*
/* source_dem -- the source digital elevation model (DEM)used for delineation of watersheds
/* flow_thresh -- or the "flow accumulation threshold" is the number of cells that represent 
/* the flow accumulation level at which streams and channels are delineated to form basis 
/* for watershed delineation
/* 
/*
/*----------------------------------------------------------------------------
/*
/* DEFAULT GRID NAMES:
/*
/* filled (filled DEM), flowdirz (flow direction grid), flowcountz (flow accumulation grid),
/* mstrms%flow_thresh% (stream network grid), mlinks%flow_thresh% (stream link real no grid),
/* mlinks%flow_thresh%I (stream link integer grid), mstrahl%flow_thresh%I (Strahler stream 
/* order integer grid), mshreve%flow_thresh%I (Shreve stream order integer grid),
/* intlnkstrah (intersection of link and Strahler stream order map), intlnkshrev 
/* (intersection of link and Shreve stream order map), mstrms'flow_thresh'l (new stream 
/* network line map) 
/*
/*----------------------------------------------------------------------------
&severity &error &routine bailout
&args source_dem flow_thresh
/* Check arguments
&call chkargs

/* Do the work
&call main

&call exit
&return /* End of streamnetwork.aml

/*--------------
&routine chkargs
/*--------------
&if [SHOW PROGRAM] <> GRID &then &do
  &type This aml must be run from GRID.
  &call bailout
&end
&do arg &list source_dem flow_thresh
  &if [NULL [VALUE %arg%]] &then &do
    &call usage
    &call bailout
  &end
&end
&if not [EXISTS %source_dem% -GRID] &then &do
  &type Input grid: %source_dem%, does not exist.
  &call bailout
&end
&else &do
  &if [LENGTH %source_dem%] > 13 &then &do
    &type Grid names must be 13 characters or less.
    &call bailout
  &end
&end

&s thresh_limit = 100000
&s datatype = [TYPE %flow_thresh%]
&if not %datatype% in {-1,-2} &then &do
  &call thresh_error
  &call bailout
&end
&else &do
  &if %flow_thresh% > %thresh_limit% &then &do
    &call thresh_error
    &call bailout
  &end
&end
&return /* End of routine chkarg


/*-----------------
&routine thresh_error
/*-----------------
&type Flow threshold values must be an integer value less than or equal to 100,000 cells
&return /* End of routine thresh_error

/*--------------------------------------------
&routine main
/*--------------------------------------------

/*--------------------------------------------
/* Set the environment
/*--------------------------------------------

&echo &brief
&terminal 9999
display 9999 2 position cr
setwindow %source_dem% %source_dem%

/*--------------------------------------------
/* Display terrain model
/*--------------------------------------------

mape %source_dem%
gridpaint %source_dem% value linear wrap gray 

/*--------------------------------------------
/* Clean data with the same flow_thresh from  
/* previous runs of this AML
/*--------------------------------------------

&if [EXISTS mstrms%flow_thresh% -GRID] &then &do
  &type Existing stream and linkage grids and coverage made previously with the same flow      
  &type accumulation thresholds will be deleted
  kill mstrms%flow_thresh% all
&end
&if [EXISTS mlinks%flow_thresh% -GRID] &then &do   
  kill mlinks%flow_thresh% all
&end
&if [EXISTS mstrms%flow_thresh%l -COVER] &then &do 
  kill mstrms%flow_thresh%l all
&end
&if [EXISTS mlinks%flow_thresh%I -GRID] &then &do
  kill mlinks%flow_thresh%I all
&end
&if [EXISTS mstrahl%flow_thresh%I -GRID] &then &do
  kill mstrahl%flow_thresh%I all
&end
&if [EXISTS mshreve%flow_thresh%i -GRID] &then &do
  kill mshreve%flow_thresh%i all
&end
&if [EXISTS intlnkstrah -GRID] &then &do
  kill intlnkstrah all
&end
&if [EXISTS intlnkshrev -GRID] &then &do
  kill intlnkshrev all
&end
&if not [EXISTS filled -GRID] &then &do 
  fill %source_dem% filled
&end
&if not [EXISTS flowdirz -GRID] &then &do 
  flowdirz = flowdirection (filled)
&end
&if not [EXISTS flowcountz -GRID] &then &do
  flowcountz = flowaccumulation (flowdirz)
&end

linesize .025
/*-----------------------------------------------
/* Display flow direction map
/*-----------------------------------------------

gridshades flowdirz 

/*-----------------------------------------------
/* Develop stream network grid
/*-----------------------------------------------

mstrms%flow_thresh% = con(flowcountz < %flow_thresh%, ~
                          setnull (flowcountz), 5)
buildsta mstrms%flow_thresh%
gridshades mstrms%flow_thresh%

/*------------------------------------------------
/* Develop a stream network segment (link) map
/*------------------------------------------------

mlinks%flow_thresh% = streamlink (mstrms%flow_thresh%, flowdirz)
buildsta mlinks%flow_thresh%
gridshades mlinks%flow_thresh%
mlinks%flow_thresh%i = int (mlinks%flow_thresh%)
buildvat mlinks%flow_thresh%i
gridshades mlinks%flow_thresh%I

/*----------------------------------------------------
/* Create a separate attribute (item) to contain link 
/* numbers for each segment of the stream network
/*----------------------------------------------------

arc tables
  sel mlinks%flow_thresh%i.vat
  additem mlinks%flow_thresh%i.vat link-id 8 8 f 2
  calculate link-id = value
quit

/*----------------------------------------------------
/* Create a line vector stream coverage from
/* the stream network link grid
/*----------------------------------------------------

mstrms%flow_thresh%l = streamline (mlinks%flow_thresh%i, flowdirz, link-id)

/*----------------------------------------------------
/* Create a "value" attribute or item for the line 
/* modeled stream coverage to allow for joining of
/* attributes between the grid and coverage stream maps
/*-----------------------------------------------------


arc tables
  sel mstrms%flow_thresh%l.aat 
  additem mstrms%flow_thresh%l.aat value 4 4 b
  calculate value = link-id
quit

/*------------------------------------------------------
/* Calculate Strahler stream order for each stream link
/* and attach the values to the stream network grid
/*------------------------------------------------------

mstrahl%flow_thresh%i = int (streamorder (mlinks%flow_thresh%i, ~
   flowdirz, strahler))
buildvat mstrahl%flow_thresh%i
intlnkstrah = mlinks%flow_thresh%I cand mstrahl%flow_thresh%I
buildvat intlnkstrah
arc joinitem mlinks%flow_thresh%I.vat intlnkstrah.vat ~
   mlinks%flow_thresh%i.vat value
gridshades mstrahl%flow_thresh%i

/*-------------------------------------------------------
/* Calculate Shreve stream order for each stream link
/* and attach the values to the stream network grid
/*-------------------------------------------------------

mshreve%flow_thresh%i = int ( streamorder ( mlinks%flow_thresh%, ~
   flowdirz, shreve ) ) 
buildvat mshreve%flow_thresh%I 
gridshades mshreve%flow_thresh%I
intlnkshrev = mlinks%flow_thresh%I cand mshreve%flow_thresh%i
arc joinitem mlinks%flow_thresh%I.vat intlnkshrev.vat ~
   mlinks%flow_thresh%i.vat value

/*-------------------------------------------------------
/* Join stream order values of the stream network grid
/* to the vector stream network coverage
/*--------------------------------------------------------

arc joinitem mstrms%flow_thresh%l.aat mlinks%flow_thresh%I.vat ~
  mstrms%flow_thresh%l.aat value
arc build mstrms%flow_thresh%l line

linecolor 5
arcs mstrms%flow_thresh%l
&return /* end of routine main

/*------------
&routine usage
/*------------
&type Usage: &r streamnetwork <source_dem_grid> <flow_threshold>
&return /* End of routine usage

/*--------------
&routine exit
/*---------------
&return /* End of routine exit

/*--------------
&routine bailout
/*--------------
&severity &error &ignore
&call exit
&return &error Bailing out of streamnetwork.aml...