N-Dim array manipulation in a loop, getting IndexError: too many indices for array

[cesardepaula]

Hello people, good morning!!

In my implementation I`m trying to manipulation N-Dim. array on a looping to calculate a covariance between 2 variables from this N-Dim. array.

To explain better, I have one raster file containing 4 bands (R, G, B, Nir). On my implementation I read raster file, extract rows, columns and number of bands. After that I read the raster file as array, according code below.

My goal on this algorithm is to calculate the covariance for R, G, B bands in relation Nir band in the same looping. For example:
Cov_1 = np.cov(R, Nir)

Cov_2 = np.cov(G, Nir)

Cov_3 = np.cov(B, Nir)

The algorithm has implemented in python is:

from osgeo import gdal
import numpy as np

ds = gdal.Open('Sun_Glint_Sample_.tif')

rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

cov = []

for i in range(bands):
    i += 1
    data_ds = ds.GetRasterBand(i).ReadAsArray(0, 0, cols, rows).ravel()
    cov = np.cov(data_ds[:, i], data_ds[:,4], bias = True)
    print('Dimension: ', data_ds.ndim)
    print('Shape: ', data_ds.shape)
    print('Array: ', data_ds[:, i])
    print('Covariance: ', cov)

Below the information printed before to calculate the covariance:

Output:
runfile('D:/CESAR_PHD/10.Lyzenga/Sun_glint_removal_v04.py', wdir='D:/CESAR_PHD/10.Lyzenga')
Dimension:  1
Shape:  (417944,)
Array:  [135 123  94 ...  31  57  77]
Dimension:  1
Shape:  (417944,)
Array:  [191 181 139 ...  49  93 128]
Dimension:  1
Shape:  (417944,)
Array:  [176 164 126 ...  42  78 107]
Dimension:  1
Shape:  (417944,)
Array:  [91 98 69 ... 12 19 33]

Below the error printed when I tried calculated the covariance:

Error:
  File "D:/CESAR_PHD/10.Lyzenga/Sun_glint_removal_v04.py", line 22, in <module>
    cov = np.cov(data_ds[:, i], data_ds[:,4], bias = True)

IndexError: too many indices for array

Thank you for help me!!

[scidam]

Try the following code:

nir_flattened = ds.GetRasterBand(bands - 1).ReadAsArray(0, 0, cols, rows).ravel()

for i in range(bands - 1):
    data_ds_flattened = ds.GetRasterBand(i).ReadAsArray(0, 0, cols, rows).ravel()
    cov = np.cov(data_ds_flattened, nir_flattened, bias=True)
    print('Covariance: ', cov) # expected cov is 2 x 2 matrix

But, I am not sure that: 1) bands counting starts from 0, 2) would it possible to read bands in arbitrary order, as I did; probably something like .seek method should be applied first.