topic Re: Substitute different equation parameters (e.g. slope) when batch processing raster calculator? in ArcGIS Spatial Analyst Questions

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

JacobNederend — Mon, 30 Jan 2017 03:57:23 GMT

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:

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

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

DanPatterson_Retired — Sat, 11 Dec 2021 11:56:43 GMT

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.

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

JacobNederend — Mon, 30 Jan 2017 15:56:21 GMT

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

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

DanPatterson_Retired — Mon, 30 Jan 2017 19:08:22 GMT

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

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

JacobNederend — Sat, 11 Dec 2021 11:56:46 GMT

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.

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

DanPatterson_Retired — Wed, 08 Feb 2017 04:03:53 GMT

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

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

JacobNederend — Sat, 11 Dec 2021 11:56:48 GMT

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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

DanPatterson_Retired — Wed, 15 Feb 2017 15:58:16 GMT

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?

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

JacobNederend — Thu, 16 Feb 2017 21:55:45 GMT

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!

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

JacobNederend — Sat, 11 Dec 2021 11:56:51 GMT

import arcpy
import numpy as np
import os
print ("Done importing...")


# Set raster directory
print('Setting raster directory...')
arcpy.env.workspace = "E:/ArcGIS_Work/NIR_rasters)"

rasterpath = "E:/ArcGIS_Work/NIR_rasters)"
outputFolder = "E:/ArcGIS_Work/calibrated_rasters)
"print('Done...')

# Get input Raster properties 
print "Loading rasters and their extent + cell size..."
nir_d1 = arcpy.Raster("E:/ArcGIS_Work/NIR_rasters/NIR_d1.tif")
lowerLeft_1 = arcpy.Point(nir_d1.extent.XMin, nir_d1.extent.YMin)
cellSize_1 = nir_d1.meanCellWidth
spatialReference_1 = nir_d1.spatialReference
nir_d2 = arcpy.Raster("E:/ArcGIS_Work/NIR_rasters/NIR_d2.tif")
lowerLeft_2 = arcpy.Point(nir_d2.extent.XMin, nir_d2.extent.YMin)
cellSize_2 = nir_d2.meanCellWidth
spatialReference_2 = nir_d2.spatialReference
nir_d3 = arcpy.Raster("E:/ArcGIS_Work/NIR_rasters/NIR_d3.tif")
lowerLeft_3 = arcpy.Point(nir_d3.extent.XMin, nir_d3.extent.YMin)
cellSize_3 = nir_d3.meanCellWidth
spatialReference_3 = nir_d3.spatialReference
nir_d4 = arcpy.Raster("E:/ArcGIS_Work/NIR_rasters/NIR_d4.tif")
lowerLeft_4 = arcpy.Point(nir_d4.extent.XMin, nir_d4.extent.YMin)
cellSize_4 = nir_d4.meanCellWidth
spatialReference_4 = nir_d4.spatialReference
nir_d5 = arcpy.Raster("E:/ArcGIS_Work/NIR_rasters/NIR_d5.tif")
lowerLeft_5 = arcpy.Point(nir_d5.extent.XMin, nir_d5.extent.YMin)
cellSize_5 = nir_d5.meanCellWidth
spatialReference_5 = nir_d5.spatialReference

# List the raster properties
lowerLeft = [lowerLeft_1], [lowerLeft_2], [lowerLeft_3], [lowerLeft_4], [lowerLeft_5]]
print lowerLeft, type(lowerLeft)

cellSize = [[cellSize_1], [cellSize_2], [cellSize_3], [cellSize_4], [cellSize_5]]
print cellSize, type(cellSize)

spatialReference = [[spatialReference_1], [spatialReference_2], [spatialReference_3], [spatialReference_4], [spatialReference_5]]
print spatialReference, type(spatialReference)

print "Done..."

# Convert Raster to numpy array
print "Converting rasters to NumPy arrays..."

