基于Python构建土地利用转移矩阵及完成Markov预测

助教课遇到这个题目，用python写了一个，放在这里备忘吧。

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Markov.py
# Created on: 2021-04-03 21:32:05
#   (generated by WCJ)
# Description: 
# ---------------------------------------------------------------------------

import os
import sys
import numpy as np
from sklearn.metrics import confusion_matrix
from PIL import Image
# sys.path.append(os.pardir)

def transition_matrix(Start_LUCC_image, End_LUCC_image):
    # 读取影像灰度值（包含空值）
    Start_LUCC_image = Image.open(Start_LUCC_image)
    End_LUCC_image = Image.open(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 = "F:\\test\FilesGenerate\\model_out_" + str(Pred_year) + ".txt"
    with open(outpath, 'w') as f:  # 写入numpy.ndarray数据
        # with open("F:\\test\FilesGenerate\\model_out.txt", '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")
        # 使用numpy.savetxt()写入数据，Data为要存的变量，因为numpy.ndarray数据无法用write()写入,数据间用'，'相隔。
        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.write("\n")  # 换行
    f.close()  # 关闭


if __name__ == '__main__':
    # 读取影像灰度值（包含空值）
    lucc_2010_img = r'F:\test\FilesGenerate\Re_guangzhou2010.tif'
    lucc_2015_img = r'F:\test\FilesGenerate\Re_guangzhou2015.tif'
    trans_matr = transition_matrix(lucc_2010_img, lucc_2015_img)
    Ptrans, P_Pred_year, Pred_year_area = markov(trans_matr, Start_year=2010, End_year=2015, Pred_year=2020)
    Save(2020, Ptrans, P_Pred_year, Pred_year_area)

单纯做Markov预测的话，可以利用一个师兄开发的CA软件中的相关模块来做：
http://www.geosimulation.cn/FLUS.html