Extract By Mask Loop Crashing

[RachelAlbritton]

My script keeps crashing, and I'm not sure why. No error messages are being printed, the program just quits.

import arcpy, os, shutil, datetime from arcpy import env from arcpy.sa import*  def outName(input,post="_Output"):     """Returns output name."""     outName=os.path.basename(input).split(".")[0]+post     return outName  #set workspaces arcpy.env.workspace = "C:\\TNC_GIS_Projects\\Surf_City" #arcpy.GetParameterAsText(0)  outputWorkspace = "C:\\TNC_GIS_Projects\\Surf_City\\Output_Folder" #arcpy.GetParameterAsText(1) arcpy.env.overwriteOutput = True  #Check out the ArcGIS 3D Analyst extension license arcpy.CheckOutExtension("3D") arcpy.CheckOutExtension("Spatial")  #Variables ObsPts = "C:\\TNC_GIS_Projects\\Surf_City\\Projected_Data\\Test_Points.shp" #arcpy.GetParameterAsText(2) footprint =  "C:\\TNC_GIS_Projects\\Surf_City\\Projected_Data\\Building_Footprints.shp" #arcpy.GetParameterAsText(3)  Elevation = "C:\\TNC_GIS_Projects\\Surf_City\\Projected_Data\\new_2010" #arcpy.GetParameterAsText(4) BareElevation = "C:\\TNC_GIS_Projects\\Surf_City\\Projected_Data\\baregrnd" #arcpy.GetParameterAsText(5) Ocean = "C:\\TNC_GIS_Projects\\Surf_City\\Projected_Data\\AtlanticOcean.shp" #arcpy.GetParameterAsText(6) FloorField = "FLOORS" #arcpy.GetParameterAsText(7)  Year = "2010" #arcpy.GetParameterAsText(8)  #Set analysis extent to elevation raster arcpy.env.extent = Elevation arcpy.env.cellSize = Elevation  #Create a temp folder to hold intermediate files. #Some of these folders will be deleted after the analysis has runs successfuly. IntermediateFiles = outputWorkspace+"/IntermediateFiles_"+Year Final_Floor_Viewsheds = outputWorkspace+"/Final_Viewsheds_"+Year SummaryTables = outputWorkspace+"/Summary_Tables_"+Year ElevAvgTables=outputWorkspace+"/ElevAvg_Tables_"+Year  if not os.path.exists(IntermediateFiles): os.makedirs(IntermediateFiles) if not os.path.exists(Final_Floor_Viewsheds) : os.makedirs(Final_Floor_Viewsheds) if not os.path.exists(SummaryTables) : os.makedirs(SummaryTables) if not os.path.exists(ElevAvgTables): os.makedirs(ElevAvgTables)  #Create Featue Layers arcpy.SetProgressorLabel("Creating feature layers") PointsFL = outName(ObsPts,"_Lyr") footprintFL = outName(footprint,"_Lyr") arcpy.MakeFeatureLayer_management(ObsPts, PointsFL) arcpy.MakeFeatureLayer_management(footprint,footprintFL)  RangeCount = int(arcpy.GetCount_management(PointsFL).getOutput(0))      arcpy.CalculateField_management(PointsFL,"OFFSETA",10,"PYTHON_9.3") arcpy.AddMessage("Observation point heights have been set to "+str(10)+" ft.")      arcpy.AddMessage("Calculating viewshed for "+str(RangeCount)+" parcel, for floor number "+str(1))      for points in range (0,RangeCount):      try:                      FID = "FID=%s" % (points)          ##Get bare earth elevation of parcel         #extract parcel raster cells         arcpy.SelectLayerByAttribute_management(PointsFL,"NEW_SELECTION",FID)         print arcpy.GetMessages(), "\n\n"         arcpy.SelectLayerByLocation_management(footprintFL,"INTERSECT",PointsFL)         print arcpy.GetMessages(), "\n\n"         outExtractByMask = ExtractByMask(BareElevation,footprintFL)         outExtractByMask.save(IntermediateFiles+"\\fp_1")#+str(1)+"_"+str(FID)[4:])         print arcpy.GetMessages(), "\n\n"              except:          print arcpy.GetMessages()  print "done" 

[StacyRendall1]

Try making the bottom part just:

for points in range (0,RangeCount):
            
        FID = "FID=%s" % (points)

        ##Get bare earth elevation of parcel
        #extract parcel raster cells
        arcpy.SelectLayerByAttribute_management(PointsFL,"NEW_SELECTION",FID)
        print arcpy.GetMessages(), "\n\n"
        arcpy.SelectLayerByLocation_management(footprintFL,"INTERSECT",PointsFL)
        print arcpy.GetMessages(), "\n\n"
        outExtractByMask = ExtractByMask(BareElevation,footprintFL)
        outExtractByMask.save(IntermediateFiles+"\\fp_1")#+str(1)+"_"+str(FID)[4:])
        print arcpy.GetMessages(), "\n\n"


