1

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.

2
  • How are you removing the extension from filename on this line OutName = str(filename) + "_cumulative" + ".tif"
    – artwork21
    May 21, 2014 at 2:30
  • I've since realized that was a problem, so now I'm trying to use os.path.basename to do this. I'm updating the code now. May 21, 2014 at 2:40

1 Answer 1

2

In the for statement that creates your loop:

for filename in os.path.basename("E:/NSFProject/Salamander_Data/NoDataToZero/HadleyGCM/"):

there is a bit of a syntax error. Because you are ending your path in a forward slash nothing is being returned for your basename. It is returning an empty string '' because there is no basename for that, so filename[-4:] does not return anything either as it is an empty string. Additionally os.path.basename just returns a string so even if you don't end on a forward slash your variable filename will simply be the first letter in the string that is returned by your basename, so in your case (without the forward slash) would be H and then sequentially loop through a-d-l-e-y-G-C-M.

Here is an example of a simple script using glob.glob to find matches of files in a folder with a .dat extension. The code then uses os.path.dirname, os.path.basename, as well as os.path.splitext to create an output file path as the second argument to the dat_to_shp function.

matches = glob.glob(r'C:\Users\Desktop\Data\MethodsData\*.dat')

for i in matches:
     dat_to_shp(i, os.path.dirname(i) + '\\' + os.path.splitext(os.path.basename(i))[0] + '.shp')
     print 'Created file {0} shapefile'.format(os.path.splitext(os.path.basename(i))[0])

For further clarification if we take the following path C:\User\Test\raster.tif, then the following python functions will have these effects:

enter image description here

Hopefully there is enough info in this example for you to get things straight with your script.

1
  • Thanks! This was very helpful. I am now encountering a new problem, but maybe that should become a new question, since you answered the original one. May 21, 2014 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.