Python: Как найти уравнение регрессии нескольких 3D-массивов данных (широтно-временное значение)?

Question

У меня есть два набора dataArrays, представляющих значение по трем координатным осям (широта, долгота, время); один набор dataArrays представляет переменную varA, один набор представляет переменную varB (пример приведен ниже).

varA
<xarray.DataArray 'varA' (time: 32, lat: 20, lon: 18)>
array([[[... 
... ]]])
Coordinates:
 * lat      (lat) float64 4.75 4.25 3.75 3.25 2.75 2.25 1.75 1.25 0.75 0.25 ...
 * lon      (lon) float64 33.25 33.75 34.25 34.75 35.25 35.75 36.25 36.75 ...
 * time     (time) datetime64[ns] 1979-01-01 1980-01-01 1981-01-01 ...

Есть 4 разных дуэта varA и varB. Я сделал график рассеяния, который показывает varA против varB для каждого широты, долготы и времени, и это для каждого дуэта (varA1 против varB1, varA2 против varB2...).; и визуализировал его на 1 фигуре (все массивы данных имеют одинаковые координаты).

for t in range(varA1['time'].size) :
    for la in range(varA1['lat'].size) :
        for lo in range(varA1['lon'].size) :
            x = varA1.values[t,la,lo]
            y = varB1.values[t,la,lo]
            plt.scatter(x,y)
for t in range(varA2['time'].size) :
    for la in range(varA2['lat'].size) :
        for lo in range(varA2['lon'].size) :
            x = varA2.values[t,la,lo]
            y = varB2.values[t,la,lo]
            plt.scatter(x,y)
... 
plt.show()

Функция полного графика рассеяния работает нормально, но теперь я хотел бы добавить линию тренда (и найти ее уравнение) для полного графика рассеяния. На самом деле, я хочу исследовать долгосрочную (со временем) и крупномасштабную (по широте и долготе) связь между varA и varB; Я знаю, что более высокие значения varA сопровождаются более высокими значениями varB (независимо от того, в каком месте или в какое время), однако я хочу получить 1 уравнение регрессии, коэффициент корреляции и среднеквадратическое отклонение (т.е. количественно определить их соотношение).

Это возможно? Я думаю, что мне нужно среднее значение всех varAs на ячейку [широта, долгота, время] и то же самое для всех varBs; так как то так:

avrR = wfdei_rain * 0
avrY = wfdei_rain * 0
for t in range(varA1['time'].size) :
    for la in range(varA1['lat'].size) :
        for lo in range(varA1['lon'].size) :
            avrA[la,lo,t] = float(sum([varA1[la,lo,t],varA2[la,lo,t],varA3[la,lo,t],varA4[la,lo,t])) / 4    
            avrB[la,lo,t] = float(sum([varB1[la,lo,t],varB2[la,lo,t],varB3[la,lo,t],varB4[la,lo,t])) / 4                                
z = np.polyfit(avrA[:,:,:],avrB[:,:,:],1)
p = np.poly1d(z)
plt.plot(x,p(x))
print('y=%.6fx+(%.6f)'%(z[0],z[1])

Это дает ошибку в функции polyfit ("ожидаемый 1D вектор для x"). Я не нахожу, как настроить polyfit для работы с DataArrays - 3D. Я нашел похожую запись: применение numpy.polyfit к xarray Dataset, но это не помогло мне.

Любые предложения о (лучшем) способе подойти к этому?

Заранее спасибо!