#! /usr/bin/env python

# All rasters must have the same cell size.  The first raster determines the cell size.
# If keyword arguments leftx=, boty=, rightx=, topy= are supplied,
#  they will be used.  Rasters can be trimmed to the bounds, but must not fall short.
# Otherwise, the first raster determines the bounds. 

def help():
  print """
   This implements the SHALSTAB slope stability model.

   This program can be downloaded from http://gis.ess.washington.edu/projects/shalstab_python
   It is based on stability.aml ( http://gis.ess.washington.edu/projects/stability/ )
   By rewriting with arcpy, we can now work directly from any arcmap-readable
   raster source, clipping files as we load.
   (At this time shalstab requires ESRI's licensed arcpy library, but we hope
    it will not be too hard to substitute public domain code.)
   We can now use datasets > 2gb and load everything as numpy arrays into memory on a
   computer with sufficient RAM.
   This program implements the SHALSTAB slope stability model, but most of the processing
   time involves calculating fuzzy flow accumulation for hillsides.
   It is sometimes used for this alone, as in http://gis.ess.washington.edu/areas/dryvalleys

   You must include one of the keyword arguments:
   a) DEM=pathname (not yet implemented.)
   b) FILLED=pathname (sink-filled DEM.  Integer or floating point.)
   c) ACC=pathname
   where "pathname" is the name of any raster
   (grid, tif, fgdb raster, etc.) recognized by arcmap and the arcpy libraries.
   The cell size of this raster must be shared by all others.
   If other keywords do not override it, the raster bounds
   will be the analysis window.
   Rasters may be in projected or geographic units.
   If the keyword is "DEM", it will be filled.  (not yet implemented.)
   If the keyword is "FILLED", fuzzy flow accumulation will be calculated.
     This is the only operation that takes significant CPU time.
     Fuzzy flow accumulation is not the same as standard ESRI flowaccumulation.
     See http://gis.ess.washington.edu/projects/stability/fuzziness
   If the keyword is "ACC", the program will proceed to stability calculations.
     It will still be necessary to use the FILLED keyword.
   If filled and acc rasters are created, shalstab will ask where to save them,
   unless overridden by other keywords.
   The sloped array is an intermediate product.  It can be saved, but we
   assume it is not really needed after the acc array is created.
   The prefix (default="") will be prepended to output rasters.
   If the prefix does not contain the string "gdb", the ".tif" will be appended.

  Default keyword arguments (Do not include quotes.  Avoid spaces.:
leftx      Left edge of analysis window (Default from input sink-filled raster.)
boty       Bottom edge ...
rightx     Right edge ...
topy       Top edge ...
rounding   Boolean.  Round coordinates to integers.  (Default False)
maxpix_str Either the name of a mask raster, or the maximum number of
            pixels to trace downhill.  If you stop processing before
            tracing fully, valley bottom contributing areas will be wrong,
            but you may save days of computer time. (Default 50000)
prefix     Prefix to output raster names.  (Default ""; Example: "e:/srtmsam/columb1/")
verbose    Verbose messages.  (Default=False)
acc        Name of input fuzzy flowaccumulation raster.  Use this to avoid recalculation.
filled     Nmae of an input sink-filled array.
xyeps      Leeway to avoid tripping up on rounding errors.  (Default 1e-7
diag       Diagnostic level for debugging.  Study the code to see what this does. (Default 0)
inslope    Input hillslope raster. (Default is "slopedeg").  If the name contains "tan",
            it will be used as tangent, and the program will calculate degrees.
            If the raster does not exist, slope will be calculated from the DEM.
filledout  Output sink-filled raster.  (Default=None)
accout     Output fuzzy flowaccumulation raster  (Default=None)
dirout     Output flowdirection raster. (Default=None)
slopedin   Input sloped raster. (Default=None)
slopedout  Output sloped DEM. (Default="sloped")
            For now, we need to convert slpray to a raster so that we can refill it,
            so we might as well save it.
dem        Input DEM  (Default=None)
deg2meters Interpret coordinates as degrees.  This overrides the implicit choice based on rightx - leftx.
maxits     Maximum number of iterations to tranverse the accumulating area.
             Default = nrows+ncols*1.5, which essentially means to do all cells.
T          Soil Transmissivity (square meters/day) (Default=65.)
depth      Soil depth in meters (Default=1.), or the name of a depth raster.
qtest      Rainfall in mm/day used to calculate critical cohesion  (Default=100.)

It is possible that the large number of keywords might interact in unfortunate ways.
No-data rows and columns will be trimmed (c. line 463).  
If everything is no-data, the program ends with a code of 2.

e.g. python -i y:/topog/python/shalstab.py FILLED=e:/srtmsam/southam1sec.gdb/filled accout=acc inslope=e:/srtmsam/southam1sec.gdb/tanslope prefix=e:/srtsam/columb1/ leftx=-76.940138888888888888 
   See also the "Programmer's Guide" and many comments embedded in the code.

  """
  raise Exception("Thanks for asking.")

# Programmer's Guide
# Functions preceding the body include:
## help():
## rowcol2world(rrow,ccol):    # rows are top down from 0
## world2rowcol(x,y):
## checkraster(inrast):
## fkeyval(keyword,default,arglist):
## skeyval(keyword,default,arglist):   # for strings
## bkeyval(keyword,default,arglist):   # for boolean

#* The raw elevation is demray
#* Flat areas are sloped.
#* Then the flowaccumulation is run.
# Intermediate rasters include dirray
#* Output files include:

# Other keyword arguments are accepted:
# No, this is implied if dy > 179.  Should degrees be converted to meters? (Default dag2meters=False)
# Bounds of analysis: leftx=, boty=, rightx=, topy= (Default determined by filled dem.)
# inslope=	Name of the input slope grid.  The default is to calculate
#               slope from the tweaked filled DEM.
# filledout=	Name of the filled DEM (with flat spots) created if the DEM=
#                keyword is used.  default=[dem]+"f"
# accout=	Name of the output flowaccumulation grid (default filled+"acc")
# dirout=	Name of output D* flowdirection (similar to ESRI flow direction)
#               default=None
# slopedout=	Name of the output tweaked elevation grid.  default=None
# maxits=	Maximum number sof times to traverse the array.
#                 default = int((nrows+ncols)*1.5).  A smaller value will cover the
#                 hillslopes and save time, but will miss parts of the streams.
# The program can operate either in the native DEM units or by converting latitude
# and longitude to meters.


from sys import argv,stdin,stdout
callname = argv[0]
version=callname+' version .13 Dec 6, 2018.' #  handling internal nodata bettr
print version
print argv
if len(argv) < 2:
  print len(argv) # diag
  help()
if argv[1] == 'help':
  print "help" # diag
  help()
