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 have a time series of rasters where the cells of each raster represents some value at a certain time.

I'd like to generate a map that highlights the magnitude of positive or negative trends in this value over time for each cell.

My somewhat naive approach is to fit a simple linear regression (X=time and Y=value) to each cell and output the array of slopes to a raster (as per example images below). This can be filtered by only exporting significant results.

time series from four individual cells

Slopes from linear regression

How else might I represent trend over time in a raster timeseries?

NB: I'm interested in general techniques not software specific instructions.

share|improve this question
1  
Regression is a good idea. You can try also (last value - first value) / years, resulting in some Y units per year. –  nadya Feb 21 '13 at 4:12
1  
@Luke Have you had any success with regression method? Or any other approach? I've started a bounty on your questions since I'm also interested in such visualization. I'll award the bounty for the answer that will be best suited for polygon data though. –  radek May 7 '13 at 20:22
3  
Regression works, the final image in my question shows the slope of the fitted line from a simple linear regression. –  Luke May 7 '13 at 21:37
add comment

3 Answers

up vote 8 down vote accepted
+50

Plotting the estimated slopes, as in the question, is a great thing to do. Rather than filtering by significance, though--or in conjunction with it--why not map out some measure of how well each regression fits the data? For this, the mean squared error of the regression is readily interpreted and meaningful.

As an example, the R code below generates a time series of 11 rasters, performs the regressions, and displays the results in three ways: on the bottom row, as separate grids of estimated slopes and mean squared errors; on the top row, as the overlay of those grids together with the true underlying slopes (which in practice you will never have, but is afforded by the computer simulation for comparison). The overlay, because it uses color for one variable (estimated slope) and lightness for another (MSE), is not easy to interpret in this particular example, but together with the separate maps on the bottom row may be useful and interesting.

Maps

(Please ignore the overlapped legends on the overlay. Note, too, that the color scheme for the "True slopes" map is not quite the same as that for the maps of estimated slopes: random error causes some of the estimated slopes to span a more extreme range than the true slopes. This is a general phenomenon related to regression toward the mean.)

BTW, this is not the most efficient way to do a large number of regressions for the same set of times: instead, the projection matrix can be precomputed and applied to each "stack" of pixels more rapidly than recomputing it for each regression. But that doesn't matter for this small illustration.


# Specify the extent in space and time.
#
n.row <- 60; n.col <- 100; n.time <- 11
#
# Generate data.
#
set.seed(17)
sd.err <- outer(1:n.row, 1:n.col, function(x,y) 5 * ((1/2 - y/n.col)^2 + (1/2 - x/n.row)^2))
e <- array(rnorm(n.row * n.col * n.time, sd=sd.err), dim=c(n.row, n.col, n.time))
beta.1 <- outer(1:n.row, 1:n.col, function(x,y) sin((x/n.row)^2 - (y/n.col)^3)*5) / n.time
beta.0 <- outer(1:n.row, 1:n.col, function(x,y) atan2(y, n.col-x))
times <- 1:n.time
y <- array(outer(as.vector(beta.1), times) + as.vector(beta.0), 
       dim=c(n.row, n.col, n.time)) + e
#
# Perform the regressions.
#
regress <- function(y) {
  fit <- lm(y ~ times)
  return(c(fit$coeff[2], summary(fit)$sigma))
}
system.time(b <- apply(y, c(1,2), regress))
#
# Plot the results.
#
library(raster)
plot.raster <- function(x, ...) plot(raster(x, xmx=n.col, ymx=n.row), ...)
par(mfrow=c(2,2))
plot.raster(b[1,,], main="Slopes with errors")
plot.raster(b[2,,], add=TRUE, alpha=.5, col=gray(255:0/256))
plot.raster(beta.1, main="True slopes")
plot.raster(b[1,,], main="Estimated slopes")
plot.raster(b[2,,], main="Mean squared errors", col=gray(255:0/256))
share|improve this answer
    
