/* REACHSHEDMODEL.AML
/* Author: Jacek Blaszczynski 03/11/99 Original Coding
/*
/*
/*------------------------------------------------------------------------------
/*
/* This AML delineates watersheds for watershed outlets located at stream 
/* junctions and sources.  The streams are modeled from Digital Elevation
/* Model from a flow accumulation map starting at cells which have the user
/* specified flow accumulation threshold.  The program first creates basic
/* hydrologic terrain analysis grid layers of a filled DEM ("filled"),
/* flowdirection map for the filled DEM ("flowdirz" - the "z" suffix
/* indicates the operation was performed on a filled DEM), and a 
/* flowaccumulation map ("flowcountz").  After creating a stream map
/* from a flow accumulation grid the program creates a line stream map
/* from the stream grid.  The endpoints of stream lines and points of
/* intersection (nodes) are identified and made into a point coverage.
/* The point coverage is then rasterized to serve as outlets for watershed
/* delineation.  The GRID "watershed" function calculates the watersheds
/* based on the outlet and the flow direction map.  
/* 
/*-----------------------------------------------------------------------------
&severity &error &routine bailout
&args source_dem flow_thresh

/* Check arguments
&call chkargs

/* Do the work
&call main

&call exit
&return /* End of reachshedmodel.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
  &if [LENGTH %source_dem%] > 13 &then &do
    &type Grid names must be 13 characters or less.
    &call bailout
  &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 display environment and
/* clean data  with the same flow_thresh
/* generated in previous runs of this AML
/*----------------------------------------

&echo &brief
&TERMINAL 9999
DISPLAY 9999 2 position cr
setwindow %source_dem% %source_dem%
mape %source_dem%
gridpaint %source_dem% value linear wrap gray

&if [EXISTS mstrms%flow_thresh% -GRID] &then &do 
  kill mstrms%flow_thresh% all
&end
&if [EXISTS mlinks%flow_thresh%i -GRID] &then &do   
  kill mlinks%flow_thresh%i all
&end
&if [EXISTS mstrms%flow_thresh%l -COVER] &then &do 
  kill mstrms%flow_thresh%l all
&end
&if [EXISTS mnodes%flow_thresh% -COVER] &then &do  
  kill mnodes%flow_thresh% all
&end
&if [EXISTS moutl%flow_thresh% -GRID] &then &do 
  kill moutl%flow_thresh% all 
&end
&if [EXISTS moutl%flow_thresh%reg -GRID] &then &do 
  kill moutl%flow_thresh%reg all
&end
&if [EXISTS msheds%flow_thresh% -GRID] &then &do 
  kill msheds%flow_thresh% all
&end
&if [EXISTS mstrahl%flow_thresh%I -GRID] &then &do
  kill mstrahl%flow_thresh%I all
&end
&if [EXISTS intlnkstrah -GRID] &then &do
  kill intlnkstrah all
&end
&if [EXISTS mshreve%flow_thresh%i -GRID] &then &do
  kill mshreve%flow_thresh%i 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
gridpaint flowdirz value linear wrap gray
&if not [EXISTS flowcountz -GRID] &then &do
  flowcountz = flowaccumulation (flowdirz)
&end

gridnodatasymbol transparent
markercolor 7
linesize .025

/*-----------------------------------------------------
/* Create a stream network grid and a vector stream map
/*-----------------------------------------------------

mstrms%flow_thresh% = con(flowcountz < %flow_thresh%, setnull (flowcountz), 5)
mstrms%flow_thresh%l = streamline (mstrms%flow_thresh%, flowdirz)
arc build mstrms%flow_thresh%l line
linecolor 5
arcs mstrms%flow_thresh%l 

/*-----------------------------------------------
/* Create a point coverage of nodes or endpoints
/* and intersections of streams to use as outlets
/* for watersheds and generate watersheds
/*-----------------------------------------------

arc nodepoint mstrms%flow_thresh%l mnodes%flow_thresh%
points mnodes%flow_thresh%
setcell %source_dem%
setwindow %source_dem% %source_dem%
moutl%flow_thresh% = pointgrid (mnodes%flow_thresh%)
buildsta moutl%flow_thresh%
gridshades moutl%flow_thresh%
moutl%flow_thresh%reg = regiongroup (moutl%flow_thresh%)
buildsta moutl%flow_thresh%reg
gridshades moutl%flow_thresh%reg
msheds%flow_thresh% = watershed (flowdirz, moutl%flow_thresh%reg)
buildsta msheds%flow_thresh%
gridshades msheds%flow_thresh%
arcs mstrms%flow_thresh%l
&echo &off

&return /* end of routine main
/*
/*------------
&routine usage
/*------------
&type Usage: &r reachshedmodel <source_dem_grid> <flow_threshold>
&type USE POLYGONSHEDS.AML TO CREATE A WATERSHED POLYGON COVERAGE !!!
&return /* End of routine usage

/*-----------
&routine exit
/*-----------
&do grd &list editshed mskgrd
  &if [VARIABLE %grd%] &then &do
    &if [EXISTS [VALUE %grd%] -GRID] &then &do
      KILL [VALUE %grd%] ALL
    &end
  &end
&end
&return /* End of routine "exit"

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