PCA主成分分析遥感影像融合

Chaoying.

已于 2022-06-06 20:44:09 修改

阅读量2.3k

点赞数 2

分类专栏：遥感与图像处理 GDAL 文章标签： python 算法图像处理遥感 GDAL

于 2022-06-06 20:43:04 首次发布

本文链接：https://blog.csdn.net/TJLCY/article/details/125153878

版权

GDAL 同时被 2 个专栏收录

14 篇文章 2 订阅

订阅专栏

遥感与图像处理

4 篇文章 3 订阅

订阅专栏

融合目的

多光谱影像具有高光谱分辨率和低空间分辨率，全色影像具有高空间分辨率和低光谱分辨率，将这两种影像进行融合就可以得到同时具有高空间分辨率和高光谱分辨率的影像

PCA主成分变换影像融合

PCA主成分分析原理

参考博客：主成分分析（PCA）原理详解

融合步骤

pca

代码

from osgeo import gdal, gdalconst, gdal_array
import numpy as np
from score import *

if __name__ == '__main__':
    pantif = './sub_image/SV1-03_20180407_L2A0001046618_1012102060010001_01-PAN.tif' # 全色影像
    muxtif = './sub_image/SV1-03_20180407_L2A0001046618_1012102060010001_01-MUX.tif' # 多光谱影像
    # 读取全色波段影像
    gdal.AllRegister()
    pands = gdal.Open(pantif, gdalconst.GA_ReadOnly)
    # 宽高
    panX = pands.RasterXSize
    panY = pands.RasterYSize
    # 坐标系
    pansrs = gdal.Dataset.GetSpatialRef(pands)
    # 读取多光谱影像
    muxds = gdal.Open(muxtif, gdalconst.GA_ReadOnly)
    muxsrs = gdal.Dataset.GetSpatialRef(muxds)
    # 数据类型
    muxBand = gdal.Dataset.GetRasterBand(muxds, 1)
    muxType = muxBand.DataType
    # 重采样
    options = gdal.WarpOptions(width=panX, height=panY, srcSRS=muxsrs,
                               dstSRS=pansrs, outputType=gdalconst.GDT_Int32, resampleAlg=gdalconst.GRIORA_Bilinear, creationOptions=['COMPRESS=LZW'])
    gdal.Warp("./pca.tif", muxds, options=options)
    # 主成分分析影像融合
    rsData = gdal_array.LoadFile("./pca.tif")
    nBands, xSzie, ySize = rsData.shape
    rsArray = None
    for i in range(nBands):
        bArray = rsData[i]
        bArray = np.reshape(bArray, (-1, 1))
        if i == 0:
            rsArray = bArray
        else:
            rsArray = np.hstack((rsArray, bArray))
    rsArrayAvg = np.mean(rsArray, axis=0)
    rsArray = rsArray - rsArrayAvg
    u, s, v = np.linalg.svd(rsArray, full_matrices=False)
    rsPCA = rsArray@v.T
    rsPCA1 = rsPCA[:, 0]
    rsPCA1Avg = np.mean(rsPCA1)
    rsPCA1Std = np.std(rsPCA1)
    # 计算全色波段的均值和标准差
    panBand = gdal.Dataset.GetRasterBand(pands, 1)
    panArray = gdal.Band.ReadAsArray(panBand)
    panArray = np.reshape(panArray, (-1, 1))
    panStd = np.std(panArray)
    panAvg = np.mean(panArray)
    # 计算变换系数k1,k2
    k1 = rsPCA1Std/panStd
    k2 = rsPCA1Avg-k1*panAvg
    # 对全色波段进行变换使其与第一主成分的标准差和均值一直
    panArray = panArray*k1+k2
    # 替换第一主成分
    rsPCA[:, 0] = panArray[:, 0]
    # 主成分逆变换
    rsArray = rsPCA@v+rsArrayAvg
    rsds = gdal.Open("./pca.tif", gdalconst.GA_Update)
    for i in range(nBands):
        bArray = rsArray[:, i]
        bArray = np.reshape(bArray, (ySize, xSzie))
        iBand = gdal.Dataset.GetRasterBand(rsds, i+1)
        gdal.Band.WriteArray(iBand, bArray)
    del rsds, pands, muxds
    print('ok!')