&args basin maskgrid classdir zfactor
/* stability.aml  Harvey Greenberg/Dave Montgomery  UW
/* FOR LIMITED DISTRIBUTION Jul  1 17:38:39 PDT 2004
/* This aml calls to C program sloparea.  (It should be verion 1.087)
/* If you are not certain that you have the latest version of this package,
/* check http://gis.ess.washington.edu/stability/
/* The argument is the name of a FILLED arc/info elevation grid.
/* We shall attempt to automatically insure that we have
/* a well-behaved grid with no sinks and with no flat spots.
/* However, human intervention and quality assurance is desirable.
/* 3/22/96: changed reclass tables so that stability classes are numbered
/* numbered 1-7, not with color numbers.
/* 6/26/96: changed the position when the mask is invoked
/* 6/28/96: when elevations are read back in, change -1 to no_data; delete some files
/*          partial processing if sloped exists
/* 3/24/97: Fix bug in calculating critical rainfall when depth ne 1.
/*	    Better handling of missing maskgrid argument.
/*	    Add labels to the VAT files.
/* 4/10/97: If critical rainfall > rain-to-saturate, then uncond. stable
/*          CAVEAT: When you change soil depth, the integrated transmissivity
/*           does not change: large and small soil depths will have the same
/*           wetness.
/* 4/11/97: Added factor of sin into qsat
/* 5/7/97:  Fixed critcoh.  Calculate critical cohesion at saturation (interesting in itself)
/* 7/14/98:  Leave qsat as floating point to avoid rounding errors with a
/*		very small value of T.
/* 8/11/98: Automatic correction for grids with XY units of feet
/* 10/15/98: This is a special version skips the sink filling.
/*		It presumes that you have already filled sinks, presumedly
/*		with the ARC/INFO fill command.
/* 6/25/2004: Welcome to the new milleneum.  This update deals with the fact that
/*            most DEMs are no longer in integer meters, and that some people
/*	      have started to run ARC/INFO under windows.  In celebration of
/*	      Paul Allen having just opened a sci-fi museum and launched an
/*	      update of the X-15 rocketplane, we are breathing new
/*	      life into a venerable aml.  This should deal well with
/*	      floating-point grids and large flat spots.  It does require
/*	      that the input grid have no sinks.  It reads the xyunits and
/*	      z units from the grid file, and converts to meters if it
/*	      encounters FEET.  Arguments can override these factors.
/* 1/27/2015: Handle all input rasters the same:
/*		Apply a Z-offset and a scale factor to create an integer
/*		raster of 0 - 2147483647.  This gives the best range of values
/*		to sloparea. 
/*	      Also prevent some flat areas from becoming no-value.
/* 1/28/2014: Pull out all all archaic automatic Z-scaling code.  If the vertical
/*		and horizontal units do not match, the user must specify that.

&if [null %basin%] &then; &do
  &type The name of a filled elevation grid is a required argument
  &type The name of a mask grid is an optional second argument.
  &type '  If "grid" is the second argument, a maskgrid will be created'
  &type '  from a polygon cover named "boundary" in the current or parent directory.'
  &type '  If "compute" is the second argument, a watershed will be computed.'
  &type '  An automatically created watershed grid will be named shedg' using a grid named "mouth".
  &type '  Enter # for the mask grid if you want to skip it and enter more arguments.'
  &type The name of a directory containing reclass files is an optional third argument
  &type '  The default value is ./ i.e. the current directory.'
  &type The zfactor is the optional fourth argument.  If it is omitted, it will
  &type '  default to one.  If horizontal units are meters, use .3048 for feet.
  &type '  Use 0.1 if you have decimeters.'
  &type e.g. &run stability elv shedg ../
  &type 'Look at the aml text for mor documentaion.'
  &return stability.aml is ending before it begins.
&end
&if [show program] <> 'ARC' &then
  &return This aml must be run from the ARC prompt
display 0  /* so a grid canvas will not pop up
&if %maskgrid% eq # OR %maskgrid% eq ' ' &then
  &delvar maskgrid
&if [null %classdir%] &then
  &set classdir ./
&type Looking for reclass files on directory %classdir%
&if [type %zfactor%] ge 0 &then /* not a number
  &set zfactor = 1
&type zfactor is %zfactor%
&if ^ [exists %classdir%cohord.reclass -file] &then
  &return The files cohord.reclass and unstord.reclass must exist on directory %classdir%.
&if [var maskgrid] &then; &do
  &ty Masking data to %maskgrid%
  &if %maskgrid% eq 'grid' &then; &do
     &ty Atttempting the create shedg
     grid
     setwindow %basin%
     setcell %basin% 
     &if [exists shedg -grid] &then
       &return ERROR: a grid named shedg already exists
     &if [exists ../boundary -cover] &then
         shedg = polygrid(../boundary)
     &else; &do
         &if [exists boundary -cover] &then
            shedg = polygrid(boundary)
         &else
            &return ERROR: could find a boundary file from which to mask.
     &end
     &s maskgrid shedg
     QUIT /* back to arc
  &end
  &if %maskgrid% ne 'mouth' and ( ^ [exists %maskgrid% -grid] ) &then
    &return The grid %maskgrid% does not exist in this workspace.
&end
&if [exists sloped -grid] AND [exists flux -grid]  &then; &do
  &ty The grid "sloped" exists, so we are skipping some steps
  GRID
  &goto gotsloped
 &end
&if ^ [exists %basin% -grid] &then
  &return The grid %basin% does not exist in this workspace.
&describe %basin%
&goto continue
&label badunit
  &type The grid %basin% lists XY units of %prj$units% and Z units of %prj$zunits%.
  &type If you know that they actually match, enter a 1.
  &type If you know the the correction factor, enter it.
  &type e.g. If XY is meters and Z is feet, the factor is .3048
  &type Enter Q if you do not have sufficient information to assure that
  &type that you will get meaningful results.
  &s zfactor = [query 'zfactor?' ]
  &if [type %zfactor%] ge 0 &then /* not a number
     &return Exiting because of insufficient zfactor info.
&label continue
&s zoffset  %grd$zmin%
&s zrange = %grd$zmax% - %grd$zmin%
&if [exist tmpscalegrd -grid] &then ; kill tmpscalegrd
&s inflator [trunc [calc 2147483647 / %zrange%]]
GRID
 tmpscalegrd = int(( %basin% - %zoffset% ) * %inflator% + 0.5)
QUIT
&s ingrid tmpscalegrd
&type  min, max, range, factor, %grd$zmin% %grd$zmax% %zrange% %inflator%
/*&RETURN TESTING
&type Do not hold your breath.  This may take some time.
&type integer grid is %ingrid%
gridimage %ingrid% # filled BIL

/* This creates files filled.bil and filled.hdr
/* A BIL image can easily be manipulated by a C program, then converted
/* back into a grid.
/* The program sloparea starts (usually) by scaling everything
/* up by a factor of 5000.  Elevations are incremented one unit until
/* there are no flat spots.  Water is directed downhill not by the
/* TOPMODEL/ARCINFO "8-pointer" method, but by distributing flow to up to
/* eight neighbors proportional to dz/dy.
/* Output files include a 32-bit file of flux per contour unit and a
/* 32-bit file of "fixed" elevation.

&sys sloparea filled 1

/* copy the hdr files.  Now all bil files are 32-bit.
/* I got lazy, and assume you are running under a linux emulator.
&sys cp filled.hdr filledcum.hdr
&sys cp filled.hdr filledsloped.hdr
&sys cp filled.blw filledcum.blw
&sys cp filled.blw filledsloped.blw

&if  [exist tmpgrd -grid] &then ; kill tmpgrd all
imagegrid filledsloped.bil tmpgrd
imagegrid filledcum.bil flux

GRID   /* enter the grid module
/* until 6/26/96, we would set the mask grid here.  This had the advantage
/* of protecting us from DEM garbage outside the watershed (as in the Chehalis)
/* but has left us prone to edge effects.
  setmask %basin%
  sloped = float(tmpgrd) / %inflator% + %zoffset%
  kill tmpgrd all

&label gotsloped  /* for partial processing if the grid "sloped" exists
&if [var maskgrid] &then ; &do
 &if %maskgrid% eq 'mouth'  &then; &do
   &if ^ [exist mouth -grid] &then
      &Return ERROR: The grid named mouth does not exist here.
   &if  [exist shedg -grid] &then
      &Return ERROR: The grid named shedg already exists.
   flowdir = flowdirection(sloped)
   shedg = watershed(flowdir,mouth)
   &s maskgrid shedg
 &end
&end
slopedeg = slope(sloped,degree)
tanslope = tan(slopedeg div deg)
sinslope = sin(slopedeg div deg)
aobs = con( sinslope eq 0,123456789,flux / sinslope )
&if [var maskgrid] &then
  setmask %maskgrid%

/* set parameters
&s T = 65     /* integrated transmissivity in square meters / day
&set depth 1    /* meter soil depth
&s phi 33     /* angle of friction
&s bd = 2.0   /* bulk density
&set coh 2  /* kPascals
&set qtest = 100 /* steady-state rainfall for testing critical cohesion
&s tanphi [tan [angrad %phi%]]
&s low = [calc [calc [calc %bd% - 1] / %bd%] * %tanphi%]

wet1 = aobs / %T% / 1000  /* the thousand is a meter/mm thing.
/* cossq = sqr(sinslope / tanslope)
cossq = con(tanslope eq 0,1.,sqr(sinslope / tanslope))
/* OLD, why 0?  qsat = con(flux > 99999,0,int(%T% * 1000 / flux + 0.5))
qsat = %T% * 1000 / aobs
qsatc = reclass (qsat,%classdir%unstord.reclass) /* same reclass table as qust
&ty Rainfall-for-saturation is on qsat, and a reclassed version is on qsatc.
&do coh &list 2 4 8
 qnt%coh% = (1000 * %coh% /  ( 9807 * %depth% *  %tanphi% * cossq )  +  %bd%  * ( 1 - ( tanslope / %tanphi% )) ) * %T% * 1000 / aobs

/* OLD  qnt%coh% = (1000 * %coh% + %bd% * 9807 * cossq * (%tanphi% - tanslope)) / (9807 * %depth% * wet1 * %tanphi% * cossq)
  /* values for critical rainfall:
  /*   9999    stable at any cohesion if pore pressure equals zero
  /*   10000   stable when saturated if pore pressure equals zero
  /*   9000    steadystate rainfall >= 9000 mm/day
  /*   other:  steadystate rainfall
  qust%coh% = con(tanslope < %low%,9999,qsat < qnt%coh%,10000,qnt%coh% < 0,0,qnt%coh% > 9000,9000,qnt%coh%)
  kill qnt%coh% all
  qust%coh%c = reclass (qust%coh%,%classdir%unstord.reclass)
  &ty Rainfall-for-instability (with root strength = %coh% KN/mē) is
  &ty on qust%coh%, and a reclassed version is on qust%coh%c
&end

/* critical cohesion
/* We have eliminated the rounding to integer of critcoh%q% 5/7/97. 
/* This is not a big deal, as values are in Pascals, not kiloPascals
&set q = %qtest% /* wet1 has a factor built in, so that q is mm
cnt%q% = 9807 * %depth% * cossq * %tanphi% *  (%q% * wet1 - %bd% * (1 - tanslope / %tanphi%))
/* same as: cnt%q% = 9807 * %depth% * cossq * (%bd% * (tanslope - %tanphi%) + %q% * wet1 * %tanphi%)
cntsat = 9807 * %depth% * cossq * %tanphi% *  (qsat * wet1 - %bd% * (1 - tanslope / %tanphi%))
cohsat = con(cntsat < 0,0,cntsat > 90000,90000,cntsat)
kill cntsat
cohsatc = reclass (cohsat,%classdir%cohord.reclass)
&ty Critical cohesion at saturation is on cohsat,
&ty and a reclassed version is on cohsatc.
/* OLD:critcoh%q% = con(tanslope < %low%,0,cnt%q% < 0,0,cnt%q% > 90000,90000,int(cnt%q% + 0.5))
critcoh%q% = con(tanslope < %low%,0,cnt%q% < 0,0,cnt%q% > cohsat,cohsat,cnt%q% > 90000,~ 
90000,cnt%q%)
critcoh%q%c = reclass (critcoh%q%,%classdir%cohord.reclass)
&ty Critical cohesion (for a q of %qtest%) is on critcoh%q%,
&ty and a reclassed version is on critcoh%q%c.
&ty The temporary files wet1 and cossq can be killed.
&set coh = 0; &s phi = 45
&s tanphi [tan [angrad %phi%]]
&s low = [calc [calc [calc %bd% - 1] / %bd%] * %tanphi%]
qnt0 = %bd% * (1 - (tanslope / %tanphi%)) / wet1
qust45 = con(tanslope < %low%,9999,qsat < qnt%coh%,10000,qnt%coh% < 0,0,qnt%coh% > 9000,9000,qnt%coh%)
kill qnt0
qust45c = reclass (qust45,%classdir%unstord.reclass)
&type A qust45 is rainfall-for-instability with phi = 45°, and no root strength.
QUIT /* back to ARC
/* Here we try to avoid confusion by attaching labels to the VATs of
/* classified grids.  However, if %classdir%unstord.reclass and
/* %classdir%cohord.reclass are changed, the legends may not be valid.
tables
&do coh &list 2 4 8
  additem qust%coh%c.vat label 24 24 c
  sel qust%coh%c.vat
  resel value eq 1
  &if [show number select] ne 0 &then
    move 'Unconditionally\unstable' TO label
  nsel;resel value eq 2
  &if [show number select] ne 0 &then
    move '0 - 50 mm/day' TO label
  nsel;resel value eq 3
  &if [show number select] ne 0 &then
    move '50 - 100 mm/day' TO label
  nsel;resel value eq 4
  &if [show number select] ne 0 &then
    move '100 - 200 mm/day' TO label
  nsel;resel value eq 5
  &if [show number select] ne 0 &then
    move '200 - 400 mm/day' TO label
  nsel;resel value eq 6
  &if [show number select] ne 0 &then
    move ' > 400 mm/day' TO label
  nsel;resel value eq 7
  &if [show number select] ne 0 &then
    move 'Unconditionally\stable' TO label
  nsel;resel value eq 8
  &if [show number select] ne 0 &then
    move 'Stable at\this cohesion' TO label
&end 
additem critcoh%q%c.vat label 24 24 c
sel critcoh%q%c.vat
resel value eq 7
&if [show number select] ne 0 &then;    move 'Unconditionally\stable' TO label
  nsel;resel value eq 6
  &if [show number select] ne 0 &then;    move '0 - 1 kPascal' TO label
  nsel;resel value eq 5
  &if [show number select] ne 0 &then;    move '1 - 2 kPascal' TO label
  nsel;resel value eq 4
  &if [show number select] ne 0 &then;    move '2 - 5 kPascal' TO label
  nsel;resel value eq 3
  &if [show number select] ne 0 &then;    move '5 - 10 kPascal' TO label
  nsel;resel value eq 2
  &if [show number select] ne 0 &then;    move '> 10 kPascal' TO label
quit /* tables
&ty normal termination of the stability program.