import string,  os, time, types, pdb
from math import sin,cos,tan,radians,ceil
import arcpy
arcpy.CheckOutExtension("Spatial") # needed for slope
import numpy as np
radian = radians(1.)
# See https://gis.stackexchange.com/questions/75528/understanding-terms-in-length-of-degree-formula
# for using Taylor series for degree-to-meter on ellipsoid with m1-p3.
m1 = 111132.92     # latitude calculation term 1
m2 = -559.82       # latitude calculation term 2
m3 = 1.175         # latitude calculation term 3
m4 = -0.0023       # latitude calculation term 4
p1 = 111412.84     # longitude calculation term 1
p2 = -93.5         # longitude calculation term 2
p3 = 0.118         # longitude calculation term 3
start=time.time()
timelist=[]	# for tallying diagnostic times.
# So variables are implicitly global unless they are changed in a function?
try:
  uff = os.getenv("PATH").index("ArcGISx64")
  print "The path is set to run in 64-bit mode."
  uffda = os.getenv("PATH")[uff+9:]
  vnumb = uffda[:uffda.find(';')]	# Check vnumb after we check rows*cols
except:
  print "The PATH seems to be set for 32 bit",os.getenv("PATH")
  print "This progam might explode on big data."
  print r'Perhaps you want something like "export PATH=/cygdrive/c/Python27/ArcGISx6410.6:$PATH"'
  print r' or "/cygdrive/c/Program Files/ArcGIS/Pro/bin/Python/envs/arcgispro-py3'":$PATH'
  print 'Continue anyway? ("Y" or "y")'
  test = sys.stdin.readline().lower()
  try:
    uff = test.index('y')
    vnumb = '032-bit'
  except:
    raise Exception("export PATH=/cygdrive/c/Python27/ArcGISx6410.6:$PATH and restart.")

def rowcol2world(rrow,ccol):	# rows are top down from 0
  """ converts numpy row and col to world coordinates.
  """
  global cellsize, leftx, boty, rightx, topy, rounding, spatialRef
  return [leftx + (ccol+0.5)*cellsize,topy - (rrow+0.5)*cellsize]

def world2rowcol(x,y):
  """ converts numpy row and col to world coordinates.
  """
  global cellsize, leftx, boty, rightx, topy, rounding, spatialRef
  return [int((topy-y)/cellsize), int((x-leftx)/cellsize)]

def checkraster(inrast):
  """Checks any raster and extent and cellsize.
     If cellsize has been set, we check consistency, else we set it.
     If leftx has been set, we check it is not too far left, else we set it.
  """
  global cellsize, leftx, boty, rightx, topy, rounding, spatialRef
  print "Checking",inrast,"for cellsize and bounds.",cellsize #,leftx,boty,rightx,topy
  try:
    desc=arcpy.Describe(inrast)
  except:
    raise Exception("Raster "+inrast+" not found")
# bad idea #  arcpy.env.snapRaster = inrast	# added 9/14 to try to fix alignment.
  xtent=desc.Extent

  testy = desc.spatialReference
  testproj = testy.PCSName
  if spatialRef != None:
#   print "Comparing spatialReference"
    if testproj != spatialRef.PCSName:
      print "mismatch for",inrast, testproj,"v.",spatialRef.PCSName #No error.  Might be trivial datum difference
  else:
    spatialRef = testy  # If not previously defined.

  testy=desc.meanCellWidth
  if cellsize != None:		# or if 'cellsize' in locals():
    print "Comparing",testy,"with", cellsize
    if round(testy/cellsize,4) != 1:
      raise Exception("This has cellsize "+`testy`+",not "+`cellsize`)
  else:
    cellsize=testy	# if not previously defined # xdim in old program
  if diag==99:
    pdb.set_trace()

  testy=xtent.XMin
  if rounding:
    testy=round(testy)
  if leftx == None:
    leftx=testy
  else:
    if (testy-leftx) > xyeps:
      raise Exception(inrast+" does not extend to the left enough:"+`leftx`+" < "+`testy`)
  
  testy=xtent.XMax
  if rounding:
    testy=round(testy)
  if rightx == None:
      rightx=testy
  else:
    if (rightx-testy) > xyeps:	# Changing this to no-error.  For now, just for east, too
      print inrast+" does not extend to the east enough.  Changing rightx to",testy
      rightx= testy

  testy=xtent.YMin
  if rounding:
    testy=round(testy)
  if boty == None:
      boty=testy
  else:
    if (testy-boty) > xyeps:	# Changing this to no-error.  For now, just for south
      print inrast+" does not extend to the south enough.  Changing boty to",testy
      boty = testy

  testy=xtent.YMax
  if rounding:
    testy=round(testy)
  if topy == None:
      topy=testy
  else:
    if (topy-testy) > xyeps:
      print inrast+" does not extend to the north enough.  Changing topy to",testy
      topy=testy
  return True ## ["Returning "+inrast,cellsize,leftx]

def fkeyval(keyword,default,arglist):
  l=len(keyword)
  keyword=string.lower(keyword)
  for arg in arglist:
    if arg[:l]==keyword:
      return float(arg[l:]) # a bad argument will raise an exception
  return default

def skeyval(keyword,default,arglist):	# for strings
  l=len(keyword)
  keyword=string.lower(keyword)
  for arg in arglist:
    if string.lower(arg[:l])==keyword:
      return arg[l:]
  return default

def bkeyval(keyword,default,arglist):	# for boolean
  l=len(keyword)
  keyword=string.lower(keyword)
  for arg in arglist:
    if arg[:l]==keyword:
      if string.lower(arg[l:]) == "true":
        return True
      elif string.lower(arg[l:]) == "false":
        return False
      else:
        raise Exception("Bad boolean argument in",arg)
  return default
###################### end functions ##################################
arcpy.env.eompression="LZW"	# [raster].save respects this only for integers.
## Everything is imported as a numpy array

################ Define the default values for input ###################
cellsize = None
spatialRef = None
dem = None	# Need only if we have no filled or acc
filled = None	# Need only if we have no acc

