I've put together a python script that does some calculations on rasters by converting them to NumPy arrays, doing the math, then converting them back to Raster and applying all the same spatial reference information, etc.
I just need to run it iteratively now on about 400 rasters.
I tried to use glob.glob to iterate through the rasters with a *tif wildcard, but couldn't get the output name to work right.
Now I'm trying to use os.path.basename with a conditional statement to filter through the files in the folder and only run my raster calculation on tif files.
When I try to run the script now, I get a "Process finished with exit code -1" message, and my print(outname) statement doesn't seem to register.
Apologies in advance if this is more of a python than an ArcPy question.
Also, thanks for your help.
import os import arcpy import numpy as np #Enable overwriting output data arcpy.env.overwriteOutput = True def cumulativecalculation(filename): #Set geoprocessing variables inraster = filename des = arcpy.Describe(inraster) sr = des.SpatialReference ext = des.Extent ll = arcpy.Point(ext.XMin, ext.YMin) #Convert GeoTIFF to numpy array a = arcpy.RasterToNumPyArray(inraster) #Flatten for calculations a.flatten() #Find unique values, and record their indices to a separate object a_unq, a_inv = np.unique(a, return_inverse=True) #Count occurences of array indices a_cnt = np.bincount(a_inv) #Cumulatively sum the unique values multiplied by the number of #occurences, arrange sums as initial array b = np.cumsum(a_unq * a_cnt)[a_inv] #Divide all values by 10 (reverses earlier multiplication done to #facilitate accurate translation of ASCII scientific notation #values < 1 to array) b /= 10 #Rescale values between 1 and 100 maxval = np.amax(b) b /= maxval b *= 100 #Restore flattened array to shape of initial array c = b.reshape(a.shape) #Convert the array back to raster format outraster = arcpy.NumPyArrayToRaster(c, ll, des.meanCellWidth, des.meanCellHeight) #Set output projection to match input arcpy.DefineProjection_management(outraster, sr) #Setting the OutName outname = filename + "_cumulative" + ".tif" #Save the raster as a TIFF outraster.save("E:/NSF Project/Salamander_Data/New_Cumulative_Output/" + outname) print (outname) #Loop through files in directory for filename in os.path.basename("E:/NSF Project/Salamander_Data/NoDataToZero/HadleyGCM/"): if filename[-4:] == ".tif": cumulativecalculation(filename)
Updated Code:
#Set the output filename
intername1 = os.path.splitext(os.path.basename(i))[0]
filename = intername1[5:]
finalname = filename + "_cumulative" + ".tif"
#Save the raster as a TIFF
outraster.save("E:\\NSF Project\\Salamander_Data\\New_Cumulative_Rasters\\" + str(finalname))
#Announce file creation
print("Created: " + str(finalname))
#Loop through files in directory
geotiffs = glob.glob("E:\\NSF Project\\Salamander_Data\\NoDataToZero\\HadleyGCM\\*tif")
for i in geotiffs:
cumulativecalculation(i)
Everything seems to be working beautifully with the code now, except it "breaks" after processing 10 rasters. I think this must be a problem with one of the input files, but I'm not sure.
Error: "File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy__init__.py", line 1860, in NumPyArrayToRaster return _NumPyArrayToRaster(*args, **kwargs) TypeError: Cannot create raster for numpy array."
I have about 400 to get through, so I'd at the very least like to keep errors like this from derailing the process, even if I end up with some missing data. Not sure how to accomplish that though.