• 日期: 2028-05-25
  • 标签: 常规

Python 计算多阶邻接度指数 (MAI) - 城市扩张分析

该脚本使用 Python 计算多阶邻接度指数 (MAI) 来分析城市扩张。它利用 GDAL 库读取栅格数据，并通过创建缓冲区来确定原始城市区域与新建城市区域之间的关系。

代码

import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt
from tqdm import trange

# 定义函数：将栅格数据转换为二进制矩阵
def read_img_data(filename):
    dataset = gdal.Open(filename)  # 打开文件
    im_width = dataset.RasterXSize  # 栅格矩阵列数
    im_height = dataset.RasterYSize  # 栅格矩阵行数
    im_data = dataset.ReadAsArray(0, 0, im_width,
                                  im_height)  # 获取栅格数据矩阵
    del dataset
    return im_data

# 定义函数：获取栅格数据投影信息
def read_img_proj(filename):
    dataset = gdal.Open(filename)  # 打开文件
    im_proj = dataset.GetProjection()  # 获取地图投影信息
    return im_proj

# 定义函数：获取栅格数据仿射矩阵信息
def read_img_geotrans(filename):
    dataset = gdal.Open(filename)  # 打开文件
    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵，左上角像素的地理坐标和像素分辨率
    return im_geotrans