################ Read arguments ## There might be a cleaner way #############
leftx      = fkeyval("leftx=",None,argv)
boty       = fkeyval("boty=",None,argv)
rightx     = fkeyval("rightx=",None,argv)
topy       = fkeyval("topy=",None,argv)
rounding   = bkeyval("rounding=",False,argv)	# The way rounding is implemented, it can shift output by one cell.
maxpix_str = skeyval("maxpix=","50000",argv)	# It could be mask raster
prefix     = skeyval("prefix=","",argv)	# Use this on simple file names; omit if name contains "/"
verbose    = bkeyval("verbose=",False,argv)
acc        = skeyval("ACC=","no",argv)
filled     = skeyval("FILLED=","no",argv)
xyeps      = fkeyval("xyeps=",1e-7,argv)	# was 0 # We need xyeps=1 to make Paraty line up.
diag       = fkeyval("diag=",0,argv)	# diagnostic level
inslope      = skeyval("inslope=","slopedeg",argv)
filledout    = skeyval("fillout=",None,argv)
accout       = skeyval("accout=",None,argv)
dirout       = skeyval("dirout=",None,argv)
slopedin    = skeyval("slopedin=","",argv) # because sloping can take days for a 2 GB array.
slopedout    = skeyval("slopedout=","sloped.tif",argv) # For now, we need to convert slpray to a raster so that we can refill it, so we might as well save it.
print "Arguments for leftx,boty,rightx,topy",leftx,boty,rightx,topy
lenprefix = len(prefix)
if lenprefix:
  prefixstripped = prefix.strip('/')
  if lenprefix > len(prefixstripped):	# output goes to a different directory.
    if not os.access(prefixstripped,os.F_OK ):	# directory does not exist.
      os.mkdir(prefixstripped)
    os.chdir(prefixstripped)
arcpy.env.workspace=os.getcwd()	# moved here march 19, 2018
# Unlike the Raster-to-Other_-Format tool, [raster].save will not compress floating-point rasters.
#  ==EXIST
if os.access("uncompressed",os.F_OK ):	# directory exist.
  print '"uncompreesed" directory exists'
else:
  print "Creating uncompressed in",arcpy.env.workspace
  os.mkdir("uncompressed")	# march 19, 2018
compresslist=""
if filled == "no":
  dem      = skeyval("DEM=","no",argv)	# otherwise None
  if dem == "no":
    raise Exception("No valid DEM or ACC raster")
  elif dem != None:
    checkraster(dem)
    minz = float(arcpy.GetRasterProperties_management (dem,"MINIMUM")[0]) # -1.
    nodata = round(minz-2)
    print nodata,"is the DEM nodata value."
else:
  checkraster(filled)
#minz does not see nodata, but is useful?
  minz = float(arcpy.GetRasterProperties_management (filled,"MINIMUM")[0]) # -1.
  nodata = round(minz-2)
  print nodata,"is the DEM nodata value."
## Depending on what the input file is:
if acc == "no":	# We will be creating fuzzy accumulation, and need to save it.
  if not accout:
    if filled:
      accout = filled+"facc"
    else:
      accout = dem+"facc"
else:
  checkraster(acc)	# first raster sets the bounds, unless set before.
  #minz = float(arcpy.GetRasterProperties_management (acc,"MINIMUM")[0]) -1.
  nodata = round(minz-2)
  print "minz set to",minz

nrows=int(round((topy-boty)/cellsize))
ncols=int(round((rightx-leftx)/cellsize))
if nrows*ncols >= 536870912:
  if vnumb < '10.5.1':	# I can do no better than read PATH, which does not distinguish 10.5 and 10.5.1
    print 'Warning: Rasters of more than 2GB will fail (without raising exceptions) with arcpy prior to 10.5.1.'
    print 'Understood?'
    uff = sys.stdin.readline()
  print "Rasters are larger than 2GB."
cleanup = nrows*ncols > 25123123	# arbitary cutoff.  Delete intermediate arrays if big.
deg2meters =  (topy-boty) < 179.
deg2meters = bkeyval("deg2meters=",deg2meters,argv)	# This is implicit in rightx - leftx.  This switch overrides
if deg2meters:
  print "degrees will be converted to meters. Rounding will be turned off."
  rounding = False
else:
  print "degrees will NOT be converted to meters."
maxits       = fkeyval("maxits=",int((nrows+ncols)*1.5),argv)
print "Actual        leftx,boty,rightx,topy",leftx,boty,rightx,topy
print "Rows and cols:", nrows,ncols
# Use code from memory.py to see if there is enough memory.
ll0 = arcpy.Point(leftx,boty)
ll0eps = arcpy.Point(leftx+xyeps,boty+xyeps)
######### shalstab model parameters.  Needs work to allow arrays.
phi=33.
for arg in argv:
  if arg[:4]=='phi=':
    try:
      phi=float(arg[4:]) # If it is not a number, it is a raster.
    except:
      rastername=arg[4:]
      checkraster(rastername)
      phi=arcpy.RasterToNumPyArray(rastername,ll0eps,ncols,nrows,nodata)
bd=2.0
for arg in argv:
  if arg[:3]=='bd=':
    try:
      bd=float(arg[3:]) # If it is not a number, it is a raster.
    except:
      rastername=arg[3:]
      checkraster(rastername)
      bd=arcpy.RasterToNumPyArray(rastername,ll0eps,ncols,nrows,nodata)
coh= -99.	# flag to use a list of 2, 4, 8
for arg in argv:
  if arg[:4]=='coh=':
    try:
      coh=float(arg[4:]) # If it is not a number, it is a raster.
    except:
      rastername=arg[4:]
      checkraster(rastername)
      coh=arcpy.RasterToNumPyArray(rastername,ll0eps,ncols,nrows,nodata)
T      = fkeyval("trans=",65.,argv)        # Soil transmissivity (square meters/day)
depth  = fkeyval("depth=",1.,argv)        # Soil depth
qtest    = fkeyval("qtest=",100.,argv)        # mm/day rainfall
print "Using model values: phi",phi,"; bd",bd,"; T",T,"; depth",depth,"; coh",coh,"; qtest",qtest
#print "Calling checkraster\n",checkraster.__doc__
area = cellsize*cellsize
edge = cellsize
print "saving diagnostics to",prefix+"python.txt"
t=open(prefix+"python.txt","w")
###if diag==7:	# testing optimization code
###  tt=open(prefix+"rowlists.txt","w")
t.write(version+'\n')
t.write("Using model values: phi"+`phi`+"; bd"+`bd`+"; T"+`T`+"; depth"+`depth`+"; coh"+`coh`+"; qtest"+`qtest`+"\n")
stdout.write("nrows,ncols: %i %i %i cells\n" % (nrows,ncols,nrows*ncols,))
t.write("nrows,ncols: %i %i %i cells\n" % (nrows,ncols,nrows*ncols,))

if deg2meters:
  stdout.write("%.2f square degrees\n" % (nrows*ncols*cellsize*cellsize))
  t.write("%.2f square degrees\n" % (nrows*ncols*cellsize*cellsize))
  edge = cellsize*111319.49 # y dimension
  area = edge*edge*math.cos(radian*(rowcol2world(nrows/2,0)[1])) # accurate for the middle.
if acc != "no":
  accray=arcpy.RasterToNumPyArray(acc,ll0eps,ncols,nrows,nodata)
