I am starting with six rasters. One is basically a DEM (topographic elevation). The other five rasters represent water surface elevation for different recurrence intervals (10%, 4%, 2%, 1%, and 0.2%). The raster I want to create will represent the percent-annual-chance of flooding for each cell.
For example, I might have a topo elevation of 601 ft for one cell. The 10%, 4%, 2%, 1%, and 0.2%, water surface elevation rasters have values of 590, 595, 600, 605, and 615 ft, respectively. In order to find the percent-annual-chance of flooding for this cell, I need to graph a line using the known values. On the X-axis I put the natural log of the percents, and I put the elevation values on the Y-axis. I use numpy’s polyfit function to find the m and b of that line, and with a little algebra find that the percent-annual-chance of flooding at this cell is 1.9%.
I have written a few python scripts that can do this, but they only work for small datasets.
One works by turning the topo raster into points, extracting the WSEL values from those rasters to the points, and then running the calculation using UpdateCursor. This works great when my rasters are small and have low resolution, but it takes a VERY long time (non-feasible) to run for the size of a county.
Another script will turn the rasters into numpy arrays and interate through each row/column combination to perform the calculation. This runs much faster than using points, however it requires that all rasters are clipped to the exact same extent. Also, with the county-size extent and small cellsize (as small as 3 feet), my computer’s memory cannot handle turning all of the rasters into arrays. It will try, but eventually I get a memory error.
Has anyone else done raster processing like this? It would be great to be able to just use Raster Calculator / Map Algebra, but I haven’t been able to figure out how to use the numpy polyfit function with it.
I'm also attaching an image from my guidance document to help explain the calculation.