I have a number of GeoTiffs for which I need to apply an equation (Y=Aexp(BX)) to each band. There are 5 rasters, so 15 bands total. I would essentially like to split the rasters apart, do the math then recomposite them the way they were. Below is my code, but the rasters it produces are all in correct. The Red band of image one is in the red band of image five instead, for example. I'm still new to Python and learning as I go so the looping is very rough and I would like to eventually set it up so that I don't need to individually load each raster.
The purpose of this script is to avoid the need to apply calibration equations manually using map algebra. As you might guess, doing that 15 times takes a long time and since I am using different A & B coefficients all the time the workload is too much to carry on like that.
*Note: I edited out the original script for all 5 rasters to save space. *Edit: solved the first script by moving the 'calibrated_data.append()' statement out one indentation. Still looking for a solution to the second script using lists instead of loading files individually as in script #2.
<!-- language: lang-py -->
import arcpy
import numpy as np
import os
import time
# 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
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
# List the raster properties
lowerLeft = [lowerLeft_1, lowerLeft_2, lowerLeft_3, lowerLeft_4, lowerLeft_5]
cellSize = [cellSize_1, cellSize_2, cellSize_3, cellSize_4, cellSize_5]
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)
print "Done..."
# all the images (in numpy arrays in a list)
dates = [d1]
# all coefficients and intercepts in lists for easy looping
print "Configuring coefficients..."
A = np.array([[0.0248262,0.0260896,0.0233547],dtype=np.float32)
B = np.array([[0.0000609, 0.0000637, 0.0000648],dtype=np.float32)
print "Done..."
print "Starting calibration calculations..."
calibrated_data = []
for date,numpyarray in enumerate(dates): #with enumerate now dates will be 1,2,3,4,5
arrayshape = numpyarray.shape
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
Y_rescaled[i,:,:] = (Y[i,:,:]-Y[i,:,:].min())/(Y[i,:,:].max()-Y[i,:,:].min())
calibrated_data.append(Y_rescaled)#_rescaled) # this has two members in it, the first is however many dates I have.
print "Done..."
print calibrated_data, type(calibrated_data)
print "Saving reflectancemaps as 32-Bit GeoTiff"
nir_calibrated1 = arcpy.NumPyArrayToRaster(calibrated_data[0], lowerLeft_1, cellSize_1)
arcpy.DefineProjection_management(nir_calibrated1,spatialReference_1)
nir_calibrated1.save(outputFolder+"//cal_nir_d1.tif")
print "Done..."
print "time elapsed: {:.2f}s".format(time.time() - start_time)
Edit: Here is a new attempt at shorter scripting with loops. The only puts out one singleband file so there's still something wrong with how I assign transformed arrays to image channels
import arcpy
import numpy as np
import os
import time
import glob
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 coefficients for Y=A*exp(B*raster_channel)
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 "Calibrating numpy arrays..."
calibrated_data = []
# List rasters in rasterpath
print "Listing rasters and their properties..."
rasters = arcpy.ListRasters("*", "TIF")
for raster in rasters:
ras = arcpy.Raster(rasterpath+raster)
lowerLeft = arcpy.Point(ras.extent.XMin, ras.extent.YMin)
cellSize = ras.meanCellWidth
spatialReference = ras.spatialReference
LL = []
LL.append(lowerLeft)
CS = []
CS.append(cellSize)
SR = []
SR.append(spatialReference)
print LL,CS,SR
numpyarray = arcpy.RasterToNumPyArray(ras, nodata_to_value=-10000)
numpyarray.astype(float)
arrays = []
arrays.append(numpyarray)
for date,numpyarrays in enumerate(arrays):
arrayshape = numpyarray.shape
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
Y_rescaled[i,:,:] = (Y[i,:,:]-Y[i,:,:].min())/(Y[i,:,:].max()-Y[i,:,:].min())
calibrated_data.append(Y_rescaled[i,:,:])#_rescaled) # this has two members in it, the first is however many dates I have.
for arrays in calibrated_data:
nir_calibrated = arcpy.NumPyArrayToRaster(arrays, lowerLeft, cellSize)
nir_calibrated.save(outputFolder + "//" + raster[:-4] + ".tif")
print "Done..."