预处理高分四号影像

处理内容：辐射定标、物理量计算、几何校正、波段合成

辐射定标

两种方法：

1.直接辐射定标（参数为1）

2.用公式（双击图像查看辐射增益和偏移）但是只能针对一个波段

注意这里用辐射增益和偏移 得到表观辐亮度

附增益和偏移系数：https://blog.csdn.net/zhangxin4832/article/details/118334914

物理量计算

PMI计算表观反射率

E双击图像查看即lrradiance ，θ为xml文件里Satellitezenith，d参数点击查看

IRS处理计算亮温

几何校正

https://blog.csdn.net/weixin_43626557/article/details/103942202

波段重组

重采样：envi工具里resize data

波段合成

1、用代码

import numpy as np
from osgeo import gdal

def layerstaking(PMS, IRS, merge):
    try:
        # Read IRS data
        IRS_ds = gdal.Open(IRS)
        IRS_geotrans = IRS_ds.GetGeoTransform()
        IRS_cols = IRS_ds.RasterXSize
        IRS_rows = IRS_ds.RasterYSize
        IRS_array = IRS_ds.GetRasterBand(1).ReadAsArray()

        # Read PMS data
        PMS_ds = gdal.Open(PMS)
        PMS_geotrans = PMS_ds.GetGeoTransform()
        PMS_cols = PMS_ds.RasterXSize
        PMS_rows = PMS_ds.RasterYSize
        PMS_lefttop_x, PMS_lefttop_y = PMS_geotrans[0], PMS_geotrans[3]
        PMS_x_res, PMS_y_res = PMS_geotrans[1], PMS_geotrans[5]

        # Calculate the relative position of IRS in PMS
        x_offset = int((IRS_geotrans[0] - PMS_lefttop_x) / PMS_x_res)
        y_offset = int((IRS_geotrans[3] - PMS_lefttop_y) / PMS_y_res)

        # Create a new array for IRS to match the PMS dimensions
        new_IRS_band = np.full((PMS_rows, PMS_cols), -999, dtype=np.float32)

        # Place the IRS data into the corresponding position
        new_IRS_band[y_offset:y_offset+IRS_rows, x_offset:x_offset+IRS_cols] = IRS_array

        # Create a new dataset for output
        driver = gdal.GetDriverByName('GTiff')
        out_ds = driver.Create(merge, PMS_cols, PMS_rows, PMS_ds.RasterCount + 1, gdal.GDT_Float32)

        # Copy PMS bands to the new dataset
        for i in range(PMS_ds.RasterCount):
            out_band = out_ds.GetRasterBand(i + 1)
            out_band.WriteArray(PMS_ds.GetRasterBand(i + 1).ReadAsArray())

        # Add the new IRS band to the output dataset
        out_IRS_band = out_ds.GetRasterBand(PMS_ds.RasterCount + 1)
        out_IRS_band.WriteArray(new_IRS_band)

        # Set geotransform and projection
        out_ds.SetGeoTransform(PMS_ds.GetGeoTransform())
        out_ds.SetProjection(PMS_ds.GetProjection())

        print("Layer stacking completed successfully.")

    except Exception as e:
        print(f"An error occurred: {e}")

    finally:
        # Cleanup
        IRS_ds = None
        PMS_ds = None
        out_ds = None

# Example usage
PMS_path = "\\resample5.tif"
IRS_path = "\\IRS_T_rpcortho.tif"
output_path = "\\merged1.tif"
layerstaking(PMS_path, IRS_path, output_path)

2、build layer stack（需坐标系范围设置）