Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I've been using ArcGIS for a couple years, but when it comes to model builder I've only worked with very simple models. I'm hoping someone that is much smarter than myself can help me in building a complicated model.

I've got a very large table of data (about 300k rows) that contains points with lat/lon data, wind speed and date/time. There are about 60 different points with wind values that have been recorded for each hour during 2012. What I'd like to do is build a model that will:

  1. Iterate through each row and select/export the data associated with first hour (00Z) Jan 1st 2012 (about 60 rows) to a separate table, then the 2nd hour (01Z) of Jan 1st 2012 to its own separate table and so on through all of 2012.
  2. Perform Inverse Distance Weighting on all of the hourly tables with lat/lon and wind values that were created from step 1 above
  3. Reclassify the rasters so that points in the raster that have wind values of 5 or greater are reclassified to "1" and all other cell's wind values are reclassified to "0"
  4. Then iterate through the rasters created in step 3 and use something like "Collect Values" & "Cell Statistics" to sum up all of the rasters into a final summed raster for all hours in 2012.

Was going to include screenshot of my data, but I guess since this is my first post I'm not trusted yet. I've broken the date down several different ways (just date, just time, date-time, and date-time in text format) because I've had issues trying to iterate through and select the data in my original, very large table only by date.

Hopefully I provided enough information, but if not I will try to go into more detail. Thanks in advance for any assistance.

share|improve this question
    
To be honest, you may want to consider a Python script that will loop through everything. Model builder will be cumbersome with nearly 52 iterations. It's possible, but there would be a LOT of copy/pasting and setting up each and every tool. How familiar are you with Python? It's also not very clear what you need on step 4. Are you wanting nearly 52 raster summaries to be produced in the end? Or a TOTAL of ALL rasters in ONE raster? –  Spatial Pariah Apr 22 '13 at 22:10
    
@angelspatial By what process do you arrive at 52 rasters? 2011 had 8,760 hours in it. –  whuber Apr 22 '13 at 22:22
    
@whuber wow, you're right. For some reason I counted it as weeks, not hours. Yeah, don't use model builder for this - ever. –  Spatial Pariah Apr 22 '13 at 22:24
    
Right, there would be many, many rasters by the end of step 3...roughly 8,760. I opted for model builder because I'm even less familiar with writing in python. I'd be willing to try, but would probably need a lot of assistance. –  Daniel Apr 22 '13 at 22:27
    
How large do you need your output raster to be, Daniel? (How many rows and columns?) Also, do you have complete data for all 60 locations or are the data occasionally missing at some of them? Incidentally, you don't need to have 8760 rasters at the end: you will have the last raster you computed (which you can delete) together with the cumulative values from all 8760 of them: just two rasters all told. –  whuber Apr 22 '13 at 22:29

1 Answer 1

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.

Plot

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")
share|improve this answer
    
I like your answer better. –  Spatial Pariah Apr 23 '13 at 14:59
    
Thank you very much, whuber, for taking the time to give this very detailed response. This is very helpful and something I can try and experiment with. Let me give a little more background into what I'm trying to do: I'm trying to create a methodology of displaying the frequency of favorable wind conditions at different locations if I were to have archived gridded wind model data for a year. Right now, the wind data at the resolution I need is not archived so it is not available. So, I'm trying to come up with a plan of attack assuming this data was available. ..cont'd.. –  Daniel Apr 23 '13 at 16:46
    
..In order to best simulate this data, I wanted to create the rasters of wind model data for 2012 by using observed wind data and putting it in a table, then using some form of interpolation (wanted Kriging, but have had issues normalizing the data) and then reclassify them, sum them, and display the final results as percentage of favorable wind conditions at each location (total favorably modeled wind conditions in 2012 divided by total modeled wind conditions in 2012) –  Daniel Apr 23 '13 at 16:46
    
But it might be that this method is not the best and I need to rethink it. If the process of creating the 8,760 rasters takes several hours (at best) to complete I would be okay with it (whether by model or Python). It might be troublesome if the process were going to take more than a day or two though. In retrospect I should have been more clear about what I'm trying to do in my original question, so I apologize. –  Daniel Apr 23 '13 at 16:57
    
It's clear what you're trying to do, Daniel: this is precisely the kind of information needed for designing wind farms (because turbines need minimum windspeeds in order to operate). Windspeeds tend to have skewed distributions (which is why I simulated them from Gamma distributions), but that can be handled in kriging via a non-linear rexpression, such as the square root, or with a GLM. Of more import is the likelihood that winds are controlled by local topography, perhaps changing abruptly in speed and direction at ridge lines. Kriging can include this kind of information; IDW can't. –  whuber Apr 23 '13 at 17:10

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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