<think>我们已知用户有两幅TIFF格式的MCD12Q1土地覆盖数据，需要计算转换矩阵（即土地覆盖类型转移矩阵）。 用户之前遇到了文件路径错误，现在指定了文件位置：E:\bishe\shuju\wutuditiff 但用户没有提供具体的文件名，因此我们需要假设文件名或让用户输入。 然而，根据用户之前的错误，可能是2019和2020年的数据，但文件名可能是“2019.LC_Type1.tif”和“2020.LC_Type1.tif”之类的。 为了确保代码的健壮性，我们将： 1. 要求用户输入两个文件的完整路径，或者 2. 在指定目录下搜索符合条件的文件（但可能有多个文件，需要用户指定顺序或按时间排序） 考虑到用户可能希望直接运行，我们可以让用户输入两个文件名（或年份）作为输入。 但问题描述中要求计算两幅图像，因此我们设计两种方案： 方案A：用户通过输入指定两个文件路径（或相对路径） 方案B：用户输入两个年份，我们根据年份构建文件名（因为MCD12Q1数据通常按年份存储） 由于用户没有提供具体文件名，我们假设文件名模式为“年份.LC_Type1.tif”，例如“2019.LC_Type1.tif” 我们将采用方案B，让用户输入两个年份，然后构建文件路径。 另外，注意处理可能的异常（如文件不存在、读取错误等）以及确保两幅图像的地理信息一致。 步骤： 1. 获取用户输入的两个年份（如2019和2020） 2. 构建文件路径：base_dir = r"E:\bishe\shuju\wutuditiff" file1 = os.path.join(base_dir, f"{year1}.LC_Type1.tif") file2 = 同理 3. 读取两个TIFF文件，并确保它们具有相同的尺寸和投影（如果不一致则需要重采样或对齐，但用户已预处理，我们假设一致） 4. 计算转移矩阵（一个n×n的矩阵，n=17，因为MCD12Q1的IGBP分类有17类） 注意：MCD12Q1的无效值通常是255或0（根据处理方式），我们需要在计算时排除无效值。 我们将修改之前的代码，加入无效值处理。 代码实现： </think> ### 计算MODIS土地覆盖转换矩阵的完整解决方案 以下是针对您指定路径（`E:\bishe\shuju\wutuditiff`）计算两幅MCD12Q1图像转换矩阵的完整代码，包含自动文件检测和异常处理： ```python import os import rasterio import numpy as np import glob # 1. 自动检测目录下的TIFF文件 data_dir = r"E:\bishe\shuju\wutuditiff" tiff_files = glob.glob(os.path.join(data_dir, "*.tif")) + glob.glob(os.path.join(data_dir, "*.tiff")) if len(tiff_files) < 2: raise FileNotFoundError(f"目录中需要至少2个TIFF文件，当前找到: {len(tiff_files)}") print("找到的TIFF文件:") for i, f in enumerate(tiff_files): print(f"{i+1}. {os.path.basename(f)}") # 2. 用户选择文件（或自动选择最新两个） choice = input("输入选择的文件编号（如'1 2'），或回车自动选择最新两个: ") if choice.strip(): idx1, idx2 = map(int, choice.split()) file1 = tiff_files[idx1-1] file2 = tiff_files[idx2-1] else: # 按修改时间排序，选最新两个 sorted_files = sorted(tiff_files, key=os.path.getmtime, reverse=True) file1, file2 = sorted_files[:2] print(f"自动选择: {os.path.basename(file1)} 和 {os.path.basename(file2)}") # 3. 安全读取图像数据 def read_lc_tiff(file_path): """读取土地覆盖TIFF并处理无效值""" with rasterio.open(file_path) as src: data = src.read(1) # MCD12Q1无效值处理 (0-16为有效分类) data[(data < 0) | (data > 16)] = 255 return data, src.profile try: img1, meta1 = read_lc_tiff(file1) img2, meta2 = read_lc_tiff(file2) # 4. 验证图像一致性 if img1.shape != img2.shape: raise ValueError(f"图像尺寸不匹配: {img1.shape} vs {img2.shape}") if meta1['crs'] != meta2['crs']: print("警告: 坐标系不一致! 结果可能不准确") # 5. 计算转换矩阵 n_classes = 17 # IGBP分类体系 transition_matrix = np.zeros((n_classes, n_classes), dtype=np.int64) # 使用向量化操作提高性能 valid_mask = (img1 != 255) & (img2 != 255) flat1 = img1[valid_mask].flatten() flat2 = img2[valid_mask].flatten() # 使用bincount快速统计 for i in range(n_classes): mask_i = (flat1 == i) counts = np.bincount(flat2[mask_i], minlength=n_classes) transition_matrix[i] = counts # 6. 输出结果 print("\n土地覆盖转换矩阵 (行: 起始类型, 列: 目标类型):") print(transition_matrix) # 7. 保存结果 output_path = os.path.join(data_dir, "transition_matrix.csv") np.savetxt(output_path, transition_matrix, delimiter=",", fmt="%d") print(f"\n结果已保存至: {output_path}") except Exception as e: print(f"处理失败: {str(e)}") # 详细错误诊断 print("\n诊断信息:") print(f"文件1: {file1} | 尺寸: {img1.shape if 'img1' in locals() else '未读取'}") print(f"文件2: {file2} | 尺寸: {img2.shape if 'img2' in locals() else '未读取'}") ``` ### 使用说明： 1. **运行流程**： - 自动扫描目录检测所有TIFF文件 - 用户选择或自动配对两个文件 - 自动处理无效值（255）和异常数据 - 输出矩阵并保存CSV文件 2. **特殊处理**： - **内存优化**：使用向量化操作替代循环，处理大型图像 - **无效值过滤**：自动排除255值（MCD12Q1标准无效值） - **坐标系检查**：警告坐标系不一致的情况 3. **输出示例**： ``` 土地覆盖转换矩阵: [[ 0 0 0 ... 0] [ 0 120 15 ... 0] ... [ 0 0 0 ... 42]] ``` ### 常见问题解决方案： 1. **文件找不到**： ```python # 手动指定文件名 file1 = r"E:\bishe\shuju\wutuditiff\2019_LC.tif" file2 = r"E:\bishe\shuju\wutuditiff\2020_LC.tif" ``` 2. **内存不足**： ```python # 分块处理大文件 block_size = 1000 # 每次处理1000x1000像素块 for i in range(0, img1.shape[0], block_size): for j in range(0, img1.shape[1], block_size): block1 = img1[i:i+block_size, j:j+block_size] block2 = img2[i:i+block_size, j:j+block_size] # 对每个块计算部分矩阵... ```