topic Re: use a for loop to interpolate points and then extract interpolation with a bounda in Python Questions

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

BenSciance — Mon, 17 Mar 2014 18:54:33 GMT

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.

Re: use a for loop to interpolate points and then extract interpolation with a bounda

curtvprice — Tue, 18 Mar 2014 03:07:45 GMT

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)

Re: use a for loop to interpolate points and then extract interpolation with a bounda

NeilAyres — Fri, 10 Dec 2021 20:30:06 GMT

You might want to include a

env.overwriteOutput = True

as well at the beginning of your script.
Good luck,
Neil

Re: use a for loop to interpolate points and then extract interpolation with a bounda

BenSciance — Tue, 18 Mar 2014 12:38:02 GMT

Thanks Neil and Curtis. I made the edits you both suggested and the loop completed with no errors.

Re: use a for loop to interpolate points and then extract interpolation with a bounda

BenSciance — Fri, 10 Dec 2021 20:30:09 GMT

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?