Python遥感实战:土地利用转移矩阵计算实战(附代码)
- 2026-03-23 20:10:17
Python遥感实战:土地利用转移矩阵计算实战(附代码) 本文基于之前公众号文章Python遥感实战：基于随机森林的土地利用土地覆盖分类的基础上，对模块化代码进行优化与更新，本期基于分类后的多时相的土地利用土地覆盖栅格数据进行土地利用转移矩阵的计算。 
尽管这个矩阵看起来很复杂，但是想要理解它很简单。由于我们希望了解土地利用变化情况，那么必然会牵扯到至少两个时间段，那么我们就知道这个矩阵是二维的，那么可以认为它是一个DataFrame 数据框 。接下来我们就需要了解矩阵中的行列含义以及下标的含义。 
这里我们查看一下不同地物对应的栅格值： 
可以看到这个位置在2010年是地物1，但是2020年变成地物2 2. 代码讲解 计算结果输出如下： 
注：上述代码已经全部打包，后期会将集成公众号更多数据分析代码和可视化代码进行整理并进一步更新打包！ 
#遥感处理#遥感土地利用转移矩阵#机器学习#逐像元分析#生态遥感 
点击蓝字，关注我们
本期分享：
土地利用转移矩阵计算实战
完整可复制代码land_trans_matrix.txt
点击关注上方“翔的学术日记”，选择加"星标"置顶重磅干货，第一时间送达
原理讲解
什么是土地利用转移矩阵？
简单来说，它是一张“人口普查表”。 只不过统计的不是人，而是土地类型。它描述了在一段特定的时间内（比如 2010 年到 2020 年），某一区域内各类土地利用类型发生变化的方向和数量。对于土地利用转移矩阵，首先重要的就是矩阵。
在地理信息科学中，土地利用转移矩阵通常用一个方阵S来表示。假设研究区内共有 n 种土地利用类型，那么矩阵 S 可以定义为：
S 代表某种地物的面积或者像元大小；每一行 i 代表 T1（初始时间）的地物类别；每一列 j 代表 T2（结束时间）的地物类别； Sij 就代表表示从 T1 时期的第 i 类土地转移到 T2 时期第 j 类土地的面积。
下面我们思考一下“横向”与“纵向”的数学逻辑。可以得到，行和表示T1 时期第 i 类土地的总面积；列和 表示 T2 时期第 j 类土地的总面积。当 i = j 时，元素 S11, S22, ..., Snn 位于矩阵的对角线上。它们代表了在研究时段内未发生变化的土地部分。
工作流程
1. 示例数据
下面是示例数据，表示2010年和2020年同一个地区的土地利用栅格数据，其中。0代表地物1，1代表地物2..... 5代表地物6，下面我们用ArcMap软件进行重新配色并打开查看：
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.nodatanodata2 = 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_idxrow_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_indicescounts = 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)
求点赞
求分享
求喜欢
点击下方账号名片 关注我们
点“阅读原文”下载文章源代码
本文来自网友投稿或网络内容，如有侵犯您的权益请联系我们删除，联系邮箱：wyl860211@qq.com 。
