I've calculated covariance between 7 rasters using R and ArcGIS. It gives me different results. Each times covariance values from R are approximately 1.874 times higher than ArcGIS values.

R - package raster, function layerStats(stack, 'cov', na.rm=TRUE) enter image description here

ArcGIS - Band Collection Statistics tool (http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000q3000000) enter image description here

Why covariance values are different in R and ArcGIS?

  • 2
    How do the means reported by ArcGIS compare to those reported by R?
    – whuber
    Sep 11, 2014 at 22:02
  • 1
    The means are identical.
    – Jot eN
    Sep 12, 2014 at 7:45
  • 1
    If na.rm=TRUE, N is equal to (n - cellStats(r, stat = "countNA") - asSample). How ArcGIS manages the missing values?
    – user32309
    Sep 12, 2014 at 10:49
  • What data type are the rasters (ints, floats, ... how many bytes)? How many cells do they have? Since the means are reported correctly, the problem is probably not related to NoData values, but it might still be of interest to know how many NoData values they have, too. Are there any cells where some of the rasters have NoData and others do have data?
    – whuber
    Sep 12, 2014 at 16:04
  • upvoting the question... wish someone had an answer
    – J. Win.
    Oct 17, 2015 at 20:44

1 Answer 1


Create some simple test rasters in R:

> m=matrix(1:9,3,3)
> m2 = matrix(c(9,2,3,4,1,5,6,8,7),3,3)

Then we can trivially compute the covariance between these matrices:

> cov(c(m),c(m2))
[1] 2.125

and I would wager doing the computation by hand would get the same result. What does layerStats do?

> D = stack(raster(m),raster(m2))

> layerStats(D, "cov", na.rm=TRUE)
        layer.1 layer.2
layer.1   7.500   2.125
layer.2   2.125   7.500

There's the same 2.125 in the cross-correlation.

Now try with an NA in there:

> m[2,2]=NA
> cov(c(m),c(m2),use="complete.obs")
[1] 2.428571
> D = stack(raster(m),raster(m2))
> layerStats(D, "cov", na.rm=TRUE)
         layer.1  layer.2
layer.1 8.571429 2.428571
layer.2 2.428571 7.500000

Again, cov agrees with layerStats if we remove the NA value.

So what does ArcGIS do in each of these situations? Save the D, import to ArcGIS, and find out...


The results from ArcGIS:

Raster m and m2:

enter image description here

Raster with an NA and m2:

enter image description here

[Now, the question is - Which way of calculating covariance between rasters is correct?]

  • 1
    +1 even though it's only half an answer, because it shows exactly the right way to approach the question: test the software on simple data with known results.
    – whuber
    Sep 12, 2014 at 15:58
  • The covariance and correlation are unchanged in ArcGIS even with a missing value? You've used the right second layer because I can see STD of layer 2 is sqrt(var(c(1,2,3,4,6,7,8,9)). How many missing values can you add to layer 2 before it changes? Does it just copy layer 1 to layer 2 if missing in layer 2? What does the docs say?
    – Spacedman
    Sep 16, 2014 at 9:20

Your Answer

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

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