Python遥感图像处理应用篇(十七)：GDAL 将归一化处理csv数据转化为多波段遥感影像

<think>嗯，用户现在让我用Python代码实现第二问中的光谱角制图（SAM）算法，主要使用gdal和numpy库。首先，我得回顾一下之前给出的算法流程，确保代码结构和步骤正确。 首先，用户的需求是实现矿物填图，输入是标准矿物光谱和卫星高光谱反射率数据。核心步骤包括数据预处理、光谱角计算、阈值分类和结果输出。我需要将这些步骤转化为Python代码。 接下来，我需要考虑如何用gdal读取和写入遥感影像数据。GDAL库处理地理空间数据很常用，可以读取多波段的高光谱数据。对于标准光谱，可能是一个包含各矿物光谱的文本文件或CSV，需要加载进内存作为参考光谱。 数据预处理部分需要归一化处理。每个像元的光谱向量和标准光谱向量都要进行归一化，即除以它们的模长。这可以通过numpy的向量运算来实现，比如计算每个像元的L2范数，然后逐元素相除。 然后，计算光谱角。这里需要计算两个归一化后的向量的点积，然后取反余弦得到角度。这一步需要高效处理，因为高光谱数据通常很大，逐像元处理可能会很慢。可能需要利用numpy的广播机制和矩阵运算来优化，避免使用循环，提升计算速度。 阈值分类部分，需要为每个像元找到最小的光谱角，并与设定的阈值比较。如果小于阈值，则分配对应的矿物类别，否则标记为未知。这里可能需要使用numpy的argmin和where函数来实现。 输出结果方面，需要将分类结果保存为GeoTIFF文件，保留原始的地理参照信息。GDAL的驱动和创建数据集的功能可以完成这一步，需要注意设置正确的投影和地理变换参数。 现在要考虑代码的结构。首先读取输入的高光谱影像，获取其元数据如行数、列数、波段数等。然后加载标准光谱库，可能是从CSV或其他格式读取，确保波段与影像数据对应。归一化处理时，可能需要对每个像元进行向量归一化，这可以通过numpy的einsum或者reshape操作来实现高效计算。 计算光谱角时，可能需要将标准光谱库转换为一个二维数组，每一行代表一个矿物的光谱。然后通过矩阵乘法计算每个像元与所有标准光谱的点积，再取反余弦得到角度矩阵。这一步是关键，需要确保矩阵运算的正确性，避免维度错误。 接下来，找到每个像元的最小角度和对应的矿物索引，应用阈值分类。最后，将结果写入输出文件，并设置与输入文件相同的地理信息。 还需要处理可能的异常情况，比如输入数据的波段顺序是否与标准光谱一致，数据是否有无效值（如NaN），这些可能需要预处理步骤，比如填充或忽略无效值。 另外，测试代码的效率和内存使用也很重要。高光谱数据通常很大，如果一次性读入内存可能会导致内存不足。可能需要分块处理数据，逐块读取、处理和写入，这样可以减少内存压力。但考虑到用户希望尽量简洁，可能暂时不处理分块，假设数据可以全部载入内存。 最后，确保代码的注释和结构清晰，方便用户理解和修改参数，如输入输出路径、阈值设置等。</think>以下是基于光谱角制图（SAM）的矿物填图Python代码实现，使用`gdal`和`numpy`库： ```python import numpy as np from osgeo import gdal def sam_classification(image_path, spectra_lib, output_path, threshold_deg=5): """ 使用光谱角制图（SAM）对高光谱影像进行矿物分类 :param image_path: 输入高光谱影像路径（已大气校正） :param spectra_lib: 标准光谱库字典，格式为{矿物名称: 光谱数组} :param output_path: 输出分类结果路径 :param threshold_deg: 角度阈值（度） """ # 读取高光谱影像 ds = gdal.Open(image_path) rows = ds.RasterYSize cols = ds.RasterXSize bands = ds.RasterCount geotrans = ds.GetGeoTransform() proj = ds.GetProjection() # 将影像数据读取为numpy数组（形状为[行, 列, 波段]） img = ds.ReadAsArray().transpose(1, 2, 0) # 转换为[rows, cols, bands] ds = None # 关闭数据集 # 准备标准光谱库 min_names = list(spectra_lib.keys()) ref_spectra = np.array([spectra_lib[name] for name in min_names]) # 数据预处理：归一化 # 计算每个像元的L2范数 img_norm = np.linalg.norm(img, axis=2, keepdims=True) img_norm[img_norm == 0] = 1e-10 # 避免除以零 img_normalized = img / img_norm # 归一化参考光谱 ref_norm = np.linalg.norm(ref_spectra, axis=1, keepdims=True) ref_normalized = ref_spectra / ref_norm # 计算光谱角（使用矩阵运算优化速度） # 计算点积：[rows, cols, num_minerals] dot_product = np.dot(img_normalized, ref_normalized.T) # 计算角度（弧度转角度） sam_angles = np.arccos(np.clip(dot_product, -1, 1)) * 180.0 / np.pi # 分类：找到最小角度及对应类别 min_angles = np.min(sam_angles, axis=2) min_indices = np.argmin(sam_angles, axis=2) # 应用阈值 classified = np.where(min_angles < threshold_deg, min_indices + 1, 0) # 0表示未知 # 写入输出结果 driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(output_path, cols, rows, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(geotrans) out_ds.SetProjection(proj) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(classified) out_band.SetNoDataValue(0) out_ds = None # 示例用法 if __name__ == "__main__": # 假设标准光谱库（示例数据，需替换为实际值） spectra_lib = { "Sdol": np.array([0.1, 0.2, 0.3, ...]), # 白云石光谱（波段需与影像对齐） "Scal": np.array([0.15, 0.25, 0.35, ...]), # 方解石 # 添加其他矿物光谱... } sam_classification( image_path="input_hyperspectral.tif", spectra_lib=spectra_lib, output_path="mineral_classification.tif", threshold_deg=5 ) ``` --- ### 代码说明 1. **数据读取** - 使用GDAL读取高光谱影像，获取空间参考信息（地理变换和投影） - 将数据转换为`[行, 列, 波段]`的numpy数组格式 2. **光谱归一化** - 对每个像元光谱和标准光谱进行L2归一化： $$ \hat{R} = \frac{R}{\|R\|_2}, \quad \hat{S}_i = \frac{S_i}{\|S_i\|_2} $$ 3. **光谱角计算** - 利用矩阵运算加速点积计算： $$ \cos\theta = \hat{R} \cdot \hat{S}_i^T $$ - 将弧度转换为角度： $$ \theta = \arccos(\cos\theta) \times \frac{180}{\pi} $$ 4. **分类与输出** - 找到每个像元的最小光谱角及其对应的矿物索引 - 根据阈值生成分类结果（0表示未知类别） - 输出为GeoTIFF文件，保留原始地理信息 --- ### 关键优化点 - **矩阵运算加速**：避免逐像元循环，利用`numpy`的广播机制 - **内存管理**：及时释放不再使用的数据（如`ds = None`） - **数值稳定性**：添加`np.clip`防止浮点误差导致`arccos`计算错误 实际使用时需确保： 1. 标准光谱与影像波段范围、数量和顺序一致 2. 输入影像已完成辐射定标和大气校正 3. 根据具体任务调整阈值（通常通过试验确定）