d1 = arcpy.RasterToNumPyArray(nir_d1, nodata_to_value=-10000)
d1.astype(float)
d2 = arcpy.RasterToNumPyArray(nir_d2, nodata_to_value=-10000)
d2.astype(float)
d3 = arcpy.RasterToNumPyArray(nir_d3, nodata_to_value=-10000)
d3.astype(float)
d4 = arcpy.RasterToNumPyArray(nir_d4, nodata_to_value=-10000)
d4.astype(float)
d5 = arcpy.RasterToNumPyArray(nir_d5, nodata_to_value=-10000)
d5.astype(float)
print "Done..."
# all the images (in numpy arrays) in a list
dates = [d1, d2, d3, d4, d5]

# all coefficients and intercepts in lists for easy looping
print "Configuring coefficients..."
A = np.array([[0.0248262,0.0260896,0.0233547],[0.0065193, 0.0068043,0.0018543],[0.0173419,0.0159686,0.0089509], [0.0166572,0.0161596,0.0109827],[0.0377285,0.0582130,0.0319386]],dtype=np.float32)
B = np.array([[0.0000609, 0.0000637, 0.0000648], [0.0001313,0.0001401,0.0001643],[0.0000772,0.0000845,0.0000900], [0.0000707,0.0000763,0.0000782],[0.0000482,0.0000431,0.0000508]],dtype=np.float32)
print "Done..."
calibrated_data = []
print "Starting calibration calculations..."
# This can be changed to however you'd like to store the values
for date,numpyarray in enumerate(dates): #with enumerate now dates will be 1,2,3,4,5
    arrayshape = numpyarray.shape
    print arrayshape
    Y = np.zeros_like(numpyarray, np.float32)
    Y_rescaled = np.zeros_like(numpyarray, np.float32)
    for i in range (0,3):
        Y[i,:,:] = A[date,i]*np.exp(B[date,i]*numpyarray[i,:,:]) # B[date,i] should be a scalar
        print Y
        # print A[date,i], B[date,i], numpyarray[i,:,:]
        Y_rescaled[i,:,:] = (Y[i,:,:]-Y[i,:,:].min())/(Y[i,:,:].max()-Y[i,:,:].min())
        print Y_rescaled
        calibrated_data.append(Y_rescaled) # this has two members in it, the first is however many dates I have.
print "Done..."
print "Saving rasters as 3-band 32-Bit GeoTiff..."
# for newarray in calibrated_data:
#     nir_calibrated = arcpy.NumPyArrayToRaster(newarray, lowerLeft, cellSize)
#     arcpy.DefineProjection_management(nir_calibrated,spatialReference)
#     nir_calibrated.save(outputFolder + "//" + "nir" + date + ".tif")


nir_calibrated1 = arcpy.NumPyArrayToRaster(calibrated_data[0], lowerLeft_1, cellSize_1)
arcpy.DefineProjection_management(nir_calibrated1,spatialReference_1)
nir_calibrated1.save(outputFolder+"//nir_cal_d1.tif")

nir_calibrated2 = arcpy.NumPyArrayToRaster(calibrated_data[1], lowerLeft_2, cellSize_2)
arcpy.DefineProjection_management(nir_calibrated2,spatialReference_2)
nir_calibrated2.save(outputFolder+"//nir_cal_d2.tif")

nir_calibrated3 = arcpy.NumPyArrayToRaster(calibrated_data[2], lowerLeft_3, cellSize_3)
arcpy.DefineProjection_management(nir_calibrated3,spatialReference_3)
nir_calibrated3.save(outputFolder+"//nir_cal_d3.tif")

nir_calibrated4 = arcpy.NumPyArrayToRaster(calibrated_data[3], lowerLeft_4, cellSize_4)
arcpy.DefineProjection_management(nir_calibrated1,spatialReference_4)
nir_calibrated4.save(outputFolder+"//nir_cal_d4.tif")

