# Why covariance values are different in R and ArcGIS?

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)`

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

Why covariance values are different in R and ArcGIS?

• How do the means reported by ArcGIS compare to those reported by `R`? Sep 11, 2014 at 22:02
• The means are identical. Sep 12, 2014 at 7:45
• 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? Sep 12, 2014 at 16:04
• upvoting the question... wish someone had an answer Oct 17, 2015 at 20:44

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)
\$covariance
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)
\$covariance
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...

EDIT

The results from ArcGIS:

Raster m and m2:

Raster with an NA and m2:

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

• +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. 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? Sep 16, 2014 at 9:20