[GIS] IDW interpolation of values from raster in ArcGIS

Question

I've been tasked with obtaining z-values for a large point set from a 2 m raster, but rather than extracting values directly from the raster, or with bilinear interpolation, I've been asked to use inverse distance weighted interpolation of cell centers in a 5 m radius.

This is beyond the capabilities of the Extract Values to Points tool or Zonal Statistics. As a workaround, I can convert the raster to points, resample to a finer resolution raster using IDW, and then extract values, but this isn't quite true IDW since it's limited to the cell size of the new raster (and by my patience for processing time).

Setting aside the fact that there are better approaches to smoothing elevations from a raster — this is the approach I've been asked to take — I could also script it in python fairly easily, but thought I might be missing something easy?

I have access to 3D Analyst, Spatial Analyst, etc in ArcGIS Desktop 10.3.

Answer 1

I've done some digging, and as far as I can tell there isn't an easy way to accomplish this. I built out a working process in Model Builder and it's quite cumbersome, but it works. I'm sure there are ways to improve this (scripting it would be the best way to go), but for anyone with a similar question, here are the steps in Model Builder:

There are four parameters:

1. Input points is the point feature class you want to interpolate via IDW
2. DEM is the raster to interpolate from
3. IDW Power is the power factor to apply, set to default of 2
4. Interpolation Radius is self-explanatory

Walking through the model:

1. Copy Features can be omitted if you don't mind adding fields to the input directly
2. Add XY Coordinates ensures there is a POINT_X and POINT_Y field for later use
3. Buffer and Clip are used to reduce the processing load in the subsequent Raster to Point -- for small rasters they can be omitted and the full raster converted to points. Use Input should be set to 'true' and linked to "Use input features for clipping geometry"
4. Spatial Join should take DEM points as the Target Feature and Input points as the Join Feature. Set Join Operation to ONE_TO_MANY and Match Option to WITHIN_A_DISTANCE
5. Add Field adds 'weight' field
6. Add Field (2) adds 'wz' field
7. Calculate Field calculates weight as
   • Expression: getw(!shape.getpart(0)!, !point_x!, !point_y!)
   • Expression Type: PYTHON 9.3 (use for all Calculate Fields)
   • Code block: def getw(point, grid_x, grid_y): #return idw weighting, defined as 1/d^p maxw = 1000 #increase for very fine grid raster p = %IDW Power% d = math.sqrt(math.pow(grid_x - point.X,2) + math.pow(grid_y - point.Y,2)) if d > 0: w = min(1/math.pow(d, p), maxw) else: w = maxw return w
8. Calculate Field (2) calculates wz as
   • Expression: !weight! * !grid_code!
9. Summary Statistics should have weight and wz both set to SUM in Statistics Field(s), with Case Field set to JOIN_FID
10. Join Field should take the following: Input Table -- Input points or Input points XY, Input Join Field -- FID, Join Table -- DEM points sum table, Output Join Field -- JOIN_FID, with Join Fields SUM_weight and SUM_wz checked
11. Add Field (3) adds z_idw field
12. Calculate Field (3) calculates z_idw as
    • Expression: !SUM_wz!/ !SUM_weight!

Because Model Builder is a bit unstable when it comes to field mapping, it helps to set this up with default point and raster fields the first time so it properly propagates field names for all the operations. Again, if you really want to improve stability and speed scripting would be preferable.