nir_calibrated5 = arcpy.NumPyArrayToRaster(calibrated_data[4], lowerLeft_5, cellSize_5)
arcpy.DefineProjection_management(nir_calibrated1,spatialReference_5)
nir_calibrated5.save(outputFolder+"//nir_cal_d5.tif")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Here is my final code. I'm sure there are improvements to be made, and I would really like to improve the looping so that I don't have to list each raster, thus the script would be more generic. But it works, and adding more rasters is a simple matter of copying and modifying a few lines.

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

JacobNederend — Wed, 15 Mar 2017 07:01:20 GMT

Upon further inspection, this code is not working correctly. The numpyarraytoraster functions are putting the wrong dates and channels into the GeoTiffs. E.g. Date1-Band1 is actually in Date3-Band1.

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

curtvprice — Thu, 23 Mar 2017 05:16:49 GMT

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&searchId=a8676672-f7e0-41eb-885c-f75654ceade5&searchIndex=0‌

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

curtvprice — Sat, 11 Dec 2021 11:56:54 GMT

Here's my approach, using arcpy map algebra. This code is not tested.

import os
import arcpy
from arcpy import env
from arcpy.sa import Exp, Int
arcpy.CheckOutExtension("spatial")

src_dir = "E:\\ArcGIS_Work\\raw_rasters"
dst_dir = "E:\\ArcGIS_Work\\cal_rasters"

# Create a list of input rasters
env.workspace = src_dir
rasters = arcpy.ListRasters("*","TIF")

# work in output workspace
env.workspace = dst_dir
env.scratchWorkspace = dst_dir

# load dictionary with calibration parameters for image "a"
# you could set up keys based on the raster name here
caldict = {"a": [2, 3.4,  3, 4.3,  4, 1.2])

# process each raster one at a time
for r in rasters:
    rpth = os.path.join(src_dir, r)
    # create layers to extract bands
    b1 = arcpy.MakeRasterLayer_management(rpth, "b1", index="1")
    b2 = arcpy.MakeRasterLayer_management(rpth, "b2", index="2")
    b3 = arcpy.MakeRasterLayer_management(rpth, "b3", index="3")
    # retrieve calibration coefficents from dict 
    b1a, b1b, b2a, b2b, b3a, b3b = caldict["a"]
    # calculate calibrated rasters for each band
    b1c = Int(b1a * Exp(b1b * b1))
    b2c = Int(b2a * Exp(b2b * b2))
    b3c = Int(b3a * Exp(b3b * b3))
    # save rasters to disk to pass to composite bands tool
    b1c.save()
    b2c.save()
    b3c.save()
    # save using same raster name as input
    arcpy.CompositeBands_management([b1c, b2c, b3c], r)
    # Clean up
    for tt in [b1c, b2c, b3c, b1, b2, b3]: 
        arcpy.Delete_management(tt)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

JacobNederend — Sat, 11 Dec 2021 11:56:56 GMT

Curtis, thanks so much for this! I have a further simplified version of my code above, but I like how short yours is and will try it out. I also have a problem with my own script where my rescaling step (normalize range to 0-1) produces poorer results (poorer correlation to ground truth data) than if I did everything manually in the GUI using raster calculator. How would I go about rescaling using map algebra following the workflow you present?

Here is my current workflow:

print "Beginning imports..."
import time
start_time = time.time()
import arcpy
import numpy as np
import os
import glob
print "time elapsed: {:.2f}s".format(time.time() - start_time)
print "Done importing..."
# Set Infolder/Outfolder
print 'Setting raster directory...'
arcpy.env.workspace = "E:/ArcGIS_Work/NIR_rasters"
rasterpath = "E:/ArcGIS_Work/NIR_rasters/"
outputFolder = "E:/ArcGIS_Work/calibrated_rasters"
print 'Done...'