# Do we need this?	  accray[accray==0]=eps
else:
  if filled == "no":
    print "filling DEM "
    if not filledout:
      filledout = "filled"	# 091
    arcpy.env.extent = arcpy.Extent(leftx,boty,rightx,topy)
    print nodata,"Creating the filled dem with extent",arcpy.env.extent.XMin,arcpy.env.extent.YMin
    outFill = arcpy.sa.Fill(dem)
    indesc = arcpy.Describe(dem)
    if indesc.pixelType[0] == 'F':	# floating point.  ESRI will not compress
      filled = "uncompressed/"+filledout
      compresslist+=";"+filled
    outFill.save(filled)
  else:
    print "Loading filled..."
  maxz = float(arcpy.GetRasterProperties_management (filled,"MAXIMUM")[0])
  eps = maxz/(2**24)	# (2**25 should work, but does not work on arrays)
  filray=arcpy.RasterToNumPyArray(filled,ll0eps,ncols,nrows,nodata)
# filray[filray==0.] = nodata
  print "New in version .0945, no longer converting filray=0. to nodata"
##### newer code: trimming no-data rows and cols
  oldmaxrow=nrows
  oldmaxcol=ncols
  for row in range(nrows):	# looking from north edge southward
    if filray[row].max() != nodata:
      break
  minrow=row
  if minrow == nrows-1:
    flagfile = prefix+"delete_empty_tile"	
    print "Input domain is totally nodata, flagging it with",flagfile
    f=open(flagfile,"w")
    import datetime
    f.write(datetime.datetime.fromtimestamp(time.time()).isoformat()+'\n')
    f.close()
#   raise Exception("There is no DEM data in the specified bounds.")
# fails outside function   return 98 
    sys.exit (0) # generates crashed-python popup? not in testy.py
#   uff=98
#   break
  for row in range(nrows-1,0,-1):
    if filray[row].max() != nodata:
      break
  maxrow=row+1
  if oldmaxrow-maxrow or minrow:
    stdout.write("Trimming to rows %d-%d\n" % (minrow,maxrow))
    filray=filray[minrow:maxrow]
    nrows=maxrow-minrow
  else:
    stdout.write("Processing all rows %d-%d\n" % (minrow,maxrow))
  mincol=0
  for col in range(ncols):
    if mincol:
      break	# break out of the outer loop if mincol has been found
    for row in range(nrows):
      if filray[row,col] != nodata:
        mincol=col+1	# We increment mincol so it will be a flag, even if it is zero.
        break
  mincol-=1
  maxcol=0
  for col in range(ncols-1,0,-1):
    if maxcol:
      break     # break out of the outer loop if maxcol has been found
    for row in range(nrows):
      if filray[row,col] != nodata:
        maxcol=col+1
        break
  if oldmaxcol-maxcol or mincol:
    stdout.write("Trimming to cols %d-%d\n" % (mincol,maxcol))
    filray=filray[:,mincol:maxcol]
    ncols=maxcol-mincol
  else:
    stdout.write("Processing all cols %d-%d\n" % (mincol,maxcol))
  if oldmaxrow-maxrow or mincol:	# trimming bottom or left
    leftx += mincol*cellsize
    boty += (oldmaxrow-maxrow)*cellsize
    topy -= minrow*cellsize
    ll0 = arcpy.Point(leftx,boty)
    ll0eps = arcpy.Point(leftx+xyeps,boty+xyeps)
##### end of newer code: trimming no-data rows and cols

  print "Sloping DEM and running fuzzy flow accumulation"
    #Try np.spacing() instead.  This can be f32, while np.nextafter is double
# print "maxz is",maxz,", eps[silon] is",eps
  stdout.write("maxz is %f, eps[silon] is %f\n" % (maxz,eps))
  t.write("maxz is %f, eps[silon] is %f\n" % (maxz,eps)) 
  if slopedin:
    print "Loading sloped array",slopedin
    slpray=arcpy.RasterToNumPyArray(slopedin,ll0eps,ncols,nrows,nodata)
  else:
    if filray.dtype == 'float32':
      slpray = filray.copy()
    else:
      print "Making a floating-point copy of the integer filled aray."
      slpray=np.array(nrows*ncols*[1.0], dtype='f32').reshape(nrows,ncols)
      for i in range(nrows):
        for j in range(ncols):
          slpray[i,j]=filray[i,j]

  print "calculating stability for..."	# 5/24/2016 kicked 247 lines to the right.
  nn=8
  nj = [1,1, 0,-1, -1,-1,  0, 1]# with rows numbered N->S, start with eastward cell and work clockwise
  ni = [0,1, 1, 1,  0,-1, -1,-1]
  nij=[[0,1],[1,1], [1,0],[1,-1], [0,-1],[-1,-1], [-1,0], [-1,1]]
  nb = [1,2,4,8,16,32,64,128]
  nbback = [16,32,64,128,1,2,4,8]
  nii=[[0,1],[1,1],[1,1], [1,0] , [0,-1], [-1,-1],[-1,-1],[-1,0]]
  njj=[[1,1],[1,0],[0,-1],[-1,-1],[-1,-1],[-1,0] ,[0,1] , [1,1]]
  dfactor = [1,1.4142,1,1.4142,1,1.4142,1,1.4142]	# It gets reset if degrees
# dir2dr&dir2dc were not used in shalstab before 10/18/2016, but had to be introduced from proffix after we caclulated flowdirection
  dir2dr=129*[0]	
  dir2dr[1] =   0
  dir2dr[2] =   -1
  dir2dr[4] =   -1
  dir2dr[8] =   -1
  dir2dr[16] =  0
  dir2dr[32] =  1
  dir2dr[64] =  1
  dir2dr[128] = 1
  dir2dc=129*[0]
  dir2dc[1] =   1
  dir2dc[2] =   1
  dir2dc[4] =   0
  dir2dc[8] =   -1
  dir2dc[16] =  -1
  dir2dc[32] =  -1
  dir2dc[64] =  0
  dir2dc[128] = 1
  eightpointdirs = 0
  # or myraster=arcpy.Raster(filled)
  #    maxz = myraster.maximum
  
  area = cellsize*cellsize
  edge = cellsize
  print "Total elapsed time:",time.time() - start
  print "Building accray"
  timelist.append("Building accray, tottime "+`time.time() - start`)
  if deg2meters:
    edge = cellsize*111319.49 # y dimension
    area = edge*edge*math.cos(radian*(rowcol2world(nrows/2,0)[1])) # accurate for the middle.
# This seems to take 9 hours for nrows,ncols = (36721, 18719)
#   accray=np.array([], dtype='f32')
#   for row in range(nrows):
#     edgex=edge*math.cos(math.radians(rowcol2world(row,0)[1])) # for tiago: 0.98727 to 0.98960
#     accray=np.append(accray,ncols*[edge*edgex]) 	# append is not a method of the array.
#   accray=accray.reshape(nrows,ncols)

