6

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 '14 at 22:02
  • 1
    The means are identical. – Jot eN Sep 12 '14 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 '14 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 '14 at 16:04
  • upvoting the question... wish someone had an answer – J. Win. Oct 17 '15 at 20:44
5

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:

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 '14 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 '14 at 9:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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