栅格的二元回归

有时候我们会用到残差趋势法，例如以植被覆盖度为因变量 、以气温和降水为自变量，逐像元建立二元线性回归模型 ，逐像元得到回归方程的系数；其次，利用气温和降水数据以及回归模型的系数，建立模型模拟得到气候影响下的植被覆盖度的预测值；最后，基于遥感影像获得的植被覆盖度观测值与基于回归模型模拟得到气候影响下的预测值做差值计算，得到的结果即为植被覆盖度残差，表示了人类活动对植被覆盖的影响。今天分享一下栅格的二元回归系数计算方法。

1 二元回归系数计算

二元回归系数系数的计算公式网上书上有计算公式，这里不再赘述。这里介绍一下Python的numpy库计算相关系数，使用sklearn.linear_model.LinearRegression类，示例如下。 下面的示例中，首先构造一个x1和x2数组，每个数组有6个元素，其次构造y = a·x1 + b·x2 + c。其中a = 3.14，b = 2.78, c = 1.27，使用LinearRegression来反演a b c，与给定的的系数abc是一致的，验证了LinearRegression进行二元回归的可行性。

import numpy as np
import sklearn
from sklearn.linear_model import LinearRegression

x1 = np.array([1,2,3,4,5,6])
x2 = np.array([3,2,1,3,3,2])

# 让a = 3.14， b = 2.78, c = 1.27构造一个y

y = 3.14 * x1 + 2.78 * x2 + 1.27 
x = np.vstack([x1, x2]).T
regr = sklearn.linear_model.LinearRegression()
regr.fit(x, y)
a, b = regr.coef_
c = regr.intercept_
print(a,b,c)
# 3.1399999999999997 2.78 1.2699999999999996

2 计算栅格的二元回归

2.1 数据准备

以植被覆盖度、降水和温度的计算为例，首先准备植被覆盖度、气温和降水数据，例如10个月的植被覆盖度、10个月的气温和10个月的降水。回归计算需要多个数据，只用一期的影像无法计算。

待计算的数据需要有相同的行列数和相同的期数。
通常植被覆盖度、降水、温度的数据是单波段的，首先需要将单波段的影像按照相同的时间顺序合并成一个多波段的，例如一个月一个波段或一年一个波段，具体看使用的数据。这个可以使用envi的layer stacking来实现，也可以使用下面的代码实现。在这里的示例中，3月份为第1波段，4月份为第2波段。。12月份为第10波段。 气温和降水分别做成这样的数据。

from osgeo import gdal
import os
"""
多个单波段tif合并成一个tif文件
"""
#修改路径
tifDir = r"E:\data\一元回归\降水"  #tif路径 单波段
outtif = r"E:\data\一元回归\降水.tif"

NP2GDAL_CONVERSION = {
  "uint8": 1,
  "int8": 1,
  "uint16": 2,
  "int16": 3,
  "uint32": 4,
  "int32": 5,
  "float32": 6,
  "float64": 7,
  "complex64": 10,
  "complex128": 11,
}
tifs = [i for i in os.listdir(tifDir) if i.endswith(".tif")]
#获取投影波段数等信息
bandsNum = len(tifs)
dataset = gdal.Open(os.path.join(tifDir,tifs[0]))
projinfo = dataset.GetProjection()
geotransform = dataset.GetGeoTransform()
cols,rows=dataset.RasterXSize,dataset.RasterYSize
datatype=dataset.GetRasterBand(1).ReadAsArray(0,0,1,1).dtype.name
gdaltype=NP2GDAL_CONVERSION[datatype]
dataset=None
#创建目标文件
format = "GTiff" #tif格式
#format = "ENVI"  # ENVI格式
driver = gdal.GetDriverByName(format)
dst_ds = driver.Create(outtif,cols, rows,bandsNum, gdaltype,options=['COMPRESS=LZW'])
dst_ds.SetGeoTransform(geotransform)
dst_ds.SetProjection(projinfo)
#写入文件
info = set()
for k in range(bandsNum):
    ds = gdal.Open(os.path.join(tifDir,tifs[k]))
    X,Y = ds.RasterXSize,ds.RasterYSize
    info.add("%s,%s"%(X,Y))
    if(len(info) != 1):
        dst_ds = None
        ds = None
        print("%s 列数行数不一样：%s,%s"%(tifs[k],X,Y))
        raise Exception("有影像行列数不一致")
    data = ds.GetRasterBand(1).ReadAsArray()    ##读取第一波段
    ds = None
    dst_ds.GetRasterBand(k+1).WriteArray(data)
    dst_ds.GetRasterBand(k+1).SetDescription("hahha_%s"%k)
    print("波段 %s ==> 文件 %s"%(k+1,tifs[k]))