# we ought to seed it.   accray=np.array(ncols*nrows*[0.], dtype='f32')
    accray=np.array(ncols*nrows*[0.], dtype='f32').reshape(nrows,ncols)
    for row in range(nrows):
      lat=radian*rowcol2world(row,0)[1]
#     edgex=edge*math.cos(math.radians(rowcol2world(row,0)[1])) #  not calibrated for ellipsoid
      edge = (m1 + (m2*cos(2*lat)) + (m3*cos(4*lat)) + (m4*cos(6*lat)))*cellsize # on ellipsoid
      edgex = ((p1*cos(lat)) + (p2*cos(3*lat)) + (p3*cos(5*lat)))*cellsize
      accray[row] =ncols*[edge*edgex] 	# append is not a method of the array.
  else:
    accray=np.array(nrows*ncols*[area], dtype='f32').reshape(nrows,ncols)
  dirray=np.array(nrows*ncols*[0], dtype='uint8').reshape(nrows,ncols)
  flagray=np.array(nrows*ncols*[False], dtype='bool').reshape(nrows,ncols)
# flagray=np.array(nrows*ncols*[0], dtype='uint8').reshape(nrows,ncols)	# This would take no more space.
  
  print "Total elapsed time:",time.time() - start
  print "Calculating directions on edges of filled DEM"
  timelist.append("Calculating directions on edges of filled DEM, tottime "+`time.time() - start`)
  # filray will have elevations and nodata 
  for i in range(nrows):
    for j in range(ncols):
      if i==0:
        flagray[i,j]=True
      elif j==(ncols-1):
        flagray[i,j]=True	
      elif i==(nrows-1):
        flagray[i,j]=True
      elif j==0:
        flagray[i,j]=True # else False
    
  ###
  # BEGIN HEART OF SLOPER 
  # Never pour to a cell that was higher in filray,
  #
  print "Total elapsed time:",time.time() - start
  timelist.append("sloping on the filled array tottime "+`time.time() - start`)
  
  rangenn = [5,6,7,4,3,0,1,2]	# better than range(nn) because it looks back at finished cells
  bigits = 0
  totfixed = 0
  fixes = 999
  todo=(nrows-2)*(ncols-2)
  itstart = time.time()
  rowlist = range(1,nrows-1)
  if diag==6:
    pdb.set_trace()
  if not slopedin:
    print "sloping on the filled array"
    while fixes > 0:
      fixes = 0
      bigits+=1
      for i in rowlist:	# skipping edge rows and columns
        for j in range(1,ncols-1):
          if not flagray[i][j]: 	# No assigned direction
            if filray[i,j] == nodata:
              flagray[i,j] =  True
              fixes+=1
              continue     # to next j   Done.
            elv = filray[i,j]
            for k in rangenn:	# Does this introduce a bias that leads to checkerboards?
              ii = i + ni[k]
              jj = j + nj[k]
              neighbfil = filray[ii,jj]
              if (neighbfil == nodata):
                flagray[i,j] =  True
                fixes+=1
                break	# to next j   Done.
              if neighbfil < elv:
                flagray[i,j] =  True
                fixes+=1
                break	# to next j   Done.
      if fixes == 0:
        break
      totfixed+=fixes
      print "checking it",bigits,"tweaked",fixes," total:",totfixed,"of",todo	#," time",time.time() - itstart
      print "flagged non-flat cells in",bigits,"its."," tottime",time.time() - itstart
#     pdb.set_trace()
    
    bigits = 0
  # totfixed = 0  6/16/2016.  See that most are done.
    fixes = 999
    maxjump = eps * 1000.	# A larger increment is probably an error.
#   print "maxjump is",maxjump
    stdout.write("maxjump is %f\n" % maxjump)
    t.write("maxjump is %f\n" % maxjump)
    while fixes > 0:
      itstart = time.time()	# for large grids, it might be wise to tally completed rows, as with flow accumulation.
      fixes = 0
      bigits+=1
      if diag==11:
        print "bigits is",bigits
        print slpray[115:121,59:61]
        pdb.set_trace()
      for i in rowlist:    # skipping edge rows and columns
        for j in range(1,ncols-1):
          if not flagray[i][j]:      # No assigned direction
            elv = filray[i,j]	# It would not have been incremented   # slpray[i,j]
  #          pours=[]
            lowest = 999999.
            for k in rangenn:
              ii = i + ni[k]
              jj = j + nj[k]
              if flagray[ii,jj]: # why do
  #              pours.append(k)
                neighbfil = filray[ii,jj]
                if neighbfil == elv:
                  if slpray[ii,jj] < lowest:
                    lowest = slpray[ii,jj]
                    lowii = ii
                    lowjj = jj
            if lowest < 999999.:
              slpray[i,j] = lowest + np.spacing(slpray[ii,jj])	# This insures an f32 input
              flagray[i,j] = True
              fixes+=1
      totfixed += fixes
##### Optimization code ########################
      if not bigits%15:
        donelist=[]
        if diag==2:
          pdb.set_trace()
        for rr in rowlist:
          done=True
          for cc in flagray[rr]:
            if not cc:      # anything in this row not done
              done=False
              break
          if done:
            donelist.append(rr)
        for rr in donelist:
          rowlist.remove(rr)
        if len(donelist):
          print "Deleting      ",len(donelist),"done rows,         ",len(rowlist),"  rows left" 
          timelist.append("Deleting done rows, leaving "+`len(rowlist)`)
        t.write("Incrementing it %4i tweaked %i,  time %.1f\n"%(bigits,fixes,time.time() - itstart))
        stdout.write("Incrementing it %4i tweaked %i,  time %.1f\n"%(bigits,fixes,time.time() - itstart))
        rowlist.reverse()
      timelist.append("it "+`bigits`+ " time "+`time.time() - itstart`)
      if fixes==0:	# ready to drop out of the while loop
        print "Are we done?"
  #     pdb.set_trace()
    stillbad = 0
    for i in range(1,nrows-1):
      for j in range(1,ncols-1):
        if not flagray[i][j]:
          stillbad += 1
    #     print "stillbad",i,j
    if stillbad:
      print stillbad,"False records still not fixed."
      pdb.set_trace()
  ################################ july 13.  Can a post-tweaker fix slpray?###########
  ################################ October 14.  No it cannot.  Just fill the array.
    if slopedout:
#     if prefix:
#       if slopedout.find("/") < 0:	# not a pathname
# no need!slopedout = prefix+slopedout
      print "Saving sloped DEM as",slopedout," No special adjustments.  Nodata is",nodata,"compression is",arcpy.env.compression
      slopedrast = arcpy.NumPyArrayToRaster(slpray,ll0,cellsize,cellsize,nodata)
      slopedrast.save("uncompressed/"+slopedout)
      compresslist+=";uncompressed/"+slopedout
      arcpy.DefineProjection_management("uncompressed/"+slopedout, spatialRef)
      arcpy.BuildPyramids_management("uncompressed/"+slopedout)
