use a for loop to interpolate points and then extract interpolation with a boundary

1469
4
03-17-2014 11:54 AM
BenSciance
New Contributor II
I'm looking to get some input into interpolating multiple point feature classes using a for loop.

Whenever I specify to run only 1 feature class in the for loop, the script completes.  BUT, when I specify all feature classes, it loops through each step with the first feature class fine and when moving to the next feature class, I get an error saying "workspace or data source is read only" and "The operation was attempted on an empty geometry" or "no spatial reference exists" at the interpolation step.

I noticed that for each interpolation and extraction that works on a single feature class, a new raster grid is saved in my env.workspace path using default naming parameters, as well as in the savepaths that I have defined myself.  Is it possible that the reason why I get the "workspace is read only" error is due to python reading the default file that is saved to my workspace automatically? 

I've also tried using the ".save" function to store the new raster grids, which you will see are commented out because they did not work.  As an alternative I used the Copy Raster tool to save the output raster grids to new geodatabases.


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

#Set up environmental workspace
env.workspace = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\SGaugeTest.gdb"

#Set up save path
#clipsavepath = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauge_Clips.gdb"

krigsavepath = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\Krig_Test.gdb"

Extractsavepath = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\Krig_Mask_Test.gdb"

# Define constant variables
BasinBound = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauges.gdb\S_BasinBoundary"


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


fcs=arcpy.ListFeatureClasses("a*")
print "Listing feature classes"

for fc in fcs:
    print (fc)
    try:
        #outclip = clipsavepath + "/" + fc + "_" + "clip"
        outkrigsave = krigsavepath + "/" + fc + "_" + "krig"
        outmasksave = Extractsavepath + "/" + fc + "_" + "krigmask"
        zfield = "Wlevel"

        #Clip gauge data with basin boundary
        #print 'clipping gauge data'
        #arcpy.Clip_analysis(fc,BasinBound,outclip)

        # Run kriging Interpolate
        print 'interpolating gauge data'
        outkrig = Kriging(fc,zfield, KrigingModelOrdinary ("Spherical",0.0008333333333),0.0008333333333)
        #outkrig.save=(outkrigsave)

        # Copy krig to gdb
        print 'copying interpolation grid to new gdb'
        arcpy.CopyRaster_management(outkrig,outkrigsave)


        # Extract krig interpolation by the corresponding basin mask
        print 'extracting interpolation with basin mask'
        outmask = ExtractByMask(outkrig,BasinBound)
        #outmask.save = (outmasksave)

        # Copy krig to gdb
        print 'copying interpolation grid to new gdb'
        arcpy.CopyRaster_management(outmask,outmasksave)


    except:
    # If an error occurred print the message to the screen
        print arcpy.GetMessages()

print 'Script is complete!!'




Here is a copy of the error messages
*** Remote Interpreter Reinitialized  ***
>>> 
Listing feature classes
a734184
interpolating gauge data
copying interpolation grid to new gdb
extracting interpolation with basin mask
copying interpolation grid to new gdb
a734072
interpolating clipped gauge data
Executing: Kriging C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauges.gdb\a734072 Wlevel C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauges.gdb\Kriging_a7345 "Spherical 0.000833" 0.0008333333333 "VARIABLE 12" #
Start Time: Mon Mar 17 15:44:56 2014
ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only.
No spatial reference exists.
ERROR 010029: Unable to create the raster __1000000.
Failed to execute (Kriging).
Failed at Mon Mar 17 15:44:56 2014 (Elapsed Time: 0.00 seconds)


Thanks for taking the time to read through this.
Tags (2)
0 Kudos
4 Replies
curtvprice
MVP Esteemed Contributor
I think the save should work - if you use the right syntax. .save() is a method, the single argument is an output path for the raster to be saved to.

outkrig.save=(outkrigsave)


should be

outkrig.save(outkrigsave)
0 Kudos
NeilAyres
MVP Alum
You might want to include a
env.overwriteOutput = True

as well at the beginning of your script.
Good luck,
Neil
0 Kudos
BenSciance
New Contributor II
Thanks Neil and Curtis.  I made the edits you both suggested and the loop completed with no errors.
0 Kudos
BenSciance
New Contributor II
I've returned back to this script to use it again and unfortunately the same error message as before keeps occurring.

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).
Failed at Thu Apr 03 13:23:32 2014 (Elapsed Time: 1.00 seconds)
script complete



Here is the code with the suggestions you have made:
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\SW_Gauges\SW_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\SW_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("a732847*")
bounds="C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\SW_Gauges\SW_Gauges.gdb\SW_BasinBoundary"

print 'beginning fc for loop'
for fc in fcs:
    print fc
    try:
        zfield = "Wlevel"
        savekrig = krigpath +"/" + "SW_" + fc
        savemask = maskpath +"/" + "SW_" + 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'


Any other possible ideas as to why I keep getting these errors?

I'm stumped as to why the computer attempts to use the grids that are saved by default ("Extract_SW_a2" in the error code block) for the functions in the loop when I specifically call for the grid that was just saved (savekrig).

All grids that I call for in the loop have defined coordinate systems and are not opened in arcmap, arccatalog, or any other program....

I think the default grids that get created and used in the functions are what's causing these errors.  But Why?
0 Kudos