@whubber Great answer (shall I say - as usual! :). Thank you. –  radek May 14 '13 at 7:57
    
Thanks @whuber. Outputting MSE required just a single line of code to be changed in my numpy script. –  Luke May 20 '13 at 8:29
add comment

What you are describing is "Change Detection". There are many techniques for change detection using rasters. Probably the most common is image differencing where you subtract one image from another to produce a third. Though, it depends upon the type of data you are trying to compare. From your image, it looks like you're comparing changes in slope over time (unless this area is subject to major land works, this isn't likely to change much). However, if you are comparing land class changes over time, you might use a different approach.

I came across this article by D. Lu et al. in which they compare different methods of change detection. Here's the abstract:

Timely and accurate change detection of Earth’s surface features is extremely important for understanding relationships and interactions between human and natural phenomena in order to promote better decision making. Remote sensing data are primary sources extensively used for change detection in recent decades. Many change detection techniques have been developed. This paper summarizes and reviews these techniques. Previous literature has shown that image differencing, principal component analysis and post-classification comparison are the most common methods used for change detection. In recent years, spectral mixture analysis, artificial neural networks and integration of geographical information system and remote sensing data have become important techniques for change detection applications. Different change detection algorithms have their own merits and no single approach is optimal and applicable to all cases. In practice, different algorithms are often compared to find the best change detection results for a specific application. Research of change detection techniques is still an active topic and new techniques are needed to effectively use the increasingly diverse and complex remotely sensed data available or projected to be soon available from satellite and airborne sensors. This paper is a comprehensive exploration of all the major change detection approaches implemented as found in the literature.

share|improve this answer
1  
I've upvoted your answer, however that paper seems to primarily review bi-temporal change detection (of which image differencing is a form) and the multitemporal methods it does cover make it difficult to (easily) identify direction/magnitude/significance of any change. My concern with bi-temporal change detection is that it does not incorporate all of the information contained in a time-series. For example, a large change from first to last year might mean nothing if the value varies considerably between years anyway. –  Luke Feb 26 '13 at 2:53
add comment

Another approach would be to calculate per cell statistics from multiple rasters. Specifically, use Cell Statistics (Spatial Analyst) to make these calculations. I think the most appropriate overlay statistic for your analysis would be to use the "range" statistic. The result will be a raster based map showing the difference between the largest and smallest values within each cell.

Cell Statistics - Range

enter image description here

share|improve this answer
2  
I've upvoted your answer, however using range alone does not allow me to identify direction (positive/negative) of any change. This could be worked around using the "HighestPosition" and "LowestPosition" tools though. My concern with this approach though is that it does not incorporate all of the information contained in a time-series. For example, a large change from min to max value might mean nothing if the value varies considerably between years anyway, though this could checked using the Cell Statistics tool with the standard deviation statistic. –  Luke Apr 22 '13 at 2:28
    
A hybrid approach may be worth looking into. For each time series, calculate the Max or Majority and Min or Minority. Use image differencing (Raster Calculator) between the two time series as @Fezter mentioned to determine the trends over time. I agree that it would be appropriate to also report error (st dev) with your time series map. –  Aaron Apr 22 '13 at 12:33
    
None of the approaches you suggest, Aaron, actually take time into account, so how could their results be construed as any objective measure of change over time? –  whuber May 7 '13 at 22:18
    
@whuber, let's assume there are ten raster datasets. If each raster dataset represents values at a corresponding time period, and cell statistics (e.g. range) are calculated on the collection of rasters, how then, is time not being taken into account? –  Aaron May 7 '13 at 23:29
2  
Because the statistics you propose carry no information about the sequence of values. For instance, the three sequences 1, 2, ..., 10; 10, 9, ..., 1; and 1, 10, 2, 9, 3, 8, 4, 7, 5, 6 all have identical cell statistics but represent the full gamut of possible trends over time, from straight up to straight down to no trend at all. None of your proposals distinguishes among these. –  whuber May 8 '13 at 12:45
add comment

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.