print "done"

I have had trouble wrapping stuff in try...except statements before. In reality you want to know if something throws an exception...

[ChrisSnyder · Accepted solution]

I could be wrong, but I think the issue is this line:

outExtractByMask.save(IntermediateFiles+"\\fp_1")

As I recall, the .save() doesn't allow you to "move" the output raster to a different workspace than what arcpy.env.workspace is set to... Rather it simply "renames and makes permanent" the scratch raster. So your code would have to be just

outExtractByMask.save("fp_1")

...and the "fp_1" grid would be saved in where ever arcpy.env.workspace is set to at the time.

[RachelAlbritton]

I could be wrong, but I think the issue is this line:

outExtractByMask.save(IntermediateFiles+"\\fp_1")

As I recall, the .save() doesn't allow you to "move" the output raster to a different workspace than what arcpy.env.workspace is set to... Rather it simply "renames and makes permanent" the scratch raster. So your code would have to be just

outExtractByMask.save("fp_1")

...and the "fp_1" grid would be saved in where ever arcpy.env.workspace is set to at the time.

Chris - your right. What's odd is I haven't ran into this issue before and it exists other places in my script. I will be sure to change it.
I also noticed that even though I'm wanting to save the file with a specific name, its automatically naming the file, which means its being overwritten every time. Is there a way to prevent this?

[ChrisSnyder]

The way I get it, all the SA commands now generate a scratch raster in whatever workspace you have defined using arcpy.env.workspace. The scratch name is auto generated and my understanding is that there is no way to directly control it. Upon doing a .save() the scratch raster is renamed (and even reformatted if you specify a file extention) to your specified name and made to be permanant... Are you saying that upon running .save('my_name'), your raster will not retain the name you gave it?

[RachelAlbritton]

Got it - I still had to put a full directory to where it was being saved. Seems to be working now. Thanks!

[BenSciance]

I'm experiencing the same issue using the .save() after extracting by mask within a for loop.

The output grid is being saved in my env.workspace with default naming parameters.

Can you elaborate on when you "put a full directory to where it was being saved" which fixed it?

[BenSciance]

The same issues keep happening with my code.  Here is the code I'm working with:

import arcpy
from arcpy import env
from arcpy.sa import*

arcpy.env.workspace=r"C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\Sylhet_Gauges\Sylhet_Gauges.gdb"

env.overwriteOutput=True
arcpy.env.overwriteOutput=True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

krigpath= "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\Sylhet_Gauges\Krig4MinDiff.gdb"
maskpath= "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Minimize_DEM_Analysis\Mask4MinDiff.gdb"

print 'listing feature classes'
fcs=arcpy.ListFeatureClasses("a*")
bounds="Sylhet_BasinBoundary"

print 'beginning fc for loop'
for fc in fcs:
    print fc
    try:
        zfield = "Wlevel"
        savekrig = krigpath +"/" + "Syhlet_" + fc
        savemask = maskpath +"/" + "Sylhet_" + fc

        print ('kriging {0}'.format(fc))
        outKrig=arcpy.sa.Kriging(fc,zfield,KrigingModelOrdinary ("Spherical",0.0008333333333),0.0008333333333)
        outKrig.save(savekrig)


        print ('extracting {0} by {1}').format(fc,bounds)
        outmask=arcpy.sa.ExtractByMask(saveKrig,bounds)
        outMask.save(savemask)

    except:
        print arcpy.GetMessages()

print 'script complete'

This code works sometimes and completes a few iterations and then fails at the kriging step with this error message:

Executing: Kriging C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\Sylhet_Gauges\Sylhet_Gauges.gdb\a733547 Wlevel C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\Sylhet_Gauges\Sylhet_Gauges.gdb\Kriging_a7331 "Spherical 0.000833" 0.0008333333333 "VARIABLE 12" #
Start Time: Thu Apr 03 13:57:44 2014
ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only.
The operation was attempted on an empty geometry.
ERROR 010029: Unable to create the raster __1000000.
Failed to execute (Kriging).
Failed at Thu Apr 03 13:57:45 2014 (Elapsed Time: 1.00 seconds)

or the extraction step with the same error message:

extracting a732847 
Executing: ExtractByMask C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\krig4MinDiff.gdb/SW_a732847 by C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\SW_Gauges.gdb\SW_BasinBoundaryC:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\SW_Gauges.gdb\SW_BasinBoundary C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\SW_Gauges.gdb\Extract_SW_a2
Start Time: Thu Apr 03 13:23:31 2014
ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only.
The operation was attempted on an empty geometry.
ERROR 010302: Unable to create the output raster: C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\SW_Gauges.gdb\Extract_SW_a2
ERROR 010067: Error in executing grid expression.
Failed to execute (ExtractByMask).

In both cases, the grid is called that was saved by default, even though I call the grids that were just supposedly saved using the .save() function.  My guess is that these default grids are read only and do not have defined geometry? hence these error messages?

I'm stumped.