# Area precipitation using IDW?

I'm having the following problem with IDW automation. I have a catchment and several gauging stations measuring daily precipitation amounts. I would like to interpolate the precipitation amounts from these stations (points) to the whole catchment. In the next step, I would like to calculate average precipitation value for the catchment. I can do it for one day (1 value for each station). The problem is that I have 30 years of daily precipitation measurements in each station, and so I need to do this process thousands of times. I have created a simple model in model builder which could do this process once, but I am not able to force it to do it repeatedly.

How do I automate this process?

Here is the modelbuilder scheme: • Can you attach visual image of the model. it seems you need to use ittiration – Dror Har Gil Sep 9 '13 at 9:25
• Not an answer, but if your ultimate aim is to estimate average catchment precipitation, I'd be tempted to just average the different time series from your monitoring sites and cut out the IDW altogether. In my experience, simple IDW is a poor estimator of rainfall. Of course, this doesn't help if you need gridded data for other purposes. Good luck! – JamesS Sep 9 '13 at 13:54

There exists an extremely efficient solution to this problem that avoids computing ten thousand grids of interpolated data.

Using IDW to estimate precipitation tends to give poor results, so I will provide a more general answer that applies to better procedures such as Kriging as well as IDW.

What is common to these interpolation procedures is that they are linear preditors. This means that the value Z'(x,y) they predict at a given location (x,y) is a linear combination of the values Z(1), Z(2), ..., Z(n) at the monitoring stations and the weights in that linear combination are determined solely from the n locations and (x,y) (and do not depend on the observations). If we write these weights as w(i), i=1, ..., n, then explicitly

``````Z'(x,y) = w(1)*Z(1) + w(2)*Z(2) + ... + w(n)*Z(n)
``````

and each w(i) is a function of the locations (x,y), (x(1),y(1)), ..., (x(n),y(n)).

As an example, the IDW-p weights are

``````w(i) = d(i)^(-p) / [d(1)^(-p) + ... + d(n)^(-p)]
``````

where d(i) are the distances from (x,y) to (x(i),y(i)) (and typically p = 2). (It turns out you won't even need to know this formula, which is good because analogous formulas for better linear interpolators, such as various Kriging interpolators, are much more complicated!)

Consider, now, the problem of finding the average within a catchment (or any other predefined fixed region). Let's call the region R. This average will be found by predicting Z(x,y) separately for each grid cell within the region and averaging those predictions. Let the cell centers have coordinates (u(1),v(1)), ..., (u(m),v(m)) where m is the number of cells in the region. The average Z'(R) is defined to be

``````Z'(R) = [Z'(u(1),v(1)) + Z'(u(2),v(2)) + ... + Z'(u(m),v(m))] / m
``````

This, too, is a linear combination of the observations. Write w(i,j) for the weight applied to Z(i) when computing Z'(u(j),v(j)). Then the weight applied to Z(i) in the average is seen algebraically to equal

``````w'(i) = [w(i,1) + w(i,2) + ... + w(i,m)] / m
``````

That is, it's just the average weight within the region.

Now suppose the stations obtain repeated observations synchronously over time. Let's introduce time explicitly into the notation by writing Z(i,t) for the observation at station i at time t. Because the stations do not move and the same fixed grid is used to represent the regions, none of the weights change over time. Therefore the catchment average at each time t is still obtained with the same linear combination w'(1)*Z(1,t) + ... + w'(n)*Z(n,t). If we would like to average these over a period spanning observation times t(1), t(2), ..., t(s), then the average algebraically works out to

``````w'(1)*Z'(1) + ... + w'(n)*Z'(n)
``````

where Z'(i) is the time average at station i,

``````Z'(i) = [Z(i,t(1)) + Z(i,t(2)) + ... + Z(i,t(s))] / s.
``````

Therefore

The average in each catchment is obtained as a suitable weighted linear combination of the averages at each station

and

The weights depend only on the locations of the stations and the set of grid cells used to represent the catchment.

All we have to do is find these n weights and the rest is easy! This analysis has collapsed the calculation of 30*365 grids (of perhaps millions of points each) into the calculation of n weights, which is a huge reduction in effort (because it's almost always the case that the number of grid cells used to represent a region greatly exceeds the number of monitoring stations).

As we have seen, the final weights ultimately derive from a grid of weights, one per grid cell, computed for each monitoring station. These are the weights associated with station i at each cell center (u(j),v(j)). Here's the punch line: you can find all these weights with just n simple calculations. Artificially set the value at station i to 1 and set all the other values to 0. Run your interpolation routine (IDW or Kriging or whatever) on these synthetic data. As we have seen, the interpolated value at (u(j),v(j)) is

``````w(1,j)*Z(1) + w(2,j)*Z(2) + ... + w(n,j)*Z(n)
= w(1,j)*0 + ... + w(i-1,j)*0 + w(i,j)*1 + w(i+1,j)*0 + ... + w(n,j)*0
= w(i,j).
``````

You just read the weights right off the grid! This leads the the following simple recipe:

1. For i = 1 to n, do the following:

• Set Z(i) = 1 and Z(k) = 0 for k <> i.
• Run the interpolation routine, creating a grid of the weights w(i,j) (where j indexes the grid cells).
2. Average the w(i,j) within each catchment. (This is a zonal mean calculation, one for each station i.) This produces a collection of weights w(i,R) where i ranges over the stations and R ranges over the catchments.

3. Average the observations at each station, producing Z'(i) for each i.

4. The time average in catchment R is the linear combination w(1,R)*Z'(1) + ... + w(n,R)*Z'(n).

Total effort: n interpolations, n zonal means, n time-averages of the observations, and one linear combination for each region.

Things get complicated when some of the monitoring stations are missing data. The same principles and analysis apply, but different sets of weights have to be computed for different groupings of active stations. This, properly speaking, is actually a problem of statistical imputation and perhaps is better handled using more sophisticated imputation procedures, including those that capitalize on the spatial and temporal correlations that are surely present in the data (using methods of geostatistics and time series analysis).

• +1. As usual a great answer that will require some digesting. If I understand correctly, this method provides a clever “shortcut” to the average that would be obtained from interpolation, for a given set of station locations and a given grid. I can see that this is useful in itself, but in terms of estimating catchment average rainfall, is this better than simply averaging the station averages? I guess this method is allowing for the uneven spatial sampling in the raw data, but how well it does this presumably depends on how good your method is for estimating weights? – JamesS Sep 9 '13 at 17:18
• @JamesS Yes, it all depends on how good your interpolation method is. However, averaging the station averages is actually cruder than IDW-2, because it does not use the spatial information at all! – whuber Sep 9 '13 at 17:21