Using Spatial Analyst, we can expect each iteration of IDW, thresholding, and accumulation of results to take from a fraction of a second to a few seconds. Multiplying by 24 * 365 = 8760 hours in a year suggests this calculation will take hours at the least.
It can be much faster. Faster is not always better, but here--because there are significant issues concerning the interpolation method, since IDW is unlikely to be accurate for such data--it will be real advantage to be able to carry out the calculations quickly on a realistic scale in order to see the results and experiment with other options.
To illustrate this and to estimate how quickly it can really be done, I implemented the full-scale calculation in R using simulated data. (It is a simple matter to read your data into R and export the resulting grid in ESRI format using, say, the Raster package.) After the usual preliminaries of getting the data into the right format (values are stored in a 60 by 8760 array and point locations in a 60 by 2 array of (x,y) coordinates), handling missing values, and so on, the entire computation comes down to two lines of code:
thresh <- threshold * xy %*% one # Compute the denominators
a <- matrix(rowMeans(xy %*% z >= thresh), n.x, n.y)
xy has precomputed the IDW weights (unnormalized). one indicates where data are not missing and z is the matrix of all the data. The two matrix multiplications xy %*% one and xy %*% z compute the denominators and numerators of all the 8760 IDW calculations at once. The comparison in the second line creates the 8760 grids of ones and zeros and rowMeans averages them for the entire year. The result, a, is a single grid of dimensions n.x and n.y.
By reducing this large number of linear interpolations to a few matrix operations, the fastest possible speed and efficiency can be attained. Moreover, this approach is general: literally any linear interpolation method (such as Kriging) can be implemented in this fashion.
On my workstation (running this on a single core even though massive parallelization is possible) it takes just under 60 seconds to process one-quarter of the grid. By reducing the grid resolution little detail is lost and the entire year's data can be processed in seconds.
To illustrate, here are the results for 60 points and 8760 hours on a 125 by 310 grid (half resolution; 58 seconds total time). For simplicity it uses all available data for the interpolation. If smaller local neighborhoods are desired, they would be involved in one (very fast) precomputation step, where the weights xy are computed, but otherwise nothing else would change.
The windspeeds during the first hour are plotted as circles whose areas are proportional to the average speed during the year; the averages range from 1 to 5.9 m/sec.
There is a price to pay: this method precomputes some values and stores them in RAM. This quarter-size calculation required 6.75 GB RAM. A full-size calculation would therefore request 27 GB RAM. If you don't have that much, though, it is a simple matter to carve up the grid into smaller rectangles, do the calculation for each rectangle separately, and mosaic the pieces afterwards. The total calculation time would not appreciably increase.
Complete code follows.
set.seed(23)
n.points <- 60
n.hours <- 24 * 365
n.x <- 620/2 # Number of columns in grid
n.y <- 250/2 # Number of rows in grid
threshold <- 5 # Wind speed threshold
#
# The data: point locations.
#
points <- matrix(rnorm(n.points*2, 10^5, 10^4), ncol=2)
points <- points[order(points[,1]), ] # Sort by x-coordinate for testing
#
# Expand the extent of the points by about 10% to define the extent of the grid.
#
ranges <- apply(points, 2, range)
d.r <- diff(ranges)
x <- seq(ranges[1,1] - .05*d.r[1], ranges[2,1] + 0.05*d.r[1], length.out=n.x)
y <- seq(ranges[1,2] - .05*d.r[2], ranges[2,2] + 0.05*d.r[2], length.out=n.y)
#
# Obtain the inverse squared distances to the points.
#
distance.to <- function(u, x, y, eps=10^(-8)) {
e <- eps * (max(x) - min(x) + max(y) - min(y))^2
return(1 / (e + (x-u[1])^2 + (y-u[2])^2))
}
system.time(
xy <- apply(points, 1,
function(u) outer(x, y, function(x0, y0) distance.to(u, x0, y0)))
)
#
# The data: an array of values, one per column.
#
p.x <- points[,1]
z.scale <- rgamma(n.points, 5, scale=((p.x - min(x))/10^4 + 2)/50)
z <- matrix(rgamma(n.points * n.hours, 10, scale=z.scale), nrow=n.points)
z[z > 10] <- NA # Introduce some missing values
#
# Perform the interpolation, thresholding, and averaging.
#
system.time({
one <- !is.na(z) # Indicates the valid values
z[is.na(z)] <- 0 # Replace missing with 0
thresh <- threshold * xy %*% one # Compute the denominators
a <- matrix(rowMeans(xy %*% z >= thresh), n.x, n.y)
})
#
# Plot the results.
#
z.sum <- rowSums(z)
z.max <- max(z.sum)
z.root <- sqrt(z.sum / z.max) * 3
filled.contour(x=x, y=y, a, color.palette=terrain.colors,
plot.axes={axis(1); axis(2);
points(points, col="Black", cex=z.root)},
main="Proportion of Year at or above 5 m/sec",
xlab="Easting")
Best Answer
If you are using ArcGIS you would find these tools in the Spatial Analyst toolbar > Interpolation. You will input your featureclass containing the wind data, your Z value will be your wind speed, and then choose your output FC which is where you will save it.
Some choices have to be made as to your interpolation type. IDW is best for a dense sample set, SPLINE is best for smooth data. Quoted below is what ESRI has to say on the matter.
Here is a previous article on the matter: How do you decide what interpolation method to use for resampling raster data?
Here is a good paper on choosing the correct method: http://webapps.fundp.ac.be/geotp/SIG/interpolating.pdf