Python遥感实战:土地利用转移矩阵计算实战(附代码)

import numpy as npimport rasterioimport pandas as pddef calculate_luc_matrix(tif_path_t1, path_t2, class_names=None):    """    计算土地利用转移矩阵    """    with rasterio.open(tif_path_t1) as src1, rasterio.open(path_t2) as src2:        # 检查两个 TIF 的形状是否一致        if src1.shape != src2.shape:            raise ValueError("两个输入影像的行列数不一致，请先进行裁剪或重采样！")        data1 = src1.read(1)        data2 = src2.read(1)        # 获取影像定义的 nodata 值        nodata1 = src1.nodata        nodata2 = src2.nodata        # 逻辑：只有当像素值等于影像定义的 nodata 时才剔除        # 如果 0 是有效类，确保 nodata1 != 0 即可        mask = (data1 != nodata1) & (data2 != nodata2)        # 提取有效像元        d1_valid = data1[mask].astype(np.int64)        d2_valid = data2[mask].astype(np.int64)    # 提取所有唯一的类别代号（包括 0）    unique_classes = np.unique(np.concatenate([d1_valid, d2_valid]))    # 建立映射：Class ID -> 矩阵索引 (0, 1, 2...)    class_to_idx = {val: i for i, val in enumerate(unique_classes)}    n = len(unique_classes)    # 初始化空矩阵    matrix = np.zeros((n, n), dtype=np.int64)    # 核心统计逻辑：利用 bincount 快速计算 (比循环快得多)    # 将二维转换关系编码为一维索引： index = row_idx * n + col_idx    row_indices = np.array([class_to_idx[v] for v in d1_valid])    col_indices = np.array([class_to_idx[v] for v in d2_valid])    combined_indices = row_indices * n + col_indices    counts = np.bincount(combined_indices, minlength=n*n)    matrix = counts.reshape((n, n))    # 构建 DataFrame 标签    if class_names is None:        class_names = {}    labels = [class_names.get(c, f"Class_{c}") for c in unique_classes]    df_matrix = pd.DataFrame(        matrix,         index=[f"T1_{l}" for l in labels],         columns=[f"T2_{l}" for l in labels]    )    return df_matrix# --- 调用示例 ---mapping = {0: "地物1", 1: "地物2", 2: "地物3", 3: "地物4", 4:'地物5', 5:'地物6'}df = calculate_luc_matrix("\2010.tif", "\2020.tif", class_names=mapping)