### 高光谱DN值到反射率的转换 要将高光谱的数字数值（Digital Number, DN）转换为反射率，通常需要经过两个主要步骤：辐射定标和大气校正。以下是具体的实现方法。 #### 辐射定标 辐射定标的目的是将DN值转换为辐亮度（Radiance）。这一步骤依赖于卫星传感器的具体参数，例如增益系数和偏移量。假设已知这些参数，则可以使用以下公式完成转换： \[ L = G \cdot DN + B \] 其中： - \( L \) 是辐亮度， - \( G \) 是增益系数， - \( B \) 是偏移量， - \( DN \) 是原始数字数值。 此过程可以通过 `numpy` 实现如下代码片段： ```python import numpy as np def dn_to_radiance(dn_values, gain, bias): """ Convert Digital Numbers to Radiance. Parameters: dn_values (np.ndarray): Array of digital numbers. gain (float): Gain coefficient from the sensor metadata. bias (float): Bias value from the sensor metadata. Returns: np.ndarray: Array of radiance values. """ return gain * dn_values + bias ``` 通过调用该函数即可获得辐亮度矩阵[^1]。 --- #### 计算反射率 计算反射率需要进一步考虑太阳天顶角、地球与太阳之间的距离以及大气影响等因素。反射率 (\(ρ\)) 的通用公式为： \[ ρ = π \cdot \frac{L}{ESUN} \cdot \cos(\theta_s) \] 其中： - \( ESUN \) 是太阳外辐射通量常数， - \( θ_s \) 是太阳天顶角（单位为弧度）。 具体实现可参考以下代码： ```python from math import pi, cos def calculate_reflectance(radiance, esun, solar_zenith_angle): """ Calculate reflectance from radiance using sun-related parameters. Parameters: radiance (np.ndarray): Array of radiance values obtained earlier. esun (float): Solar exo-atmospheric irradiance constant. solar_zenith_angle (float): Sun zenith angle in degrees. Returns: np.ndarray: Array of reflectance values. """ # Convert solar zenith angle to radians and compute cosine theta_s_rad = solar_zenith_angle / 180.0 * pi cos_theta_s = cos(theta_s_rad) # Compute reflectance reflectance = (pi * radiance) / (esun * cos_theta_s) return reflectance.clip(0, 1) # Ensure valid range [0, 1] ``` 注意，在实际应用中可能还需要调整某些参数以适应特定场景或数据集的要求[^2]。 --- #### 使用GDAL读取和保存数据 为了加载和存储地理空间数据文件，推荐使用 GDAL 库来处理栅格图像。下面是一个完整的流程示例： ```python from osgeo import gdal # 打开输入影像并获取波段数据 dataset = gdal.Open('input_image.tif') band = dataset.GetRasterBand(1) dn_array = band.ReadAsArray() # 假设已知增益(gain=0.0003), 偏移(bias=-0.1)，太阳外辐射通量(esun=1975.6)，太阳天顶角(solar_zenith_angle=30°) gain = 0.0003 bias = -0.1 esun = 1975.6 solar_zenith_angle = 30 # 转换为反射率 radiance = dn_to_radiance(dn_array, gain, bias)[^1] reflectance = calculate_reflectance(radiance, esun, solar_zenith_angle)[^2] # 创建新的GeoTIFF文件用于保存结果 driver = gdal.GetDriverByName('GTiff') output_dataset = driver.Create('output_reflectance.tif', dataset.RasterXSize, dataset.RasterYSize, 1, gdal.GDT_Float32) output_band = output_dataset.GetRasterBand(1) output_band.WriteArray(reflectance) # 设置投影信息和其他元数据 output_dataset.SetProjection(dataset.GetProjection()) output_dataset.SetGeoTransform(dataset.GetGeoTransform()) # 关闭资源 output_dataset.FlushCache() del output_dataset ``` 以上代码展示了如何利用 Python 中的 NumPy 和 GDAL 完成从 DN 到反射率的整个工作流[^3]。 ---