【python】批量实现modis数据的辐射定标，大气校正及地形校正

批量实现modis数据的辐射定标，大气校正及地形校正

1.定义辐射定标函数

def radiance_cal(band, gain, bias, scale_factor):
    rad = gain * (band - bias) * scale_factor
    return rad

2.定义反射率计算函数

def reflectance_cal(rad, esun, d):
    refl = np.pi * rad * d ** 2 / (esun * np.cos(np.deg2rad(sza)))
    return refl

3.获取增益、偏移和比例因子等信息

meta = modis_file.GetMetadata('SUBDATASETS')
gain = float(meta['gain'])
bias = float(meta['bias'])
scale_factor = float(meta['scale_factor'])
esun = float(meta['esun'])
d = float(meta['distance'])
sza = float(meta['sza'])

4.可视化输出结果

xsize = modis_file.RasterXSize
ysize = modis_file.RasterYSize
gt = modis_file.GetGeoTransform()
xmin, xres, _, ymax, _, yres = gt
xmax = xmin + xsize * xres
ymin = ymax + ysize * yres
xx, yy = np.meshgrid(np.linspace(xmin, xmax, xsize), np.linspace(ymin, ymax, ysize))
xy = np.column_stack((xx.flatten(), yy.flatten()))
xy = np.flipud(xy)
values = np.flipud(reflectance).ravel()
grid = griddata(xy, values, (xx, yy), method='cubic')
plt.imshow(grid, cmap='gray')
plt.show()

5.批量执行 MODIS 数据的辐射校正、大气校正和地形校正
完整代码

import os
import numpy as np
from osgeo import gdal

# 设置 MODIS 文件路径
data_dir = "/path/to/modis/data/"

# 设置输出文件路径
out_dir = "/path/to/output/dir/"

# 加载大气校正和地形校正模块
import spectral
import spectral.io.envi as envi

# 加载 MODIS 文件列表
file_list = os.listdir(data_dir)

# 循环处理每个 MODIS 文件
for file_name in file_list:
    if file_name.endswith(".hdf"):
        # 设置输出文件名
        out_name = file_name.replace(".hdf", ".tif")

        # 定义函数，用于执行辐射校正
        def radiance_correction(file_name, out_dir, out_name):
            # 加载 MODIS 数据
            hdf_file = gdal.Open(os.path.join(data_dir, file_name), gdal.GA_ReadOnly)
            subdatasets = hdf_file.GetSubDatasets()

            # 选择辐射数据集
            radiance_subdataset = None
            for i in subdatasets:
                if "Radiance" in i[0]:
                    radiance_subdataset = i[0]
                    break
            if radiance_subdataset is None:
                raise Exception("Radiance subdataset not found in HDF file.")

            # 加载辐射数据
            radiance_data = gdal.Open(radiance_subdataset, gdal.GA_ReadOnly)
            radiance_band = radiance_data.GetRasterBand(1)
            radiance_array = radiance_band.ReadAsArray()

            # 加载定标数据
            calibration_subdataset = None
            for i in subdatasets:
                if "Reflectance" in i[0]:
                    calibration_subdataset = i[0]
                    break
            if calibration_subdataset is None:
                raise Exception("Reflectance subdataset not found in HDF file.")
            calibration_data = gdal.Open(calibration_subdataset, gdal.GA_ReadOnly)

            # 加载定标参数
            slope_band = calibration_data.GetRasterBand(1)
            slope_array = slope_band.ReadAsArray()
            intercept_band = calibration_data.GetRasterBand(2)
            intercept_array = intercept_band.ReadAsArray()

            # 执行辐射校正
            radiance_array = radiance_array * slope_array + intercept_array

            # 保存结果
            out_driver = gdal.GetDriverByName("GTiff")
            out_file = out_driver.Create(os.path.join(out_dir, out_name), radiance_band.XSize, radiance_band.YSize, 1, radiance_band.DataType)
            out_file.GetRasterBand(1).WriteArray(radiance_array)

            # 释放资源
            out_file.FlushCache()
            out_file = None
            radiance_band = None
            radiance_data = None
            hdf_file = None
            gdal.UseExceptions()

        # 执行辐射校正
        radiance_correction(file_name, out_dir, out_name)

        # 加载辐射校正结果
        radiance_file = os.path.join(out_dir, out_name)
        radiance_data = gdal.Open(radiance_file, gdal.GA_ReadOnly)
        radiance_band = radiance_data.GetRasterBand(1)
        radiance_array = radiance_band.ReadAsArray