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.
The average in each catchment is obtained as a suitable weighted linear combination of the averages at each station
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
You just read the weights right off the grid! This leads the the following simple recipe:
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).
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.
Average the observations at each station, producing Z'(i) for each i.
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).