计算两个栅格的相关系数

有时候我们需要计算两个栅格的相关系数，判断相关性，例如计算NDVI和降水的相关系数，NDVI和温度的相关系数。今天分享一下计算两个栅格相关系数的计算方法。

1 相关系数计算

相关系数的计算公式网上书上有计算公式，这里不再赘述。这里介绍一下Python的numpy库计算相关系数，使用np.corrcoef()函数，示例如下。

import numpy as np
x1 = np.array([9.6,17.1,64.8,40.9,136.3,182.5,78.3,3.7,26.,0.4])
x2 = np.array([5.,13.,18.1,23.5,29.2,27.5,23.3,16.5,8.2,-0.7])
np.corrcoef(x1, x2)[0, 1]
# 0.7946519586442787

这个计算结果对不对呢？用相同的数据再用excel计算一下。excel中使用CORREL函数计算两列数据的相关系数。可以发现与numpy计算的结果是一致的。因此可以放心的使用np.corrcoef()来计算相关系数。

2 计算栅格的相关系数

2.1 数据准备

以降水和温度的相关系数计算为例，首先准备气温和降水数据，例如10个月的气温和10个月的降水。相关系数计算需要多个数据，只用一期的影像无法计算。

待计算的数据需要有相同的行列数和相同的期数。
通常ndvi、降水、温度的数据是单波段的，首先需要将单波段的影像合并成一个多波段的，例如一个月一个波段或一年一个波段，具体看使用的数据。这个可以使用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等库来操作，此外使用numba加速循环。代码在下面贴出，注释十分的丰富，很利于理解。只需要在最后的main函数中修改输入输出路径即可。

import numpy as np
from osgeo import gdal
import time
from numba import jit

@jit
def compute(arr1,arr2,src_nodta):
    """
    计算相关系数，c为通道数，h为行数，w为列数
    :param arr1: 影像1的数据，np数组，shape为[c,h,w]
    :param arr2: 影像2的数据，np数组，shape为[c,h,w]
    :param src_nodta: 忽略值，数字
    :return: 相关系数图像，np数组，shape为[h,w]
    """
    band1 = arr1[0]
    out = band1 * 0 - 2222
    rows,colmns = out.shape
    for row in range(rows):
        for col in range(colmns):
            if src_nodta is None:
                x1 = arr1[:, row, col]
                x2 = arr2[:, row, col]
                corr = np.corrcoef(x1, x2)[0, 1]
                out[row, col] = corr
            else:
                if band1[row, col] != src_nodta:
                    x1 = arr1[:, row, col]
                    x2 = arr2[:, row, col]
                    corr = np.corrcoef(x1, x2)[0, 1]
                    out[row, col] = corr
    return out

def yiyuanhuigui(imgpath1, imgpath2, outtif):
    """
    计算两个影像的相关系数
    :param imgpath1: 影像1，多波段
    :param imgpath2: 影像2，与影像1的波段数相同、行列数相同
    :param outtif: 输出结果路径
    :return: None
    """
    # 读取影像1的信息和数据
    ds1 = gdal.Open(imgpath1)
    projinfo = ds1.GetProjection()
    geotransform = ds1.GetGeoTransform()
    rows = ds1.RasterYSize
    colmns = ds1.RasterXSize
    data1 = ds1.ReadAsArray()
    print(data1.shape)
    # 读取影像2的数据
    ds2 = gdal.Open(imgpath2)
    data2 = ds2.ReadAsArray()
    src_nodta = ds1.GetRasterBand(1).GetNoDataValue()

    # 创建输出图像
    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    dst_ds = driver.Create(outtif, colmns, rows, 1,gdal.GDT_Float32)
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(projinfo)

    # 删除对象
    ds1 = None
    ds2 = None

    # 开始计算相关系数
    out = compute(data1,data2,src_nodta)
    # 写出图像
    dst_ds.GetRasterBand(1).WriteArray(out)

    # 设置nodata
    dst_ds.GetRasterBand(1).SetNoDataValue(-2222)
    dst_ds = None

if __name__=="__main__":
    tif1 = r"E:\data\一元回归\降水sub.tif"
    tif2 = r"E:\data\一元回归\温度sub.tif"
    outtif = r"E:\data\一元回归\corr3.tif"
    t0 = time.time()
    yiyuanhuigui(tif1, tif2, outtif)
    print(time.time() - t0)

2.3 结果示例

2.4 示例数据获取

本文的示例数据获取方法：公众号【老王GIS】回复关键词【xishu】获取链接。

3 结束语

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