### a good time for calculate raster domain?  Or do it in python?
  ################### It might be worth replacing the next 5 lines. or not.
    fillrast = arcpy.sa.Fill("uncompressed/"+slopedout) # python fill, as in http://gis.stackexchange.com/questions/211568/numpy-flood-fill-dem-at-specific-seed-location may not be mature.
    try:
      fillrast.save("uncompressed/refilled.tif")	# I think this should work > 2GB
      compresslist+=";uncompressed/refilled.tif"
    except:
      print 'Failed to replace existing', prefix+"uncompressed/refilled.tif"
    arcpy.DefineProjection_management(prefix+"uncompressed/refilled.tif", spatialRef)
    
    # END HEART OF SLOPER  Now calculating fuzzy flow accumulation. */
  print "Elapsed tottime:",time.time() - start
  if diag==6:
    pdb.set_trace()
  print "Now add up the areas"
  timelist.append("Now add up the areas, tottime "+`time.time() - start`)
  if slopedin:
    fillrast = slopedin
    flagray=np.array(nrows*ncols*[True], dtype='bool').reshape(nrows,ncols) # not set in sloping operation
  dirrast = arcpy.sa.FlowDirection(fillrast)	############################### HUHHHHHHHHHHH ##
  dirray=arcpy.RasterToNumPyArray(dirrast,ll0eps,ncols,nrows,nodata)
  
  for i in range(0,nrows):
    flagray[i,0] = False		# We never run flows off the edge cells
    flagray[i,ncols-1] = False		# We never run flows off the edge cells
    for j in range(0,ncols):
      if (i==0) or (i==nrows-1):
        flagray[i,j] = False          # We never run flows off the edge cells
#     if slpray[i,j] == nodata:
      if filray[i,j] == nodata:
        if diag==667:
          print "filray is nodata",i,j
          pdb.set_trace()
        accray[i][j] = 0.	#eps	# 11/6/2017  0.	# else area;  was fray
        flagray[i,j] = False	# no need to be done
  
# pdb.set_trace()
  
  totalcount = (nrows - 2) * (ncols - 2)
  print "Calculating flow for",totalcount,"cells."
  rowlist = range(1,nrows-1)	# For optimizing
  bigits = 0
  thiscount = 99
  rangenn = [5,6,7,4,3,0,1,2]	# in case in ended the last loop reversed.
###  maxlist = []	# list of rows with highest elevations.  for optimization.
  while thiscount > 0:
    itstart = time.time()
    bigits+=1
### if diag==7:	# testing optimization code
###   print bigits,thiscount,len(rowlist)
#     pdb.set_trace()
    thiscount = 0
    for i in rowlist:	#range(1,nrows-1):		# skipping the edge cells.
      if deg2meters:
        lat = radian*(i * cellsize + boty)
        edge = (m1 + (m2*cos(2*lat)) + (m3*cos(4*lat)) + (m4*cos(6*lat)))*cellsize # on ellipsoid
        edgex = ((p1*cos(lat)) + (p2*cos(3*lat)) + (p3*cos(5*lat)))*cellsize
        diagonal = (edge**2 + edgex**2) ** .5
        dfactor = [edgex,diagonal,edge,diagonal,edgex,diagonal,edge,diagonal]
      if diag == 99:
        pdb.set_trace()
###   if bigits==1:
###     maxlist.append([slpray[i].max(),i])	#for optimization
      for j in range(1,ncols-1):
        if flagray[i,j]:	# Needs to be done
          elv = slpray[i,j]
          ready = True
          for k in rangenn:
            ii = i + ni[k]
            jj = j + nj[k]
#           if (ii < 0) or (jj < 0) or (ii == nrows) or (jj == ncols):
#             pdb.set_trace()		# why a problem?
            if (slpray[ii,jj] > elv) and flagray[ii,jj]:
              ready = False
              break  # not ready for this cell */
          if ready:  # no unaccounted-for uphill flows */
            if diag == 12:
              pdb.set_trace()
            thiscount+=1
            tweight = 0.	# 
            cum = accray[i,j]
            maxdrop = -9999888
            pourk = -1
            weights=[]
            ks=[]
            for kk in rangenn:	# We could reuse k, but this is okay.
              iii = i + ni[kk]
              jjj = j + nj[kk]
              if (iii < 0) or (jjj < 0) or (iii == nrows) or (jjj == ncols):
                pdb.set_trace()	# does not happen
              if elv > slpray[iii,jjj]:
                weight=(elv + eps - slpray[iii,jjj])/dfactor[kk]
                weights.append(weight)
                ks.append(kk)
                tweight += weights[-1]
                if weight > maxdrop:
                  maxdrop = weight	# need this?
                  pourk = kk
            nks=len(ks)
            if diag == 14:
              pdb.set_trace() # why neg?
            if nks:	# pours at least one place
              for indx in range(len(ks)):	# all the lower cells
                kk = ks[indx]	# 1-8 k values
                iii = i + ni[kk]
                jjj = j + nj[kk]
                if (iii < 0) or (jjj < 0) or (iii == nrows) or (jjj == ncols):
                  pdb.set_trace()	# does not happen
                accray[iii,jjj] += cum * weights[indx]/tweight 
              dirray[i][j] = nb[pourk]
            else:
# Just pour as directed by esri # next 6 lines: Octo 14, 2016
              eightpointdirs+=1		# Count problem cells for diagnostics
              kk=dirray[i,j]
              iii = i + dir2dr[kk]
              jjj = j + dir2dc[kk]
              if (iii < 0) or (jjj < 0) or (iii == nrows) or (jjj == ncols): #off the edge. okay
                continue
              accray[iii,jjj] += cum
  
            flagray[i][j] = False	 # -bigits /* flag this cell done */
    totalcount = totalcount - thiscount
#   print " it",bigits,thiscount,"more done,",totalcount,"to go. time",time.time() - itstart
### if bigits==1:
###   maxlist.sort()
###   maxlist.reverse() 
###   optimize_mode=True
### elif bigits < nrows/25:
###   if bigits%10 == 5:
###     rowlist = []
###     for ilist in range(100 + bigits*8):	# a tunable parameter
###       rowlist.append(maxlist[ilist][1])
###     rowlist.sort()
###     timelist.append("Setting rowlist to "+`len(rowlist)`+" rows")
###     if diag==7:	# testing optimization code
###       pdb.set_trace()
###       tt.write("\nit %i\n" % bigits)
###       for indx in rowlist:
###         tt.write("%i " % indx)
### else:	# Through with the fancy stuff
###   if optimize_mode:
###     optimize_mode=False
###     rowlist = range(1,nrows-1)
###     timelist.append("Setting rowlist to "+`len(rowlist)`+" rows")
    
    if not bigits%20:
