实现对应像元的线性回归（demo：假设是以连续三年的两种遥感器的栅格影像（假设是2*2大小的）为例）

import numpy as np
from scipy import stats, linalg


#main函数


g = np.array([[[1, 3], [2, 8]], [[6, 2], [6, 9]], [[8, 1], [10, 11]]])   #第一种传感器
print("g:", g)
print("g第一个图像:", g[0])
print("g第一个图像第一行:", g[0, 0])
print("g第一个图像第一行第一列的元素:", g[0, 0, 0])
a = g.shape[0]
print("a(0维-图像数):", a)
b = g.shape[1]
print("b(1维-行数):", b)
c = g.shape[2]
print("c(2维-列数):", c)
n = b*c
print("n(一个图像的像元数)：", n)
m = np.array([[[2, 6], [4, 10]], [[7, 5], [8, 12]], [[9, 2], [12, 15]]])   #第二种传感器
print("m:", m)
print("m(第一个图像):", m[0])



#定义数组，为了把栅格数组转化为能做逐像元线性回归
x = np.array(list([0.0 for i in range(a)]for i in range(n)))
print("x初始:", x)
y = np.array(list([0.0 for i in range(a)]for i in range(n)))
print("y初始:", y)
lin = np.array(list([[[0.0 for i in range(c)]for i in range(b)]for i in range(5)]))
print("lin初始:", lin)

#把三维b行c列的数组，把对象中对应的第一个元素提取出来，一一对应转换为二维数组
l=0
for k in range(c):
    for j in range(b):
        #print("l:", l)
        for i in range(a):
            #x.append(g[i, j, k])
            print("j:", j)
            x[l, i] = g[i, k, j]
            print("gx:", x)
            y[l, i] = m[i, k, j]
            print("my:", y)
        l = l + 1

##对转换出的二维数组线性回归+把回归结果slope等，放在栅格三维数组里
    j = 0
    k = 0
for i in range(x.shape[0]):
    res = stats.linregress(x[i], y[i])
    slope, intercept, r_value, p_value, std_err = stats.linregress(x[i], y[i])
    print("slope:", slope)
    print("intercept:", intercept)
    print("r_value:", r_value)
    print("p_value:", p_value)
    print("std_err:", std_err)
    lin[0, j, k] = slope
    lin[1, j, k] = intercept
    lin[2, j, k] = r_value
    lin[3, j, k] = p_value
    lin[4, j, k] = std_err
    k = k+1
    if k > c-1:
        j = j+1
        k = 0
print("lin:", lin)