Substitute different equation parameters (e.g. slope) when batch processing raster calculator?

[JacobNederend]

This is a cross-post from the GIS StackExchange because it gained no traction there.

I have 15 multiband rasters for a total of 45 single-band images (5 dates x 3 sensors x 3 bands per sensor for which I need to apply a unique calibration equation. All the equations follow the format Y=A*exp(B*X), and I have a spreadsheet with all the corresponding A & B parameters for each.

Is there a feasible way for me to apply these functions iteratively using ModelBuilder to speed up the process? I am also willing to learn some ArcPy to make it work. Essentially, how would I specify the A * B variables in the calculator, and then tell it to search the spreadsheet to make the right substitutions? I couldn't turn up any forum or support site posts about this.

My current plan is:

1. Rename all the rasters to something logical like: "sensor_date.tif"
2. Split into temporary single-bands using "Make Raster Tool"
3. Feed output into raster calculator and output calibrated singlebands as .tif
4. Feed new singlebands into raster calculators to find several different indices (NDVI, NDRE, SAVI, etc.)

[DanPatterson_Retired]

could be done in the raster calculator, but use RasterToNumPyArray and use numpy to do the heavy lifting

>>> X = np.random.randint(1,5, size=(3,5,5))
>>> A = 1.1
>>> B = 0.5
>>> Y = A*np.exp(B*X)
>>>
>>> Y
array([[[ 8.1,  3.0,  3.0,  3.0,  8.1],
        [ 1.8,  3.0,  8.1,  1.8,  3.0],
        [ 1.8,  4.9,  1.8,  4.9,  3.0],
        [ 1.8,  8.1,  1.8,  4.9,  8.1],
        [ 1.8,  8.1,  4.9,  4.9,  1.8]],

       [[ 3.0,  1.8,  1.8,  4.9,  1.8],
        [ 1.8,  4.9,  4.9,  8.1,  3.0],
        [ 1.8,  8.1,  3.0,  3.0,  1.8],
        [ 1.8,  1.8,  3.0,  3.0,  3.0],
        [ 8.1,  8.1,  4.9,  3.0,  1.8]],

       [[ 4.9,  4.9,  8.1,  3.0,  8.1],
        [ 8.1,  3.0,  4.9,  3.0,  3.0],
        [ 8.1,  4.9,  3.0,  1.8,  3.0],
        [ 3.0,  1.8,  1.8,  1.8,  3.0],
        [ 3.0,  1.8,  3.0,  8.1,  1.8]]])
>>>
>>> Y[0,...]
array([[ 8.1,  3.0,  3.0,  3.0,  8.1],
       [ 1.8,  3.0,  8.1,  1.8,  3.0],
       [ 1.8,  4.9,  1.8,  4.9,  3.0],
       [ 1.8,  8.1,  1.8,  4.9,  8.1],
       [ 1.8,  8.1,  4.9,  4.9,  1.8]])
>>> # Slice to get the individual band results and use NumPyArrayToRaster
>>> # to get out to your favorite grid format if needed.
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

In the example above, the sequence is bands (3) X and Y direction shape (5, 5) So the above would represent a 3 band raster with a 5x5 extent.  If your raster has the bands swapped, (ie 5,5,3) then you may have some swapping to do, or just slice them out individually and process.

An yes, It took about a minute to produce the results, so you could make a list of A, B coefficients and process within a function.

[JacobNederend]

Thanks for the suggestion Dan. I like this idea, just need to learn lots more about python and how to make it happen.

[DanPatterson_Retired]

See my blog for examples of working with numpy

as for the 'ins' and 'outs' to arcmap/pro

RasterToNumPyArray—Help | ArcGIS Desktop 

NumPyArrayToRaster—Help | ArcGIS Desktop 

[JacobNederend]

import arcpy import os
arcpy.env.workspace = "E:\\ArcGIS_Work\\ModelBuilder_Testing\\ModelBuilder_NIR"
newdir = "E:\\ArcGIS_Work\\ModelBuilder_Testing\\mb_singlebands"
# Create a list of rasters in the workspace
rasters = arcpy.ListRasters("*","TIF")
list1 = [] print rasters
# Copy band_1 of each raster in list to create temporary singleband image. Keep same name but use basename property to exclude file extension
for i in rasters:
    print 
   strip_ext = arcpy.Describe(i)
   temp_rast = arcpy.MakeRasterLayer_management(i, os.path.join(newdir, strip_ext.basename+"_b1"),"#", "#", "1") # Makes temporary raster of band 1
   b1_rast = arcpy.CopyRaster_management(temp_rast, os.path.join(newdir, strip_ext.basename+"_b1"+".tif"), "#", "#", "-10000", "NONE", "NONE", "16_BIT_UNSIGNED", "NONE", "NONE") #save temporary singleband to disk list1.append(b1_rast)‍‍‍‍‍‍‍‍‍‍‍‍

This code successfully extracted each band from the multiraster into singleband tiffs. I haven't figured out the numpy array part yet but could probably save disk space by feeding the temporary rasters to it instead of using copyraster_management.

[DanPatterson_Retired]

skip the copyraster and use rastertonumpyarray as in the link I posted

[JacobNederend]

Still stuck with automating this. Arcpy exceeds available memory almost immedia which I think is due to the 32-bit Python. However, the rastertonumpyarray works at least. Still not sure how to substitute variables and loop over files in the directory. Code as it is now:

# Raster to Numpy Array tool for applying calibration equations

import arcpy
import numpy
import os

# Set raster directory
arcpy.env.workspace = "E:/ArcGIS_Work/ModelBuilder_Testing/ModelBuilder_NIR_clip"
rasters = arcpy.ListRasters("*","TIF")
rasterpath = "E:/ArcGIS_Work/ModelBuilder_Testing/ModelBuilder_NIR_clip/"
# Get input Raster properties
for raster_name in rasters:
    print(rasterpath+raster_name)
    raster = arcpy.Raster(rasterpath+raster_name)
    print (type(raster))
    
  
# Convert Raster to numpy array
    arrays = arcpy.RasterToNumPyArray(raster, nodata_to_value=-10000)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

Jacob... you seem stuck on this automating stuff without confirming that you can do a single case using numpy arrays (to them, and from them).  You might want to complete that portion first otherwise revert back to the manual method, since you would have been done by now... or is this going to be an ongoing venture?

[JacobNederend]

You are definitely right, I need to get it working on a single-raster basis first. I've had success there using GDAL and Rasterio. This will be an on-going task, with twice as many rasters to process in the immediate short-term and many more over the coming season. 

Thanks for the continued feedback!

[curtvprice]

could probably save disk space by feeding the temporary rasters to it instead of using copyraster_management.

This isn't technically correct. Although the documentation seems to suggest that temp rasters are not written to disk, they usually are as part of map algebra processing. I demonstrated this in this recent post:  https://community.esri.com/people/curtvprice/blog/2017/03/03/temporary-rasters-in-arcpy?sr=search&se...‌