###   stdout.write(" it %i, %d cells done in %i rows. time %.1f seconds\n" % (bigits,thiscount,len(rowlist),time.time() - itstart))
########### optimization by skipping done rows ########################
      stdout.write(" it %i, %i more done.  time %.3f\n" % (bigits,thiscount,time.time() - itstart))
      if bigits > nrows/16:
        donelist=[]
        if diag==2:
          pdb.set_trace()
        for rr in rowlist:
          done=True
          for cc in flagray[rr]:
            if cc:	# anything in this row not done
              done=False
              break
          if done:
            donelist.append(rr)
        for rr in donelist:
          rowlist.remove(rr)         
        if len(donelist):
          print "Deleting      ",len(donelist),"done rows,         ",len(rowlist),"  rows left" 
          timelist.append("Deleting "+`len(donelist)`+" done rows, leaving "+`len(rowlist)`)
        rowlist.reverse()	# after fri2, try this every tenth, and start later
    timelist.append("it "+`bigits`+" "+`thiscount`+" more done, "+`totalcount`+" to go. time "+`time.time() - itstart`)
    if bigits > maxits:
      maxval = ceil(accray.max())
      accray[flagray]=maxval	# Flag unprocessed cells with the largest value.
      print "Max iterations exceeded; leaving with",totalcount,"unprocessed cells, which are flagged as",maxval
      break
  for lin in timelist:
    t.write(lin+'\n')	# python.txt
  
  print "Elapsed time:",time.time() - start, "Saving grids"
# pdb.set_trace()
  
  #raise Exception("Finished adding areas.  Leaving for now")
  if accout:
#   if prefix:
#     if accout.find("/") < 0:	# not a pathname
#       accout = prefix+accout
    if (accout.find('gdb') < 0) and (accout[-4:] != '.tif'):
      accout = accout+'.tif'
    print "Saving flowaccumulation raster as",accout
    accrast = arcpy.NumPyArrayToRaster(accray,ll0,cellsize,cellsize,nodata)
    accrast.save("uncompressed/"+accout)
    compresslist+=";uncompressed/"+accout
    try:
      arcpy.DefineProjection_management("uncompressed/"+accout, spatialRef)
    except:
      print "ERROR! Could not save","uncompressed/"+accout
  if dirout:
    if prefix:
      if dirout.find("/") < 0:	# not a pathname
        dirout = prefix+dirout
    print "Saving flowdirection raster as",dirout
    dirrast = arcpy.NumPyArrayToRaster(dirray,ll0,cellsize,cellsize,nodata)
    dirrast.save(dirout)	# integer.  This should compress.
    arcpy.DefineProjection_management(dirout, spatialRef)
    arcpy.BuildPyramids_management(dirout)
  if cleanup:
    del dirray
  print eightpointdirs,"cybergenetic sinks fixed by using ESRI flowdirections."
print "Using model values: phi",phi,"; bd",bd,"; T",T,"; depth",depth,"; coh",coh,"; qtest",qtest
tanphi = np.tan(radian*phi)
low = (1. - (1./bd)) * tanphi		# Stable at any cohesion.
##slopedeg = slope(sloped,degree)	# insert the goodies from myslope.py
##slope    = skeyval("slope=",None,argv)	# slprast in first threefoot test
##if slope:	# sloped DEM
##  checkraster(slope)
##  slpray=mknumpray(slope,ll0eps,0,slyce4)
##else:
# pdb.set_trace()
try:
  del filray
except:	# We skipped the creation of this array.
  uff=1
try:	# hmm, "try" does not intercept errors in the next line.
# checkraster("slopedeg")
  checkraster(inslope)
  print "Assuming that an existing slope grid is what we want"
  if inslope.find("deg") >= 0:
    print 'degrees'
    degray = arcpy.RasterToNumPyArray(inslope,ll0eps,ncols,nrows,0)
    degray[degray==0]=eps
    radray = degray*radian
    tanray = np.tan(radray)
    if cleanup:
      del degray
  elif inslope.find("tan") >= 0:
    print 'We are starting with the tangent'
    tanray = arcpy.RasterToNumPyArray(inslope,ll0eps,ncols,nrows,0)
    tanray[tanray==0]=eps
    radray = np.arctan(tanray)	
    print "Got radray from tanray"
  else:
    print "Could not indentify",inslope,"as degrees or tangent or what!"
    pdb.set_trace()
except:
  if deg2meters:
    print "Calculate slopes with the right latitude factor."
    tanray = np.array(nrows*ncols*[0.],dtype='f32').reshape(nrows,ncols)
    xmatrix = np.array([1,0,-1, 2,0,-2, 1,0,-1],dtype='f32').reshape(3,3)	# fixed 10/24/17
    ymatrix = np.array([1,2,1, 0,0,0, -1,-2,-1],dtype='f32').reshape(3,3) 
#    pdb.set_trace()	# something odd happens here.
    for i in range(1,nrows-1):
      lat = radian*(i * cellsize + boty)
      edge = (m1 + (m2*cos(2*lat)) + (m3*cos(4*lat)) + (m4*cos(6*lat)))*cellsize # on ellipsoid
      edgex = ((p1*cos(lat)) + (p2*cos(3*lat)) + (p3*cos(5*lat)))*cellsize
      for j in range(1,ncols-1):
        mx = ((xmatrix * slpray[i-1:i+2,j-1:j+2]).sum() / (8.*edgex))**2.0
        my = ((ymatrix * slpray[i-1:i+2,j-1:j+2]).sum() / (8.*edge))**2.0 # corrected error from .0951
        tanray[i,j] = (mx+my)**0.5	# fixed 4 lines 10/24/17

# something odd here    pdb.set_trace()

    radray = np.arctan(tanray)	
    tanrast = arcpy.NumPyArrayToRaster(tanray,ll0,cellsize,cellsize,nodata)
    tanrast.save(prefix+"uncompressed/tangent.tif")
# QA check 10/25 shows this good, but it is built off sloped, not the raw DEM
    arcpy.DefineProjection_management("uncompressed/tangent.tif", spatialRef)
    compresslist+=";uncompressed/tangent.tif"
#   raise Exception("goodies from myslope.py not inserted yet.  pause")
  else:
    try:
      arcpy.gp.Slope_sa(filled,"uncompressed/slopedeg.tif","DEGREE") # or should I do it myself
      compresslist+=";uncompressed/slopedeg.tif"
      degray =arcpy.RasterToNumPyArray("uncompressed/slopedeg.tif",ll0eps,ncols,nrows,0)
      degray[degray==0]=eps
      radray = degray*radian
      tanray = np.tan(radray)
    except:
      print "Slope_sa says that it failed, but did it?"
      pdb.set_trace()
