Python：如何沿轴执行两个 numpy 3D 数据集的线性回归？

Question

我有两个特定区域的数据集：第一个是降雨量，第二个是该区域的植被测量 (npp)。因此，前两个维度 (x,y) 表示地理位置。第三个维度是时间（8 个时间步长）。我想要做的是对 8 个降雨值与 8 个植被值的每个位置进行线性回归。结果应该是几个二维数组，其中为每个位置计算 p 值、r²、斜率和理想情况下的残差，或者所有值都在 3D 数组中。

nppList = glob.glob(nppPath+"*.img")
rainList = glob.glob(rainPath+"*.img")

nppImg = [gdal.Open(i) for i in nppList]
rainImg = [gdal.Open(i) for i in rainList]

nppFiles = [i.ReadAsArray() for i in nppImg]
rainFiles = [i.ReadAsArray() for i in rainImg]

# get nodata
nppNodata = nppImg[1].GetRasterBand(1).GetNoDataValue()
rainNodata = rainImg[1].GetRasterBand(1).GetNoDataValue()

# convert to float and set no data
nppStack = nppStack.astype(float)
nppStack[nppStack == nppNodata] = np.nan
rainStack = rainStack.astype(float)
rainStack[rainStack == rainNodata] = np.nan

# instead of range(0,8) there should be the rainfall variable, but on a pixel base
def linReg(a):
    return stats.linregress(a, range(0, 8))

lm = np.apply_along_axis(linReg, axis=2, arr=nppStack)

我知道函数 numpy.apply_along_axis() 但这里一个函数只能应用于一个数组。我正在寻找一种可能性，可以沿一个轴在两个数组上应用函数，最好不要循环遍历数组。

Answer 1

source for scipy.stats.linregress 表示仅不支持维度大于 2 的数组（仅适用于您的 x 和 y 数据恰好位于同一数据结构中的情况）。

老实说，在您的情况下，我会使用 Python 循环——代码中最慢的部分不太可能在数据点上循环；相反，回归本身将决定速度。

在这种情况下，您可以展平位置轴，使用单个循环，然后将回归结果重新整形为 3D。比如：

n = nx * ny
frain = rainStack.reshape((n, 8))
fnpp = nppStack.reshape((n, 8))
reg_results = np.empty((n,5))
for i in range(n):
    reg_results[i] = stats.linregress(frain[i], fnpp[i])
reg_results[i].reshape((nx,ny,8))    # back to 3D