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

JacobNederend · ‎01-29-2017

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

JacobNederend · ‎03-14-2017

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.

JacobNederend · ‎03-15-2017

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.

curtvprice · ‎03-22-2017

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

JacobNederend · ‎03-23-2017

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

curtvprice · ‎03-23-2017

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)