Latest Contributions by JacobNederend

Thanks! This worked so I am marking it as the answer. It did, however, introduce a new problem. The ExG map gets promoted to 64-bit double precision, which unfortunately doesn't work with the Binary Threshold function I need to apply from the Image Analysis Window. I posted about that previously, but never did manage to implement it myself: https://community.esri.com/thread/197481-why-does-arcpy-produce-different-results-than-raster-calculatormap-algebra   ... View more

Dan, Thanks so much for this. I have not yet had time to run it (my Python skills have not yet extended to defining functions) but I trust that it is correct so I am marking it as the answer.  ... View more

I made a mistake in the formula and forgot to carry over the normalization (e.g. r/r+g+b) from my original spatial analyst script. Altering your example with the division yields a divide by zero error: numpyarrays = arcpy.RasterToNumPyArray(infile, nodata_to_value=infile.noDataValue).astype(float) # Apply the CIVE formula of 0.441r - 0.881g +0.385b +18.78745 r = numpyarrays[0,:,:] g = numpyarrays[1,:,:] b = numpyarrays[5,:,:] args = [r, g, b] print("Band 1\n{}\nBand 2\n{}\nBand 6\n{}".format(*args)) CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745 print("Result...\n{}".format(CIVE)) #---------------------------------------------------------------------------------------------------------------------------------- C:/Users/Default.Default-PC/.PyCharmEdu3.5/config/scratches/scratch_5.py:48: RuntimeWarning: invalid value encountered in divide   CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745 C:/Users/Default.Default-PC/.PyCharmEdu3.5/config/scratches/scratch_5.py:48: RuntimeWarning: divide by zero encountered in divide   CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745 ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

This didn't work, unfortunately. Numpy would be better but am not sure how I would go about it. Here's what I am thinking: # assume that I already set up a list called composites, which are 6-band GeoTiffs for which I need to use bands 1,2, and 6 # Also assume that I have lists for all the relevant geotransform properties in order to convert output back to raster #N ot sure what to set NoData to in RasterToNumpyArray nparrays = [] for i, raster in enumerate(composites):     nparrays.append(arcpy.RasterToNumPyArray(raster, nodata_to_value=raster.noDataValue).astype(float)) # Apply the CIVE formula of 0.441r - 0.881g +0.385b +18.78745 CIVE_arrays = [] for date, numpyarrays in enumerate(nparrays):     CIVE = np.zeros_like (numpyarrays, np.float32)     # **not sure how to actually run the formula from here** because I don't know how to select dimensions 0,1, and 5     for i in range (0,5) #         CIVE[i,:,:] = 0.441*[i=0]-0.881[i=2]+0.385[i=6]+18.78745   # obviously the i= are invalid syntax ... View more

Also, the NoData value is -1.797693e+308. I also found this discussion from 2015 which you were also involved in. This is almost certainly my problem, but the OP never updated with a solution.  ... View more

Thanks for the advice, Dan. I think you are right about the bit-depth because the NoData values and symbologies are identical. The Excess Green Index I produced is 32-bit floating point, whereas the CIVE and all other indices are 64-bit double-precision. Could it be possible that the decimal coefficients in the CIVE index are messing up the bit-depth? If so, how could I fix the outputs to be 32-bit floating point? Some formulas are below for comparison. CIVE = 0.441r - 0.881g +0.385b +18.78745 ExG = 2g-r-b ExR = 1.4r-b ExGR = Exg-ExR ... View more

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) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

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. ... View more

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. ... View more

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! ... View more

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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

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. ... View more

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