Lineregress numpy rasterio error

Question

I want to get the lineregress parameters of 2 numpy arrays (that both come from the same shape mask, so they should have the same shape) but I get this error message:

    ValueError                                Traceback (most recent call last)
<ipython-input-11-9eba54f5694d> in <module>()
----> 1 slope, intercept, r_value, p_value, std_err = linregress(a, b)

C:\Anaconda3\lib\site-packages\scipy\stats\_stats_mstats_common.py in linregress(x, y)
     77 
     78     # average sum of squares:
---> 79     ssxm, ssxym, ssyxm, ssym = np.cov(x, y, bias=1).flat
     80     r_num = ssxym
     81     r_den = np.sqrt(ssxm * ssym)

C:\Anaconda3\lib\site-packages\numpy\lib\function_base.py in cov(m, y, rowvar, bias, ddof, fweights, aweights)
   2493     else:
   2494         X_T = (X*w).T
-> 2495     c = dot(X, X_T.conj())
   2496     c *= 1. / np.float64(fact)
   2497     return c.squeeze()

ValueError: shapes (2,1471,769) and (769,1471,2) not aligned: 769 (dim 2) != 1471 (dim 1)

No idea of what is happening here, It seems that one matrix is flipped. I have been looking for something about this, but nothings works for me.

My code:

import rasterio
import numpy as np
from scipy.stats import linregress

stack1 = r'path/to/raster1.tif'
stack2 = r'path/to/raster2.tif'


with rasterio.open(stack1) as stck10:

    STCK1 = stck10.read()

with rasterio.open(stack2) as stck11:

    STCK11 = stck11.read() 

slope, intercept, r_value, p_value, std_err = linregress(STCK10, STCK11)

Stacks 1 and 2 comes from the same mask, so there should excatly the same size, shape, etc.

Answer 1

This error is not a result of Rasterio, it is because you are passing a 3-D array into the scipy.stats.linregress function. The x and y arguments are single sets of measurements (if both are specified). The axis reversal you are experiencing is a result of the transpose operation performed in the calculation of a covariance matrix. If you wish to get regression parameters for the entire set of data (from all bands) for both .tif files, you can flatten the resulting arrays like so:

slope, intercept, r_value, p_value, std_err = linregress(STCK10.ravel(), STCK11.ravel())

Otherwise, you can read individual bands at a time using STCK1 = stck10.read(band_index) or slice the array as needed to change your measurement sample.