Looping numpy array to raster puts the wrong bands into new rasters

Question

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

Answer 1

Solved finally! Here's the code:

print "Beginning imports..."
import time
start_time = time.time()
import arcpy
import numpy as np
import os
import glob
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 =[]

# 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)    
print raster_names
for name in raster_names:
    rasters = arcpy.Raster(rasterpath + os.sep + name)
    # rasters = arcpy.Raster(rasterpath + 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 raster in 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)
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)
    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])
    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)