&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.