/*----------------------------------------------------------------------------
/*
/*NAME: polygonsheds.aml
/*
/*P----------------------------- Purpose -------------------------------------
/*
/* Converts a grid of modeled watersheds into a polygonal coverage
/*
/* The problem with vectorizing a watershed grid modeled from terrain data and 
/* making a polygonal coverage out of it are slivers, which is something one 
/* wouldn't expect when source data is in grid or raster format.  When GRID 
/* creates watersheds it does so basing their boundaries on the location of 
/* the watershed outlet, cells flowing into, and the location of the next 
/* watershed outlet.  The program "grows" watersheds following flow direction 
/* until all the cells that flow into any particular outlet are included, that 
/* is unless they flow to another outlet.  This results in a string of usually 
/* two or three cells attached at the edge of the greater watershed area.  
/* When this is vectorized these individual cells, often attached to each 
/* other diagonally by their corners, become separate polygons which create 
/* many tiny new watersheds.  This AML provides a simple way to remove the 
/* diagonally adjacent cells by setting an area threshold for watersheds and 
/* removing ones that fall under the threshold while the map is still in grid 
/* format.
/*
/*A---------------------------- Arguments ------------------------------------
/*
/* input_shed   raster map of watersheds
/* output_shed  polygon map of watersheds
/* area_thresh  usually an area equal to 10000 sq. meters is adequate: 
/*              represents small areas of watersheds delineated by diagonally 
/*              adjacent cells which are to be removed
/*H----------------------------- History -------------------------------------
/*
/* Jacek Blaszczynski 04/21/97 Original Coding
/*                    05/14/99 Corrections and edits (adding display functions)
/* 
/*E===========================================================================
&severity &error &routine bailout 
&args input_shed output_shed area

/* Check arguments
&call chkargs

/* Do the work
&call main

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

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


&s arealimit = 10000000
&s datatype = [TYPE %area%]
&if not %datatype% in {-1,-2} &then &do
  &call area_error
  &call bailout
&end
&else &do
  &if %area% > %arealimit% &then &do
    &call area_error
    &call bailout
  &end
&end
&return /* End of routine chkargs

/*-----------------
&routine area_error
/*-----------------
&type Area value must be an integer or real value less than or equal to 20,000 sq. meters
&return /* End of routine area_error

/*-----------
&routine main
/*-----------
&echo &brief
DISPLAY 9999 2 position cr
clear
setwindow %input_shed% %input_shed%
mape %input_shed%
gridshades %input_shed%
&s editshed = [SCRATCHNAME -DIRECTORY]
%editshed% = con (zonalarea (regiongroup (%input_shed%)) > %area%, %input_shed%, -999)

&s mskgrd = [SCRATCHNAME -DIRECTORY]
%mskgrd% = con (%editshed% == -999, setnull(%editshed%), 0)

&s tempgrid = [SCRATCHNAME -DIRECTORY]
%tempgrid% = nibble(%editshed%, %mskgrd%, DATAONLY)

%output_shed% = gridpoly (%tempgrid%)
rename %tempgrid% %output_shed%c
gridshades %output_shed%c
arc build %output_shed% poly
linesize .025
linecolor 9
arcs %output_shed%
kill %editshed% all
kill %mskgrd% all
&echo &off
linecolor 1

&return /* End of routine main
 

/*------------
&routine usage
/*------------
&type Usage: &r polygonsheds <watershed_grid_to_be_filtered><output_watershed_grid><area_threshold>
&type The area threshold needs to be under 20000 sq. meters.
&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
&type When finished, please LIST watersheds.  Small areas with 
&type grid-code of -9999 were NODATA areas in the original grid.
&type These false watersheds will have to be selected out. To make
&type a grid watershed map that corresponds to polygonal watersheds
&type reselect the polygon map to exclude areas with grid-code = -9999
&type and rasterize the resultant map so that polygon-id item value
&type is the same as value for each area in the grid map.
&return /* End of routine exit

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