#if uff!=98:	 # uff would be set before "break" on an empty matrix.
if cleanup:
  try:
    del slpray
  except:	# We skipped the creation of this array.
    uff=2
sinray = np.sin(radray)
##aobs = con( sinslope eq 0,123456789,flux / sinslope )
smallest = np.min(sinray)
if not smallest:
  print "There is a zero slope value.  We will revert to the 8-pointer method there. Really?!!"
sinray[sinray==0]=eps
#  No latitude correction here for edge. Not knowing the aspect, the edge length is problematic
aobsray = accray / edge / sinray	
if cleanup:
  del accray
cossqray = np.cos(radray)**2
##qsat = %T% * 1000 / aobs
qsatray = (T * 1000.) / aobsray
#qustray=np.array(nrows*ncols*[-666.0], dtype='f32').reshape(nrows,ncols)
if type(coh) != types.StringType:
  if coh < 0.:	# flag
    cohlist = [2,4,8]
  else:
    cohlist = [coh]	# one value
else:
  cohlist = [coh]		# one array
  outname='qust'+`coh`+'.tif'	# Does this name get overridden?

for coh in cohlist:
  qntcoh = (1000. * coh /  ( 9807 * depth *  tanphi * cossqray )  +  bd  * (1 - ( tanray / tanphi )) ) * T * 1000 / aobsray
##  qnt%coh% = (1000 * %coh% /  ( 9807 * %depth% *  %tanphi% * cossq )  +  %bd%  * ( 1 - ( tanslope / %tanphi% )) ) * %T% * 1000 / aobs
# 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
########################### New flags (See e:/tiago/qust2.lyr) :
#  9000    steadystate rainfall >= 9000 mm/day
#  9100    stable at this cohesion when saturated, if pore pressure equals zero (10000 in aml)
#  9200    stable at any cohesion when saturated, if pore pressure equals zero (9999 in aml)
#  (-99    no-data)
  qustray=np.array(nrows*ncols*[-99], dtype='int16').reshape(nrows,ncols) 
  for i in range(nrows):
    for j in range(ncols):
      if tanray[i,j] < low:
        qustray[i,j]=9200
      elif qsatray[i,j] < qntcoh[i,j]:
        qustray[i,j]=9100
      elif qntcoh[i,j] < 0.:
        qustray[i,j]=0
      elif qntcoh[i,j] > 9000.:
        qustray[i,j]=9000
      else:
        qustray[i,j]=qntcoh[i,j] + .4999	# round to integer

# qustray = qntcoh.copy()
# qustray[tanray<low]=9200.
# qustray[qsatray < qntcoh]=9100	# moved this before 9200 line 3/21/2017, moved back 11/8/17
# qustray[qntcoh < 0.]=0.
# qustray[qntcoh > 9000.]=9000.
  outname=prefix+"critQ"+`coh`+".tif"

  tmprast = arcpy.NumPyArrayToRaster(qustray,ll0,cellsize,cellsize,nodata)
  tmprast.save(outname)	# integer
  arcpy.DefineProjection_management(outname, spatialRef)
# And critical cohesion
# Check this part, especially.
#cntsat = 9807 * depth * cossqray * tanphi *  (qsat * wet1 - bd * (1 - tanray / tanphi))
#cohsat =     con(                   cntsat < 0,0,cntsat > 90000,90000,cntsat)
#critcoh%q% = con(tanslope < %low%,0,cnt%q% < 0,0,cnt%q% > cohsat,cohsat,cnt%q% > 90000,90000,cnt%q%)
# moved the next line to keep things in order 12/14/2016
wet1 = aobsray / (T * 1000.)
cnt = 9807 * depth * cossqray * tanphi *  (qtest * wet1 - bd * (1 - tanray / tanphi))	# cnt%q% in aml
cntsat = 9807 * depth * cossqray * tanphi *  (T*1000./aobsray * wet1 - bd * (1 - tanray / tanphi))	# cnt%q% in aml
if cleanup:
  del cossqray
cohsat = np.where(cntsat < 0,0,cntsat)
cohsat[cohsat > 9000]=9000
tmprast = arcpy.NumPyArrayToRaster(cohsat,ll0,cellsize,cellsize,nodata) # pulled out indent 10/6/2017.  right?
outname=prefix+"uncompressed/cohsat.tif"
compresslist+=";uncompressed/cohsat.tif"
tmprast.save(outname)
arcpy.DefineProjection_management(outname, spatialRef)
critcoh = np.where(tanray < low,0,cnt)
critcoh[cnt < 0]=0
critcoh=np.where(cnt > cohsat,cohsat,critcoh)
critcoh[cnt > 9000]=9000	# still not matching aml on nodata and on flags.  Get back to it.
tmprast = arcpy.NumPyArrayToRaster(critcoh,ll0,cellsize,cellsize,nodata)
outname=prefix+"uncompressed/critcoh.tif"
tmprast.save("uncompressed/critcoh.tif")
compresslist+=";uncompressed/critcoh.tif"
tmprast.save(outname)
arcpy.DefineProjection_management(outname, spatialRef)
# Now compressed all the uncompressed output files
compresslist.strip(";")
print "Compressing the files", compresslist
arcpy.RasterToOtherFormat_conversion(Input_Rasters=compresslist, Output_Workspace=arcpy.env.workspace, Raster_Format="TIFF")
if diag==666:
  print "Odd errors follow; pause to check it out."
  pdb.set_trace()
for bigfile in compresslist.split(';')[1:]:
  big=os.stat(bigfile).st_size
  small=float(os.stat(bigfile[13:]).st_size)
  compression=small/big
  stdout.write('%s compressed to %.3f size\n' % (bigfile[13:],compression))
  t.write('%s compressed to %.3f size\n' % (bigfile[13:],compression))
import shutil
#try:
#  shutil.rmtree('uncompressed')
#  t.write('"uncompressed" was deleted\n')
#  stdout.write('"uncompressed" was deleted\n')
#except:
#  t.write('"uncompressed" was not removed because of open files.')
#  stdout.write('We presume that "uncompressed" will be removed when python closes.')
print 'We have stopped trying to delete the "uncompressed" directory because we think that attempts to delete locked files cause Error popups that freeze scripts.'
end = time.time()
elapsed = end - start
stdout.write("Elapsed time was %.0f seconds.\007\a" % elapsed)
t.write("Elapsed time was %.0f seconds.\a" % elapsed)

t.close()	# python.txt diagnostics
###if diag==7:	# testing optimization code
###  tt.close()	# rowlist.txt diagnostics
print "!!!!!!!!!!!!!! Have a gneiss data!!!!!!!!!!!\n\n\n"