# Set empty lists to fill while looping
nirs=[]
lowerLeft =[]
cellSize =[]
spatialReference =[]
nodata=[]
# List rasters in rasterpath. Note: Could use arcpy.ListRasters
# but this method is more generic to Python, thus easy to adapt
# to GDAL, OpenCV, Rasterio, etc.
print "Listing rasters and their properties..."
raster_names = []
for root, dirs, files in os.walk(rasterpath):
    for file in files:
        if file.endswith('.tif'):
            raster_names.append(file)
# raster_names = [(file for file in rasterpath if file.endswith('.tif')) for _,__,files in os.walk(rasterpath)]
print raster_names
for name in raster_names:
    rasters = arcpy.Raster(rasterpath + os.sep + name)
    nirs.append(rasters)
    lowerLeft.append(arcpy.Point(rasters.extent.XMin, rasters.extent.YMin))
    cellSize.append(rasters.meanCellWidth)
    spatialReference.append(rasters.spatialReference)
print rasters
print "Converting rasters to numpy arrays..."
nparrays = []
for i, raster in enumerate(nirs):
    nparrays.append(arcpy.RasterToNumPyArray(raster, nodata_to_value=-10000).astype(float))
print nparrays
print "time elapsed: {:.2f}s".format(time.time() - start_time)
print "Done..."
# Set coefficients for Y=A*exp(B*raster_channel)
# Lists are in format of [[date 1 channels]...[date 5 channels]] aka [[R1,G1,B1],...[R5,G5,B5]]
print "Listing coefficients..."
A = np.array([
        [0.0248262,0.0260896,0.0233547],
        [0.0065193,0.0068043,0.0018543],
        [0.0173419,0.0159686,0.0089509],
        [0.0166572,0.0161596,0.0109827],
        [0.0377285,0.0582130,0.0319386]
    ],dtype=np.float32)
B = np.array([
        [0.0000609,0.0000637,0.0000648],
        [0.0001313,0.0001401,0.0001643],
        [0.0000772,0.0000845,0.0000900],
        [0.0000707,0.0000763,0.0000782],
        [0.0000482,0.0000431,0.0000508]
    ],dtype=np.float32)
print "Done..."
# Apply the calibration equations
print "Starting array calibrations..."
calibrated_data = []
for date, numpyarrays in enumerate(nparrays):
    Y = np.zeros_like (numpyarrays, np.float32)
    Y_rescaled = np.zeros_like (numpyarrays, np.float32) # will store arrat re-scaled b/w 0-1 since reflectance cannot be negative.
    for i in range (0,3):
        Y[i,:,:] = A[date,i]*np.exp(B[date,i]*numpyarrays[i,:,:]) # B[date,i] should be a scalar
        Y_rescaled[i,:,:] = (Y[i,:,:]-Y[i,:,:].min())/(Y[i,:,:].max()-Y[i,:,:].min())
        print Y_rescaled
    calibrated_data.append(Y_rescaled) # this has two members in it, the first is however many dates I have.
print "time elapsed: {:.2f}s".format(time.time() - start_time)
print "Done..."
# Save rasters back to GeoTiff
print "Saving nparray to raster..."
for i, data in enumerate(calibrated_data):
    calibrated = arcpy.NumPyArrayToRaster(data, lowerLeft[i], cellSize[i],value_to_nodata=-9999)
    arcpy.DefineProjection_management(calibrated,spatialReference[i])
    calibrated.save(outputFolder + "//" + (raster_names[i])[:-4] + ".tif")
print "Done..."
print "time elapsed: {:.2f}s".format(time.time() - start_time)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

curtvprice — Fri, 24 Mar 2017 02:34:38 GMT

You could use map algebra to linearly scale your output from 0 to 1 (though 0 to 255 is often used).

ras1 = (ras - ras.minimum) / (ras.maximum - ras.minimum)

However you may get better performance using Rescale by Function—Help | ArcGIS Desktop

ras1 = RescaleByFunction(TfLinear(), from_scale=0.0, to_scale=1.0)