dst_ds = None

2.2 计算

数据准备好后就可以使用Python写代码来计算了。需要使用numpy、gdal、sklearn等库来操作。代码在下面贴出，注释十分的丰富，很利于理解。只需要在最后的main函数中修改输入输出路径即可。结果tif有3个波段，分别是temp系数a、rain的系数b和常数项c。

import numpy as np
from osgeo import gdal
import os
import sklearn
from sklearn.linear_model import LinearRegression
from tqdm import tqdm


def regression(X, Y):
    '''
    X为n行2列的numpy数组，n为年份数，第一列为气温，第二列为降水
    Y为目标：n个元素的numpy数组，为fvc
    '''
    regr = sklearn.linear_model.LinearRegression()
    regr.fit(X, Y)
    a, b = regr.coef_
    c = regr.intercept_
    return a, b, c


def erhuigui(fvcpath, temppath, rainpath, outtif):
    dataset = gdal.Open(fvcpath)
    projinfo = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    depth = dataset.RasterCount
    rows = int(dataset.RasterYSize)
    colmns = int(dataset.RasterXSize)

    fvc = dataset.ReadAsArray()
    src2 = gdal.Open(temppath)
    temp = src2.ReadAsArray()

    src3 = gdal.Open(rainpath)

    rain = src3.ReadAsArray()

    dst_ds1 = driver.Create(outtif, colmns, rows, 3,
                            gdal.GDT_Float32)
    dst_ds1.SetGeoTransform(geotransform)
    dst_ds1.SetProjection(projinfo)
    d0 = fvc[0]
    outA = d0 * 0 - 2222.0
    outB = d0 * 0 - 2222.0
    outC = d0 * 0 - 2222.0
    # 开始计算
    for row in tqdm(range(rows)):
        for col in range(colmns):
            if d0[row, col] != 0:
                y = fvc[:, row, col] * 1.0
                x1 = temp[:, row, col] * 1.0
                x2 = rain[:, row, col] * 1.0
                x = np.vstack([x1, x2]).T
                a, b, c = -2222, -2222, -2222
                try:
                    a, b, c = regression(x, y)
                except:
                    pass
                outA[row, col], outB[row, col], outC[row, col] = a, b, c
    print("a---------------")
    dst_ds1.GetRasterBand(1).WriteArray(outA)
    dst_ds1.GetRasterBand(2).WriteArray(outB)
    dst_ds1.GetRasterBand(3).WriteArray(outC)
	dst_ds1.GetRasterBand(1).SetNoDataValue(-2222)
    dst_ds1.GetRasterBand(2).SetNoDataValue(-2222)
    dst_ds1.GetRasterBand(3).SetNoDataValue(-2222)
    dst_ds1 = None


if __name__ == "__main__":
    # 修改路径
    fvcpath = r"R01.tif"  #植被覆盖度路径
    temppath = r"R02.tif"  #气温路径
    rainpath = r"R03.tif"  #降水路径
    outtif = r"out.tif"
    erhuigui(fvcpath, temppath, rainpath, outtif)

2.3 结果示例

3 结束语

好朋友们，以上就是今天分享的全部内容了，如有疑问和错误欢迎批评指正。