The program holds arrays in memory for (some) speed, and tile size is limited only by the RAM size of the machine on whcih it is run. (It needs to have code inserted for monitoring a freeing memory.)
The program can read a wide variety of raster formats. They need to have matching projections and cell size, but do not need matching spatial extent. shalstab.py now depends on the proprietary ESRI arcpy library, but we hope to rewrite it for open source libraries.
The routing of groundwater on hillsides requires replacement of 8-pointer flow accumulation with what we insist on calling fuzzy flow accumulation. (See illustrations.)
GIS technicians typically fill spurious "sinks" on a DEM, then run an 8-pointer flowdirection algorithm. When the algorithm encounters a flat area, it finds a good-enough path to the pour point. Fuzzy flowaccumulation does not use a flowdirection raster. It must start at the top and distribute flow proportionally to all downhill cells. Then it iterates through the array again, dealing with all the cells that have no upstream contributors. This requires an array with none of the flat spots that are created by filling sinks.
Flat spots have been dealt with by preconditioning the DEM with small slopes. In the 20th-century version of shalstab, flat spots were incrementing upstream from their pour points, one millimeter at a time. When 30m DEM were supplanted by 10m DEMs, we started to find flat spots greater then 1000 cells long, so cells could sometimes be incremented higher than the neighbor cells that flowed into them, causing errors. This problem was reduced by scaling up the (integer) elevations to cover the range of 32-bit integers.
Now that we have floating-point lidar DEMs, we increment flat spots by the smallest possible epsilon, but we still have trouble. As we increment flat cells by hundreds of epsilons, we may encounter an occasional cell that is a millimeter higher than the flat spot, and thus becomes a sink. I try to iterate through the DEM until such minor artifacts are fixed, but I have failed to build a robust algorithm that does not sometimes undertake infinite repairs. That is what happened to Oregon City. This problem can be reduced by using mapped or guessed culverts, but that does not guarantee reliable success.
To build a robust fuzzy flowaccumulation algorithm, I think I should: 1) Make one iteration in adding slopes to flat spots. 2) Fill the sinks that I have accidentally created. 3) Let ESRI calculate eight-pointer flow directions on this DEM. 4) Calculate fuzzy accumulation, but If I find a cell that has no lower neighbors, use the ESRI flowdirection.This is the version that continues to fail on complex DEMs, and which may or may not get debugged:
shalstab.py looks okay, but is not yet certified. And it contains debugging code.
Oregon City example needs to be redone with the current version.
This is a premature draft page. If you have questions, ask Harvey