I have a 9 year weekly time series (~500 raster grids of equal cell size and extent). I'm interested in obtaining the regression line slope between pixels (Imagine stacking all 500 grids on top of one another and running a linear regression between each individual pixel). Essentially, this would measure the net change between pixels through my time series.

My main question is about the coefficients defined in the first answer to Making linear regression across multiple raster layers using ArcGIS Desktop?

Is the coefficient defined as 12 due to the temporal scale of the original question being at the monthly level and, if so, should my code's coefficient be set to 7 due to the weekly temporal scale between grids?

I have been referencing these Q&As that are similar:

I wrote this python code after reading the links referenced above:

from __future__ import division
import arcpy
import os
from arcpy import env
from arcpy.sa import*

# define workspace

# enable overwriting

# check spatial analyst extension

# define output paths

# list all rasters in the workspace
# sort rasters numerically

# get the number of rasters
print n

# setup index

# define division
print 'the global seed is {0}'.format(seed)

for raster in rasters:
    print i

    print 'Raster {0} is {1}:'.format(i,raster)
    print 'the coef for raster {0} is {1}'.format(i,coef)

    # Multiple raster by coefficient
    if i==1:
        i+=1  # same as saying i=i+1
        print 'adding {0} to outSlope'.format(raster)
    if i==6:

# Save final slope grid
print 'saving final slope grid'
outSlope.save(slopePath + "/" + "floodChange_py")

print 'script is complete'

I created this code as a test which appeared to work on a 5 week subset of my data. If what I have referenced above is correct then this code will work on a time series (of equal cell size and extent) of any length. Although there were no bugs in the code, I have a feeling I may have made a mistake someplace.

On a side note I'm using PyScripter, hence the "from future import division" in the first line. This disables the automatic floor function which rounds floating point numbers to the nearest integer away from zero.


2 Answers 2


Because we are free to choose how the times are designated, let's center them at zero with unit spacing, so that for n=500 times their values are -(n-1)/2, -(n-1)/2+1, ..., +(n-1)/2-1, +(n-1)/2. The ordinary least squares formula for the slope is obtained by multiplying each of these time values by their corresponding observations, taking the average of all those values, and dividing that average by the variance of the times. Because the average time is at zero, you can compute that the variance just the average squared time,

([-(n-1)/2]^2 + [-(n-1)/2+1]^2 +  ... + [(n-1)/2]^2) / n = (n^2-1)/12.

That is where the factors of 12 come in. You can also now see why the factor of 1/(n^3-n) in the first-referenced answer appears. Its reciprocal equals n(n^2-1): the factor of n comes from the averaging of the n products and the factor of n^2-1 comes from the variance of the times.

The units in which the slope is thereby expressed are the units of the observations divided by the time interval. For weekly time series, then, the slope gives the mean change in values per week.

A better way to assess the net change of values over the course of the times is to subtract the median of the values observed during the first few weeks from the median of values observed during the last few weeks. This is conceptually simple, it is easy to compute (use local statistics), it is fast (you only need to access a small subset of the datasets), it is robust (all kinds of weird behavior during the time period will not affect it), it is flexible (you can choose the durations of the two end periods), it more accurately captures the actual start-to-end-differences, and it is easily applied to intermediate times to assess shorter-term changes.


Here is an R solution from a NDVI timeseries thread. Although, I am with @whuber, there are more computationally efficient methods that are more interpretable. I like his idea of subtracting the median from the end points. Although I would want to see this done, systematically, through the time series in order to identify extremes that illustrate a departure from the overall trend. This seems analogous to Quantile regression approaches.

Often time series exhibit stochasticity that confound linear fits and make regression approaches difficult without a specialized model. If you are staying with regression you may want to consider a mixed effects ARIMA(I) model or a polynomial regression (LOWESS) smoother to account for serial autocorrelation and converge on a more stable regression slope.

An alternative approach, that you could explore, is the residual error using an Ordinary Least Squares (OLS) regression. You would actually want the serial autocorrelation in the residuals which, negates the need to account for serial autocorrelation as a random effect. There are several papers on Krigin residual error (AKA. Regression Kriging) based on spatial autocorrelation IID assumptions. You could draw similar inference on temporal residual error. However, keep in mind that the spatial autocorrelation would be confounded in the residual error as you are inherently dealing with a spatial-temporal model.

  • 1
    +1. But we might as well go all the way and admit that this situation is more than rich enough to support a full spatial-temporal model.
    – whuber
    Commented May 28, 2014 at 20:57
  • @whuber, I have always thought this the case with these large raster stacks. Since it seems to come up often, it would be nice to formalize some advice in dealing with this problem. Commented May 28, 2014 at 21:19
  • 1
    Current wisdom (and ongoing research) indicate that the right way to approach many of these situations begins with a physical model, often given as a partial differential equation or a discrete approximation thereof. ("Physical" is meant broadly; for instance, I am reading over a paper in which surveys of tuna abundance in fisheries are modeled this way.) That precludes offering any canonical one-size-fits-all answer (although I realize that's not exactly what you were suggesting)--but it's also what makes the study of spatially extensive phenomena so varied and interesting.
    – whuber
    Commented May 28, 2014 at 21:33
  • It does seem like the most tractable way to address spatial-temporal problems is to directly solve for a given hypothesis. You can, in theory, solve this type of problem using MCMC approaches (per Cressie, etc...) but if you break down the actual state-space it becomes intractable quite quickly. The higher the spatial and temporal dimensionality the more difficult it will be to converge. For point pattern problems, I have always thought that Scan Statistics are quite elegant for spatial-temporal data. Unfortunately, this methods cannot be adapted for this type of analysis. Commented May 28, 2014 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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