Want cell linear regression values for a netCDF or multi-band raster

I have climate raster type netCDF data.

I want the slope and other regression statistics of temperature for each cell over time.

I have read that GEOV IDL can do a temporal regression to achieve what I want but it is not open source. Does anyone know of a good open source software that can accomplish what I would like? R maybe?

Here is how you could get the slope, using R and the raster package. To (also) get the intercept see help(calc)

library(raster)
# b <- brick("file.nc")
# example data:
b <- brick(system.file("external/rlogo.grd", package="raster"))

# here time is 1 to n, but you can set it something else
time <- 1:nlayers(b)
# write a function that reruns the value or values of interest
fun <- function(x) { lm(x ~ time)\$coefficients }
x <- calc(b, fun)
plot(x)

If you have NA values, you need a more complex function

fun2 <- function(x) {
d <- na.omit(cbind(x, time))
if (nrow(d) > 2) {
lm(x ~ time, data=data.frame(d))\$coefficients
} else {
NA
}
}
b[1:10] <- NA
x2 <- calc(b, fun2)

to get r^2

fun <- function(x) { summary(lm(x ~ time))\$r.squared }

I noticed that this function is a bit slow as, for each raster cell, it fits a model and returns a lot of information that is not used. For each model some computations are repeated (as the independent variable is fixed). If you want simple output like the slope or intercept you can easily shortcut things by directly computing these via linear algebra, and pre-computing some of the intermediate (constant) results.

For the case without NAs only:

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
time <- 1:nlayers(b)

LMfun <- function(x) { lm(x ~ time)\$coefficients }
system.time( xlm <- calc(b, LMfun) )

# user  system elapsed
# 7.95    0.00    7.96

# add 1 for a model with an intercept
X <- cbind(1, time)
# pre-computing constant part of least squares
invXtX <- solve(t(X) %*% X) %*% t(X)
# much reduced regression model.  is to get the slope
LAfun <- function(y) (invXtX %*% y)
system.time( xla <- calc(b, LAfun) )

# user  system elapsed
# 0.06    0.00    0.06

So this approach is about 130 times faster!

If you're familiar with Python you can use the netCDF4-python library that can read and write both netCDF 3 and 4 data to numpy arrays. For example:

from netCDF4 import Dataset

root_group = Dataset("path_to_dataset", format='NETCDF4')
print root_group  \$ netCDF4 style dump

data = root_group.variables["some_variable"][:]

Python has a large number of other libraries which you can use for regression. Personally my first stop would be scikit-learn but you should also try statsmodels (especially if you're familiar with R - the models are defined in a very similar way).

I'm not as familiar with R, but there is certainly a ncdf4 package of CRAN and you can use the lm method to fit linear models easily for example. Also there are many tutorials available.

• Thanks for your explaination. This looks pretty straight forward, definitely will give it a try. I am much more familiar with python than R. – Megan Apr 28 '15 at 16:33
• This does not really answer the question – Robert Hijmans Apr 30 '15 at 0:34

GDAL has support for netCDF files:

gdalinfo --formats|grep -i cdf
GMT (rw): GMT NetCDF Grid Format
netCDF (rw+): Network Common Data Format

So, you can open this kind of files directly in QGIS. For this page:

http://www.unidata.ucar.edu/software/netcdf/examples/files.html

I downloaded this ECMWF_ERA-40_subset.nc climate sample file. It has 17 multiband layers (each layer with 62 bands: 1054 in total) and they were totally loaded in QGIS. I rendered the p2t layer ("2 metre temperature") as pseudocolor (5 classes) and the result was: It can be observed that the image is centered in the Pacific Ocean with a bound rectangle of: 0 90 357.5 -90 (See the metadata file). This is an important data to assign a projection by using gdalwarp in QGIS (Raster -> Projections -> Warp).

Imported to GRASS, you can calculate linear regression from two raster maps (y = a + b*x) with r.regression.line.

On the other hand, with the help of Value Tool Plugin in QGIS, you can do an exploratory previous analysis of all datasets and bands: • Thanks, this is good information to have. I am currently working in ArcGIS so I can view the files and do analysis between two isolated bands/raster layers quite easily. It was the regression values between all the bands that I specifically want. – Megan Apr 28 '15 at 16:26