超简单制作栅格数据土地利用转移矩阵(Arcgis和Python实现)

由于大面积矢量数据难以快速进行数学计算，使用栅格数据一样能达到同样的结果。
数据：GLobeland 30
预处理：数据预处理在Arcgis和Python2.7平台上完成，包括批量投影转换、批量处理无效值、批量图幅拼接和掩膜提取。为便于精确计算面积，投影转换选用Albers等积投影，WGS-1984地理坐标系。
数据要求：是适用于面积计算的等积投影坐标系
数据属性如下图（可以不加type,光有代码一样）：

步骤1：打开Spatial Analyst工具-局部分析-合并
选中要处理的不同年份的栅格（注意：30m分辨率的大面积数据跑起来还是较慢，要耐心等待）

结果如下图：
转换结果在属性表可查

2.导出属性表

点击导出，导出txt文件

3.表内分析
首先将txt文件转到表中
注意点逗号

做运算单位为平方千米：栅格分辨率×栅格分辨率/1000000。下拉改变格式

选中三列，插入数据透视表

右下方选项行是起始年份，列是目标年份，值为面积。

结果如下图所示（将类型代码替换为类型即可）

4.结果分析
黄色部分：2000-2020年（类型50：湿地）土地利用类型无变化的面积38932.38。
大红部分：2000-2020年有75556.26面积的（类型30：草地）转化为湿地。
玫红色部分（最下面一行）：2020年各土地利用类型面积。
紫色部分：2000年各土地利用类型面积。

python实现

# -*- coding: utf-8 -*-
import os
import sys
import numpy as np
from sklearn.metrics import confusion_matrix
from osgeo import gdal
# reference: https://blog.csdn.net/sinat_35763722/article/details/117791706
​
def read_tif(path):
    g = gdal.Open(path)
    b = g.GetRasterBand(1)
    data = b.ReadAsArray(0, 0)
    return data
​
# 转移矩阵计算函数
def transition_matrix(Start_LUCC_image, End_LUCC_image):
    # 读取影像矩阵
    Start_LUCC_image = read_tif(Start_LUCC_image)
    End_LUCC_image = read_tif(End_LUCC_image)
    # 重排列数组。flatten是numpy.ndarray.flatten的一个函数，即返回一个一维数组。
    lucc_start_array = np.asarray(Start_LUCC_image).flatten()
    lucc_end_array = np.asarray(End_LUCC_image).flatten()
    # 构建混淆矩阵
    transition_area = confusion_matrix(lucc_start_array, lucc_end_array)
    # 提取有效区域（去除空值），得到土地利用面积转移矩阵
    trans_matr = transition_area[:len(transition_area) - 1, :len(transition_area) - 1]
    print("transition_matrix:", "\n", trans_matr)
    return trans_matr
​
# 马尔科夫预测
def markov(trans_matr, Start_year, End_year, Pred_year):
    sum_Start_year = trans_matr.sum(axis=1)  # 按行求各地类和
    sum_End_year = trans_matr.sum(axis=0)  # 按列求各地类和
    whole_area = trans_matr.sum()  # 计算总栅格数
    P_End_year = sum_End_year / trans_matr.sum()  # 计算初始状态矩阵
    # 计算一阶转移概率矩阵
    Ptrans0 = np.empty(trans_matr.shape)
    for i in range(len(sum_Start_year)):
        Ptrans0[i] = trans_matr[i] / sum_Start_year[i]
    # 计算转移概率矩阵
    n = int((Pred_year - End_year) / (End_year - Start_year))
    print(n)
    E = np.identity(len(sum_Start_year))
    for i in range(n):
        E = np.dot(E, Ptrans0)
    Ptrans = E
    # 计算预测年份状态矩阵
    P_Pred_year = np.dot(P_End_year, Ptrans)
    # 计算预测年份各地类面积
    Pred_year_area = np.array(np.around(P_Pred_year * whole_area), dtype=int)
    np.set_printoptions(formatter={'float': '{: 0.6f}'.format})
    print("转移概率矩阵:", "\n", np.around(Ptrans, decimals=6))
    print("预测年份状态矩阵:", "\n", np.around(P_Pred_year, decimals=6))
    print("预测年份各地类面积:", Pred_year, "年\n", Pred_year_area)
    return Ptrans, P_Pred_year, Pred_year_area
​
# 保存计算结果
def Save(Pred_year, Ptrans, P_Pred_year, Pred_year_area):
    outpath = r"G:\21_能量与三北防护林\2-data\statistic\model_out_semi-humid" + str(Pred_year) + ".txt"
    with open(outpath, 'w') as f:  # 写入numpy.ndarray数据
        f.write('%d' % Pred_year)
        f.write(":\n")
        f.write("转移概率矩阵:\n")
        np.savetxt(f, np.round(Ptrans, 6), delimiter="\t", fmt="%.6f")
        f.write("预测年份状态矩阵:\n")
        np.savetxt(f, np.round(P_Pred_year, 6), delimiter="\t", fmt="%.6f")
        f.write("预测年份各地类面积:\n")
        np.savetxt(f, Pred_year_area, delimiter="\t", fmt="%s")
    f.close()  # 关闭
​
​
if __name__ == '__main__':
​
    lucc_st_img = r'G:\21_能量与三北防护林\2-data\land_use_CIC\TNSP_region\semi-humid\1992.tif' # 开始年份
    lucc_ed_img = r'G:\21_能量与三北防护林\2-data\land_use_CIC\TNSP_region\semi-humid\2015.tif' # 结束年份
    trans_matr = transition_matrix(lucc_st_img, lucc_ed_img)
    print(trans_matr)
    Start_year = 2000
    End_year = 2015
    Pred_year = 2020
    Ptrans, P_Pred_year, Pred_year_area = markov(trans_matr, Start_year, End_year, Pred_year)
    Save(Pred_year, Ptrans, P_Pred_year, Pred_year_area)
​

喜欢的同学一键三连呦！