# 定义函数：写入 GeoTIFF 文件
def write_img(filename, im_proj, im_geotrans, im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影信息

    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

# 定义函数：计算多阶邻接度指数
def MAI(Original_urban_construction_land_file, New_urban_construction_land_file):
    urban_construction_land = 1  # 定义栅格中城市建设用地的值为 1

    # 读取原始城市建设用地栅格数据
    Original_urban_construction_land_data = read_img_data(Original_urban_construction_land_file)

    # 读取新建城市建设用地栅格数据
    New_urban_consturction_land_data = read_img_data(New_urban_construction_land_file)

    # 创建一个和新建城市用地栅格数据有相同行列数的空数组，用于存储计算后的多阶邻接度指数值
    MAI_array = np.full((New_urban_consturction_land_data.shape[0],
                         New_urban_consturction_land_data.shape[1]),
                        fill_value=-9999,
                        dtype=np.float64)

    # 遍历新建城市建设用地栅格数据的二维矩阵
    for height in range(New_urban_consturction_land_data.shape[0]):
        for width in range(New_urban_consturction_land_data.shape[1]):
            # 如果检测到新建城市建设用地，为其建立多级缓冲区
            if New_urban_consturction_land_data[height, width] == urban_construction_land:
                # 初始缓冲区为 3*3 的邻域窗口
                buffer_init = 3
                buffer_init_k = 1
                buffer_flag = 0

                # 当缓冲区内未检测到原始城市用地时，建立多级缓冲区，直到检测到为止
                while buffer_flag == 0:
                    # 创建一个空列表，以便后面可以计算缓冲区内原始城市用地的个数
                    buffer_urban = []

                    # 遍历缓冲区内的栅格，查找原始城市建设用地的存在情况
                    for height_height in range(buffer_init):
                        for width_width in range(buffer_init):
                            # 只检测 n*n（n>=3）的缓冲区最外层栅格
                            if (height_height == 0
                                    or height_height == buffer_init - 1) or (
                                        (height_height != 0 and height_height != buffer_init - 1)
                                    and (width_width == 0 or width_width == buffer_init - 1)):
                                # 将缓冲区内的栅格值存入之前创建的列表空间中
                                buffer_urban.append(
                                    Original_urban_construction_land_data[
                                        height - buffer_init_k + height_height,
                                        width - buffer_init_k + width_width])

                    # 如果在缓冲区内找到原始城市用地，则计算多阶邻接度指数
                    if urban_construction_land in buffer_urban:
                        Ai = buffer_urban.count(urban_construction_land)  # 计算缓冲区内原始城市用地的面积
                        A0 = len(buffer_urban)  # 计算缓冲区面积
                        N = buffer_init_k  # 计算缓冲区级数
                        MAI_cal = N - (Ai / A0)  # 根据 MAI 公式计算新建城市建设用地的 MAI 值
                        MAI_array[height, width] = MAI_cal  # 将计算的 MAI 值放入相应的空数组中
                        break  # 停止缓冲，遍历下一个栅格
                    # 如果在缓冲区内未找到原始城市用地，则 n+1，扩大缓冲区
                    if urban_construction_land not in buffer_urban:
                        buffer_init += 2
                        buffer_init_k += 1

    # 显示计算后的 MAI 栅格数据
    plt.imshow(MAI_array)
    plt.colorbar()
    plt.show()
    return MAI_array

# 原始城市建设用地栅格数据路径
Original_urban_construction_land_file_path = 'E:\MAI\test\mai\old'
# 新建城市建设用地栅格数据路径
New_urban_construction_land_file_path = 'E:\MAI\test\mai\new'

# 新建城市建设用地栅格数据文件
New_urban_construction_land_file = [
    '2005-2010.tif', '2010-2015.tif', '2015-2020.tif'
]

# 原始城市建设用地栅格数据文件
Original_urban_construction_land_file = ['2005.tif', '2010.tif', '2015.tif']

# 遍历所有数据，计算每个时期的 MAI 值
for urban_year in trange(len(New_urban_construction_land_file)):
    MAI_calcu = MAI(
        Original_urban_construction_land_file_path + '/' +
        Original_urban_construction_land_file[urban_year],
        New_urban_construction_land_file_path + '/' +
        New_urban_construction_land_file[urban_year])
    
    # 将计算出的 MAI 值写入到 GeoTIFF 文件中
    write_img(
        str(New_urban_construction_land_file[urban_year][:-4] + '_MAI.tif'),
        read_img_proj(New_urban_construction_land_file_path + '/' +
                      New_urban_construction_land_file[urban_year]),
        read_img_geotrans(New_urban_construction_land_file_path + '/' +
                          New_urban_construction_land_file[urban_year]),
        MAI_calcu)

解释

  1. 导入库: 导入必要的库，包括 NumPy、GDAL、Matplotlib 和 tqdm。
  2. 函数定义:
    • read_img_data: 读取栅格文件并将其转换为 NumPy 数组。
    • read_img_proj: 获取栅格文件的投影信息。
    • read_img_geotrans: 获取栅格文件的仿射变换矩阵。
    • write_img: 将 NumPy 数组写入 GeoTIFF 文件。
    • MAI: 计算多阶邻接度指数。
  3. 数据路径设置: 设置原始城市用地和新建城市用地的栅格文件路径。
  4. 数据文件列表: 定义原始城市用地和新建城市用地栅格文件的名称列表。
  5. 循环计算 MAI: 使用 trange 循环遍历每个时间段的城市用地数据，并使用 MAI 函数计算 MAI 指数。
  6. 输出结果: 将计算出的 MAI 值写入 GeoTIFF 文件。

使用说明

  1. 安装必要的库：pip install numpy gdal matplotlib tqdm
  2. 将代码中的数据路径更改为实际的路径。
  3. 运行脚本。

输出

脚本将为每个时间段的城市用地数据创建一个新的 GeoTIFF 文件，其中包含计算出的 MAI 值。

参考资料

附加信息

  • 该代码示例仅供参考，您可以根据您的具体需求进行调整。
  • MAI 指数可用于分析城市扩张模式，例如城市蔓延程度和城市核心区的稳定性。
  • 您可以根据不同的研究目的选择不同的缓冲区大小和计算方法。
  • 您可以使用其他 GIS 软件或库来可视化和分析计算结果。
Python 计算多阶邻接度指数 (MAI) - 城市扩张分析

原文地址: https://www.cveoy.top/t/topic/mVyN 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录