How to estimate the standard error on the mean difference between two DEMs?

Question

I have two DEMs of the same area measured on two different dates. Differencing these DEMs gives a map of the change in height between these two dates. I would like to calculate the mean change in height and the standard error on the mean over a region spanning a large number of pixels.

Calculating the mean is straight forward, but the standard error seems more complicated in that the noise is correlated between neighbouring pixels - as a result, just dividing the standard deviation of the height change in my region of interest by the root of the number of pixels gives an unrealistically low standard error.

Is there a way to calculate the standard error that takes into account this spatial autocorrelation?

Perhaps not the best route, but one you could take, would be to estimate the spatial autocorrelation length, then sample your subtracted dem at that length (or some factor e.g. 2x of it), computing your statistics on that sampled set. — Jon, Apr 24 '18 at 17:09
I am assuming that you want a local assessment (say, 5x5 focal window) of the correlated error, correct? There are methods that account for autocorrelation in the correlation, by correcting degrees of freedom. However, I am not aware of a method, or really a need, in correcting something like a RMSE for autocorrelation. You could apply a regression with correlated errors and examine the resulting RMSE but, this really seems like overkill. The issue may be that the moments are being effect by nonlinearity in the data (eg., highly skewed distributions in [x,y]). — Jeffrey Evans, Apr 24 '18 at 17:30
@JeffreyEvans Thanks for your comment, I'm trying to estimate the mean and standard error on the mean for an area covering ~16,000 pixels of my height difference raster, so it's quite a large area. I was imagining it might be possible to estimate an 'effective number of pixels', something like this: https://www.degruyter.com/view/j/mms.2011.xviii.issue-4/v10178-011-0052-x/v10178-011-0052-x.xml but in two dimensions instead of just one, which I guess is kind of what jon is suggesting — user7821537, Apr 24 '18 at 18:01
The paper that you are referencing is using the ACF which is a statistic used in serial autocorrelation. This would not equate to a spatial structure as it is sequential in nature and does not account for the matrix-neighbor characteristics in spatial autocorrelation. — Jeffrey Evans, Apr 24 '18 at 18:26

Jeffrey Evans · Answer 1 · 2018-04-24T18:39:49.677

You could most certainly examine the autocorrelated structure of the correlations using the Modified t-test, along with the F-statistic and p-values. A starting point may be the using the raster.modifed.ttest. Details on this R function are provided in this post. I have to say that I do not see where the standard error around the mean difference is going to be at all informative. A more correct statistic would be the root mean squared error (RMSE) or variants thereof.

Here is a method for calculating the local RMSE in R. It uses a spdep neighbor function so, is less than ideal on large rasters. The required object class is: SpatialPixelsDataFrame, SpatialGirdDataFrame or SpatialPointsDataFrame. You can get the correct class by reading your raster using rgdal::readRDGL or using the raster::raster function and then coercing using as(x, "SpatialPixelsDataFrame")

First we need to add required libraries and create some data.

library(gstat)
library(sp)
library(raster)
data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y

v <- fit.variogram(variogram(log(elev) ~ 1, meuse), vgm(.1, "Sph", 1000, .6))
elev <- krige(elev ~ 1, meuse, meuse.grid, v, nmax = 30)
  gridded(elev) <- TRUE
  names(elev@data) <- c("elev","var")  
elev2 <- elev
  elev2@data$elev <- elev2@data$elev + runif(nrow(elev), -3, 3) 

plot(stack(raster(elev,1),raster(elev2,1)))

Now we can specify a RMSE and se functions. If you just want the RMSE and standard error of the differences across the entire extent it is as easy as.

rmse <- function(x, y) sqrt(mean((x - y)^2, na.rm=TRUE))
se <- function(x) { sd(x, na.rm=TRUE) / sqrt(length(x[!is.na(x)])) }

rmse(elev$elev, elev2$elev  ) # global RMSE
se(elev$elev - elev2$elev)    # global Standard Error

Here is where we calculate the neighbors within the specified distance and loop through them to derive the local statistic(s).

d <- 200  
nb <- spdep::dnearneigh(sp::coordinates(elev), 0, d)
  xy.rmse <- rep(NA,nrow(elev))
  xy.se <- rep(NA,nrow(elev))
    for (i in 1:length(nb)) {
      x.var <- elev@data[nb[[i]], ][1][, 1]
      y.var <- elev2@data[nb[[i]], ][1][, 1]
      xy.rmse[i] <- rmse(x.var,y.var)
      xy.se[i] <- se(x.var-y.var)     
    }

We can add the local RMSE & SE to the original elev raster and plot the results.

elev$rmse <- xy.rmse
elev$se <- xy.se
plot(raster(elev,3),main="RMSE (d=200)")  
plot(raster(elev,4),main="Standard Error (d=200)")

How to estimate the standard error on the mean difference between two DEMs?

1 Answers1