2

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. However, I am not sure why the script is not working. Any help would be greatly appreciated.

 #-------------------------------------------------------------------------------
    # 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()
  • What error are you receiving? – Sarah Sep 20 '13 at 17:27
  • The error is "IDW failed". – Ken Dekalb Sep 20 '13 at 17:34
  • I'm not quite sure - but the command for IDW specifies that the third input should be the cell size? Idw (in_point_features, z_field, {cell_size}, {power}, {search_radius}, {in_barrier_polyline_features}) – Sarah Sep 20 '13 at 17:43
  • I deleted the r_level and left the code as gp.IDW_sa(stations, rain_amount, "1000", "2" "VARIABLE 12", "") and I still get the error "IDW failed". – Ken Dekalb Sep 20 '13 at 17:58
  • I'm guessing you got your code from here. If you read further down the changes needed are discussed? – Sarah Sep 20 '13 at 19:16
2

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')
| improve this answer | |
  • can you provide an example for importing a shapefile and executing this script? – Ken Dekalb Sep 21 '13 at 14:29
  • @KentonDekalb I'll work on that later on today – dassouki Sep 23 '13 at 16:07
0

Your 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 then what your error handler is reporting.

I don't have the cred to add this as a comment but 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)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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