Interpolation iteration using ArcPy

Question

My current project requires interpolation of state level weather data.

Approximately 5 years of daily data needs to be interpolated. I still need to figure out the best interpolation method. However, for this exercise I chose IDW. I found a script in ArcGIS forum and tried to use it.

Why is the script not working?

 #-------------------------------------------------------------------------------
 # Name: Inverse Distance Method Interpolation
 # Purpose:
 #
 # Author: 
 #
 # Created: 19/09/2013
 # Copyright: 
 # Licence: 
 #-------------------------------------------------------------------------------
 # ---------------------------------------------------------------------------
 # rainfall_idw2012.py
 # Created on: Tue Apr 08 2008 10:47:39 AM
 # ---------------------------------------------------------------------------
 # Import system modules
 import arcpy
 import sys
 import os
 import arcgisscripting
 import string
 
 arcpy.env.overwriteOutput = True
 
 # Create the Geoprocessor object
 gp = arcgisscripting.create()
 
 # Check out any necessary licenses
 gp.checkoutextension("Spatial")
 
 # Load required toolboxes...
 gp.AddToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx")
 
 #Output_variance_of_prediction_raster = ""
 modelPath = r"D:\Geodatabase\Rainfall\r2012"
 
 
 # Setting workspace
 gp.workspace = modelPath
 
 
 # Input Variables
 stations = r"D:\Geodatabase\MyGeodatabase.gdb\mr2012"
 studyArea = r"D:\Geodatabase\MyGeodatabase.gdb\boundary"
 
 
 try:
 # Get list of dau=ily events for looping
 daily_events = gp.ListFields(r"D:\Geodatabase\MyGeodatabase.gdb\mr2012", "r2012*")
 
 daily_events.reset()
 
 # Get the first reading and start loop
 
 dailyevent = daily_events.Next()
 
 while dailyevent:
 
 rain_amount = dailyevent.Name
 
 # Build output data layer name from the rainfall:
 
 r_level = modelPath + _rain_amount
 
 extract = modelPath + _rain_amount
 dailyevent = daily_events.Next()
 
 gp.AddMessage ("Current Operation: Preparing Interpolating using IDW")
 
 # Process: IDW interpolation
 
 gp.IDW_sa(stations, rain_amount, r_level, "1000", "2" "VARIABLE 12", "")
 
 
 
 # Process: Extract by Mask...
 gp.ExtractByMask_sa(r_level, studyArea, extract)
 
 except:
 gp.AddMessage("IDW failed " + gp.GetMessages())
 print "IDW failed ", gp.GetMessages()

• Sarah: I'm guessing you got your code from here. If you read further down the changes needed are discussed?
• Kenny D: You are right, I got the code from ESRI Forum. I had already incorporated those changes! Still no luck...am I missing something?

Answer 1

I've successfully used these two functions for thousands of interpolations and all were unittested. IT is based on my knowledge as well as help from stackoverflow

import numpy as np
from scipy.ndimage import map_coordinates
def oneD_interpolate(x, x_list, y_list):
 """
 interpolate in one dimension
 """
 return np.interp(x, x_list, y_list)
def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1):
 """
 interpolate in two dimensions with "hard edges"
 """
 arr = np.atleast_2d(arr)
 ny, nx = arr.shape # Note the order of ny and xy
 x1 = np.atleast_1d(x1)
 y1 = np.atleast_1d(y1)
 # Change coordinates to match array
 # IF for example you have to interpolate across an array that
 # maximizes at 5000 or is minimum at 0
 if nx == 1:
 x1 = np.zeros_like(x1.shape)
 else:
 x_stride = (xmax - xmin) / float(nx - 1)
 x1 = (x1 - xmin) / x_stride
 if ny == 1:
 y1 = np.zeros_like(y1.shape)
 else:
 y_stride = (ymax - ymin) / float(ny - 1)
 y1 = (y1 - ymin) / y_stride
 # order=1 is required for the most part; however, read the proper
 # documentation if you need anything else
 return map_coordinates(arr, np.vstack((y1, x1)), order=1, mode='nearest')

Answer 2

You are getting the non-descriptive error "IDW failed" because of your error handler. Take the code out of the error handler, re-run, then report the traceback.It would be far more useful to everyone than what your error handler is reporting.

Your code needs to get cleaned up too. There is some redundancy there, specifically in importing both arcpy and arcgisscripting.

Try this and report back your error.

# Import system modules
import arcpy, sys, os, arcgisscripting, string
arcpy.env.overwriteOutput = True
# Create the Geoprocessor object
gp = arcgisscripting.create()
# Check out any necessary licenses
gp.checkoutextension("Spatial")
# Load required toolboxes...
gp.AddToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx")
#Output_variance_of_prediction_raster = ""
modelPath = r"D:\Geodatabase\Rainfall\r2012"
# Setting workspace
gp.workspace = modelPath
# Input Variables
stations = r"D:\Geodatabase\MyGeodatabase.gdb\mr2012"
studyArea = r"D:\Geodatabase\MyGeodatabase.gdb\boundary"
# Get list of dau=ily events for looping
daily_events = gp.ListFields(r"D:\Geodatabase\MyGeodatabase.gdb\mr2012", "r2012*")
daily_events.reset()
# Get the first reading and start loop
dailyevent = daily_events.Next()
while dailyevent:
 rain_amount = dailyevent.Name
 # Build output data layer name from the rainfall:
 r_level = modelPath + _rain_amount
 extract = modelPath + _rain_amount
 dailyevent = daily_events.Next()
 gp.AddMessage ("Current Operation: Preparing Interpolating using IDW")
 # Process: IDW interpolation
 gp.IDW_sa(stations, rain_amount, r_level, "1000", "2" "VARIABLE 12", "")
 # Process: Extract by Mask...
 gp.ExtractByMask_sa(r_level, studyArea, extract)