Covariance Matrix - Python

[cesardepaula]

I am implementing an algorithm in Python where I basically read 02 raster files, transform them to array I apply the function numpy.cov() to calculate the covariance matrix between these two raster files. Since I have 02 variables, namely two raster files of the same size, I expected the result of the covariance matrix to be another 2 x 2 matrix, but the result is a huge array whose values are not consistent with the expected. I will post below the code of the implemented sample and the result obtained in Python and I will also post the result of the covariance matrix obtained by the r.covar function of QGis. By my experience the value computed by QGis is the expected result for the covariance matrix.

Algorithm in Python:

from osgeo import gdal 
import numpy as np 

gdal.AllRegister() 

ds_r = gdal.Open('Sun_Glint_Sample_Red.tif') 

ds_n = gdal.Open('Sun_Glint_Sample_Nir.tif') 

rows = ds_r.RasterYSize 
print('Rows:', rows) 
cols = ds_r.RasterXSize 
print('Cols:', cols) 
bands = ds_r.RasterCount 
print('Bands:', bands) 


data_r = ds_r.ReadAsArray(0, 0, cols, rows) 

data_n = ds_n.ReadAsArray(0, 0, cols, rows) 

cov_r_n_arr = np.cov(data_r, data_n, bias=True) 

print('Size of Covariance Matrix calculated:', np.size(cov_r_n_arr)) 
print('Covariacne Matriz between the bands Red and Nir:', cov_r_n_arr)

Results from Covariance Matrix in Python:

Output:
Rows: 587 
Cols: 712 
Bands: 1 
Size of Covariance Matrix calculated: 1378276 
Covariacne Matriz between the bands Red and Nir: [[ 9.73841290e+02  5.59548015e+02  4.04565317e+01 ... -3.73723666e+01 
  -1.00562133e+01  5.41880898e+01] 
 [ 5.59548015e+02  8.18410331e+02  4.06466679e+02 ...  1.01076194e+01 
   1.23401342e+01  5.28321274e+00] 
 [ 4.04565317e+01  4.06466679e+02  6.39376657e+02 ...  1.98574438e+01 
   8.69184762e-01 -2.48161769e+01] 
 ... 
 [-3.73723666e+01  1.01076194e+01  1.98574438e+01 ...  2.37988588e+02 
   1.27523701e+02  1.51711728e+01] 
 [-1.00562133e+01  1.23401342e+01  8.69184762e-01 ...  1.27523701e+02 
   2.43866003e+02  1.18732639e+02] 
 [ 5.41880898e+01  5.28321274e+00 -2.48161769e+01 ...  1.51711728e+01 
   1.18732639e+02  2.08042748e+02]] 

Results from Covariance Matrix using QGis:

Output:
N = 417944 
1.000000 0.886645 
0.886645 1.000000

Thank you all very much and I look forward to helping you understand what I'm doing wrong in Python.

[scidam]

Try to flatten arrays first, e.g.

np.corrcoef(data_r.ravel(), data_n.ravel())

[cesardepaula]

Thank you scidam, now the results from covariance are the same results calculated using QGis.