current community

your communities

more stack exchange communities

Stack Exchange
tour help
Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 0 down vote favorite

Here is my code:

#Generate Cumulative Thresholds
#5/15/14

import os
import glob
import arcpy
import numpy as np

#Enable overwriting output data
arcpy.env.overwriteOutput = True

#########################################################################################################
def cumulativecalculation(i):

    #Set geoprocessing variables
    inraster = i
    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)

    #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)

This works perfectly for about 10 rasters, then "breaks." I get the error:

File "E:/NSF Project/Salamander_Data/PythonScripts/Calculate_Cumulative_Raster_Iterative.py", line 73, in cumulativecalculation(i)

File "E:/NSF >Project/Salamander_Data/PythonScripts/Calculate_Cumulative_Raster_Iterative.py", line 52, in cumulativecalculation outraster = arcpy.NumPyArrayToRaster(c, ll, des.meanCellWidth, des.meanCellHeight)

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'm not sure why this particular file seems to be causing trouble. Even if I end up with missing data, I'd at least like to be able to skip over input like this in order to process the rest of the rasters (around 400).

Any ideas on how to do that, or what might be throwing a monkey wrench in arcpy.NumPyArrayToRaster ?

share|improve this question

1 Answer 1

up vote 2 down vote accepted

Put your code in a try catch block and use the continue statement to continue even if you run into problems. This will ensure that one raster not processing does not stop the entire workflow. Something like this:

        try:
             for part in row[0]:
             #place your code here
        except:
        #Clean up by deleting temporary datasets and print the information of rows that were not able to process
              arcpy.Delete_management(tempLyr)
              print 'Could not process {0}.'.format(row[0])
              print 'Error: {0}.'.format(sys.exc_info()[1])
              ct += 1
              continue

So in your case just put your for loop in a try catch continue statement like this:

try:
     for i in geotiffs:
          cumulativecalculation(i)
except:
     #Code here to report which raster did not process
     continue
share|improve this answer
    
Thanks! Very helpful. –  Vergentorix May 21 '14 at 15:14
    
By the way, it turns out this is some kind of other issue. The script will only crank through 10 rasters at a time before showing this same error. I wonder if it's some kind of memory issue rather than bad input? Anyway, the trick you showed me is very instructive. –  Vergentorix May 21 '14 at 18:44
    
For the record I did answer your question "Even if I end up with missing data, I'd at least like to be able to skip over input like this in order to process the rest of the rasters". Are you sure that you have the continue statement underneath your except: block? It should still continue forward based on the error message you are giving. If it were a memory issue it would say MemoryError and not be giving the exception you have. If I have time later today I will take a deeper look at your code. –  JamesSLC May 21 '14 at 23:51
    
I copied your code exactly, aside from making sure indentations were correct. I started getting an error "SyntaxError: 'continue' not properly in loop." And couldn't figure out how to get around that. Everything I read afterward indicated that your solution should work, so I'm not sure what that is about. Re: The script breaking after 10 iterations, I tried breaking up the input into subdirectories in different configurations. Even rasters that seemed to have caused the problem before were converted without issue as long as they were not past tenth in line. –  Vergentorix May 22 '14 at 0:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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