使用ArcGIS(ArcMap)自带的Arcpy包进行简单的python脚本自动化——以计算武汉部分地区面积气温降水等参数为例

kukeer2

已于 2024-10-09 14:18:22 修改

阅读量489

点赞数 21

文章标签： python arcgis 经验分享自动化

于 2024-10-09 14:18:21 首次发布

本文链接：https://blog.csdn.net/kukeer2/article/details/142749140

版权

因为这是第一次写csdn的文章，所以有什么文章上逻辑错误或者改进都可以私信发给我，希望大家多多评论。

本项目旨在利用 ArcGIS (ArcMap) 中的 Arcpy 包，对武汉市5种地形的气温、降水量、NPP（净初级生产力）等环境参数进行自动化计算和分析。项目的主要步骤包括从武汉的Shapefile 文件及相应的 PNG 图像中提取和计算地形的面积、气象数据及碳汇量，最后将计算结果输出到 Excel 文件中。通过自动化 Python 脚本，显著提高了数据处理的效率，减少了手动操作的复杂性。整个流程充分展示了如何将地理信息系统 (GIS) 与编程技术结合，从而实现精确的空间数据分析和环境监测

1 通过ArcMap计算操作简述（简单介绍，具体步骤不详细描述）

1.6 使用 IDW 工具进行插值，修改坐标系，再把栅格转化为点，

1.7 进行空间连接，将点数据连接到多边形数据，并计算平均值

2 python自动化

2.1 自动配准

2.2 ArcMap过程自动化（raster_processing）

2.3自动化选择地形进行计算

2.4脚本文件

首先先介绍一下任务。

我们有武汉市_省界.shp ，气温降水npp的png图和各种地形的多边形栅格图。要利用这些图去计算每个地形所对应的面积、气温、降水量、NPP量和碳汇量，再计算总碳汇量。

总的流程图如下，里面有相关的公式和具体计算流程。我就不一一解释了，有问题可以私信我。主要是需要将这个过程用ArcMap完成。所以我们先用ArcMap进行一次操作。再给大家介绍通过ArcGIS自带的Arcpy包使用python编写自动化程序，这样会节省手动的很多操作。

1 通过ArcMap计算操作简述（简单介绍，具体步骤不详细描述）

1.1 给PNG图加坐标信息

首先在ArcMap界面中导入我们的武汉市_省界.shp

再导入气温png图（其他两个图步骤一样，这里以气温为例）。并且进行手动的地理配准。

接着导出格式为tiff的栅格数据。这样我们就给降雨图加入了坐标信息。

1.2 进行重分类

1.3 将栅格转化为多边形

1.4 修改地理坐标系

1.5 将栅格转化为点

1.6 使用 IDW 工具进行插值，修改坐标系，再把栅格转化为点，

1.7 进行空间连接，将点数据连接到多边形数据，并计算平均值

这样就得到了每个地区的面积，并且得到每块栅格的温度值（华氏度，需要转化成摄氏度）。降雨和NPP同样，最后将得到的数值放入前面流程图里的公式里就能计算出每个类别地形的相应数值。但是弊端也十分明显，就是过程复杂，并且每一种图都得单独操作一遍，计算也是需要一步一步进行手动计算。所以我采用ArcGIS自带的python包Arcpy将这一过程写成用python执行的自动化过程。下面我将介绍这个过程。

在介绍这个过程之前，我得先说明一个问题就是我使用的是ArcGIS10.8所自带的Arcpy，这个包支持的python版本是2.7，已经是很古老的版本，很多语法和现在常用的python3版本有不一样的地方，所以有些语句和现在常用的python语法有些不一致，ArcGISpro好像是支持python3版本，感兴趣的朋友可以去试一下。环境配置我就不描述了，有问题的朋友可以私信我。

2 python自动化

2.1 自动配准

因为气温，降水和npp是png格式图片，png格式图片没有地理坐标系信息，所以需要将png格式的图像配准到参考tif图像的地理坐标系中。这里我们使用GDAL库，GDAL是一个用于读取和写入地理空间栅格数据格式的开源库。它提供了一套工具和API用于处理各种类型的地理空间数据。因为Arcpy库里没有自动配准的相关代码，所以这步与后面的步骤需要分开进行。

def read_image(path):
    """读取影像"""

    dataset = gdal.Open(path)
    # 栅格矩阵的列数
    im_width = dataset.RasterXSize
    # 栅格矩阵的行数
    im_height = dataset.RasterYSize
    # 波段数
    im_bands = dataset.RasterCount
    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
    im_proj = dataset.GetProjection()  # 获取投影信息
    arr = dataset.ReadAsArray()
    del dataset

    return im_width, im_height, im_geotrans, im_proj, arr

首先我们要读取和png配准的tif格式图像的信息，使用gdal.Open(path)打开指定路径的影像文件，这里path是影像文件的路径

dataset变量保存了这个影像文件的GDAL数据集对象。

im_width = dataset.RasterXSize: 获取栅格矩阵的列数，即图像的宽度。

im_height = dataset.RasterYSize: 获取栅格矩阵的行数，即图像的高度。

im_bands = dataset.RasterCount: 获取影像的波段数，即图像包含多少个波段（例如，RGB图像有三个波段）。

im_geotrans =dataset.GetGeoTransform(): 获取影像的仿射变换参数，通常用于描述图像空间坐标和地理坐标之间的转换关系。

im_proj = dataset.GetProjection(): 获取影像的投影信息，描述影像在地理空间中的投影格式

arr = dataset.ReadAsArray(): 读取影像数据信息并以数组形式返回。数组的形状通常是(波段数, 高度, 宽度)。

这个函数最后返回一个包含图像宽度、高度、仿射变换、投影信息和波段数据的元组。

def writeTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):
    """将数组写成图像"""

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2:
        im_data = np.array([im_data])
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape
        # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, im_width, im_height, im_bands, gdal.GDT_Byte)
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

这个函数writeTiff的目的是将一个数组保存为TIFF格式的图像文件。该函数首先检查im_data的形状(维度)，并根据不同的形状将数据调整为统一的格式，以便在后续处理中正确写入文件。例如，如果数据是二维的，则将其转换为单一波段的三维数组。gdal.GetDriverByName("GTiff")获取GeoTIFF文件的驱动，用driver.Create()创建一个新文件。如果数据集成功创建，为其设置仿射变换和投影信息。SetGeoTransform方法用于指定像素到地理坐标的转换参数，而SetProjection方法则为文件指定地理坐标系。通过遍历波段数，使用dataset.GetRasterBand(i + 1)选择具体波段，然后使用WriteArray将相应波段的数据写入文件。

if __name__ == '__main__':
    print('正在进行影像读取...')
    tif_path = 'Registration/Manual_registration/png.tif' #参考图像
    pngs_path = 'Registration/PNG'
    out_path = 'Run_time_folder/tif'

    dir = os.path.exists(out_path)
    if dir != True:
        os.makedirs(out_path)

    pngs = [i for i in os.listdir(pngs_path) if i.split('.')[-1] == 'png']
    for i in pngs:
        print(i)
        im_width, im_height, im_geotrans, im_proj, arr = read_image(tif_path)
        png_path = pngs_path + '\\' + i  # 配准图像
        input_dataset = gdal.Open(png_path)
        if input_dataset is not None:
            # 获取输入数据集的驱动程序
            input_driver = input_dataset.GetDriver()
            # 指定新文件名
            output_file = out_path + '\\' + i
            # 使用 CreateCopy 方法创建新的数据集副本
            output_dataset = input_driver.CreateCopy(output_file, input_dataset)
            if output_dataset is not None:
                output_dataset.SetGeoTransform(im_geotrans)
                output_dataset.SetProjection(im_proj)
                # 关闭数据集
                input_dataset = None
                output_dataset = None
                print("created successfully as {}.".format(output_file))
            else:
                print("Error: Unable to create copy.")
        else:
            print("Error: Unable to open {}.".format(png_path))

然后就是主函数，利用前面的writeTiff和read_image函数读取参考TIFF文件的地理信息，将存储在指定目录下的PNG图像文件转换为具有相同地理参照的TIFF文件。并且用os包在没有目录的时候创建该目录，确保后续生成的TIFF文件能正确保存。但是有个问题这段代码是用python3写的。。。所以可能需要大家换个配置环境,或者嫌麻烦的话可以手动配准，哈哈哈哈哈。

2.2 ArcMap过程自动化（raster_processing）

这里我们使用Arcpy的相关函数，感兴趣的同学可以去官网查看函数的使用方法了解工具语法—ArcMap | 文档 (arcgis.com)

import arcpy
import numpy as np
import os

arcpy.env.workspace = workspace
# 新建文件夹
tiff_filename = os.path.basename(input_tiff)
filename = os.path.basename(input_raster)
new_folder = "process_" + tiff_filename + "_" + filename
folder_path = os.path.join(workspace, new_folder)
if not os.path.exists(folder_path):
    os.makedirs(folder_path)


# 输出重分类后的栅格路径
output_reclass = os.path.join(workspace, new_folder, "reclass1.tif")
# 输出栅格转换为多边形的路径
output_polygon = os.path.join(workspace, new_folder, "polygon1.shp")
# 输出monthly_temperature栅格转换为点的路径
output_point_feature = os.path.join(workspace, new_folder, "monthly_temperature_point.shp")
# 输出IDW进行插值的栅格路径
output_raster_idw = os.path.join(workspace, new_folder, "Idw_shp1")
# 输出out_idw栅格转换为点文件路径
output_point_feature_idw = os.path.join(workspace, new_folder, "idw_point.shp")
# 输出连接结果路径
output_join_result = os.path.join(workspace, new_folder, "Join_Output.shp")

首先我们需要定义工作空间，并且在指定的工作空间下创建一个新文件夹，并且存储不同步骤的结果，这样可以保证不同类别生成的tif在不同的文件夹里，避免发生错误。

# 重分类
reclass_values = arcpy.sa.RemapValue([[0, "NODATA"], [1, 254, 1], [255, "NODATA"]])
out_reclass = arcpy.sa.Reclassify(input_tiff, "Value", reclass_values)
out_reclass.save(output_reclass)
 # print("重分类完成，输出路径为: " + output_reclass)

这里使用 arcpy.sa.RemapValue 创建一个重新映射对象 reclass_values，将像素值0和255转换为无数据（NODATA），并将1到254的像素值重新分类为1。再用 arcpy.sa.Reclassify 函数根据定义的 reclass_values 对 input_tiff 栅格文件进行重分类。这里指定字段为"Value"，这意味着重分类按照像素的值进行。

# 栅格转换为多边形
arcpy.RasterToPolygon_conversion(output_reclass, output_polygon, "NO_SIMPLIFY")
# print("栅格转换为多边形完成，输出路径为: " + output_polygon)

# monthly_temperature栅格转换为点
arcpy.conversion.RasterToPoint(input_raster, output_point_feature)
# print("monthly_temperature栅格转换为点，转换完成。")

arcpy.RasterToPolygon_conversion函数将重分类后的栅格文件（output_reclass）转换为一个矢量多边形图层，结果保存到output_polygon指定的路径。参数 "NO_SIMPLIFY" 表示在转换过程中不简化输出的多边形边界。

arcpy.conversion.RasterToPoint函数将input_raster栅格数据转换成点矢量文件并保存到output_point_feature路径中。每个栅格像素中心会被转换为一个点，该点保留栅格像素的值。

# 使用 IDW 工具进行插值
cell_size = 100
original_coord_system = arcpy.env.outputCoordinateSystem
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4548)
out_idw = arcpy.sa.Idw(output_point_feature, "grid_code", cell_size)
out_idw.save(output_raster_idw)
# print("IDW插值完成，输出路径为: " + output_raster_idw)

然后使用ArcPy库中的IDW（反距离加权）插值工具来处理地理空间数据，将点数据插值为栅格格式。先设置输出栅格的像素大小为100个单位，这将影响插值结果的分辨率。然后将当前环境的输出坐标系更改为空间参考编号4548。也就是CGCS2000，这一步是为了把投影坐标系转化为经纬度坐标系。用arcpy.sa.Idw函数把点文件output_point_feature及其属性grid_code进行插值。

# out_idw栅格转换为点
arcpy.conversion.RasterToPoint(output_raster_idw, output_point_feature_idw)
# print("栅格转换为点完成，输出路径为: " + output_point_feature_idw)

栅格转化为点。

field_mappings = arcpy.FieldMappings()
grid_code_field = arcpy.FieldMap()
grid_code_field.addInputField(output_point_feature_idw, "grid_code")
output_field = grid_code_field.outputField
output_field.name = "gridcode"
grid_code_field.outputField = output_field
field_mappings.addFieldMap(grid_code_field)

创建一个FieldMappings对象用于定义字段如何在空间连接中处理。创建FieldMap对象以映射指定字段。这里将output_point_feature_idw中的“grid_code”字段添加到字段映射中，并将输出字段名设置为“gridcode”。最终把这个字段映射添加到FieldMappings对象中。

arcpy.analysis.SpatialJoin(
    target_features=output_polygon,
    join_features=output_point_feature_idw,
    out_feature_class=output_join_result,
    join_type="KEEP_ALL",
    field_mapping=field_mappings
)

使用arcpy.analysis.SpatialJoin将output_polygon和output_point_feature_idw进行空间连接，创建一个新的输出要素类（output_join_result）。join_type="KEEP_ALL"表示保留所有目标要素，即使它们没有匹配的连接要素。

target_coordinate_system = arcpy.SpatialReference(4549)
arcpy.management.DefineProjection(output_join_result, target_coordinate_system)
arcpy.management.Project(output_join_result, output_feature_class, target_coordinate_system)

定义和投影坐标系

area_field = "add"
arcpy.management.AddField(output_feature_class, area_field, "DOUBLE")
arcpy.management.CalculateField(output_feature_class, area_field, "!shape.geodesicarea@squaremeters!", "PYTHON_2.7")

给output_feature_class添加一个名为“add”的双精度字段。接着计算每个要素的测地面积，并将结果存储在“add”字段中

area_field = "Add"
cursor = arcpy.SearchCursor(output_feature_class, [area_field])
area = []
for row in cursor:
    area.append(row.getValue(area_field))
del cursor
area_array = np.array(area)
arcpy.ClearWorkspaceCache_management()
arcpy.env.outputCoordinateSystem = original_coord_system

使用搜索光标（SearchCursor）遍历output_feature_class中的各行，提取“Add”字段的值并储存到一个列表area中。再将列表转换为NumPy数组area_array用于后面的计算，最后return area_array就得到每个地区的面积了。

2.3自动化选择地形进行计算

用上面自动化的过程的函数raster_processing，导入不同的文件再进行计算就可以得到想要的不同地形的各种参数。这里我们写一个可以自动化选择文件进行导入并计算不同地形的函数。

import os
from tqdm import tqdm
import raster_processing
import pandas as pd
import numpy as np

首先我们要定义三个函数

def modification(num):
    average_num = np.mean(num[num != 0])
    zero_indices = np.where(num == 0)
    num[zero_indices] = average_num
    return num


def calculation(num, temperature, rain, NPP):
    formula_A = NPP - (np.exp(0.22) * (1.25 * np.exp(0.05452 * temperature) * rain / (4.259 + temperature)) ** 0.87)
    result = num * (0.6 + (formula_A * 0.01) + 3.8) * 3.6666
    return result


def process_and_calculate(workspace_path, path, temperaturetif, raintif, NPPtif):
    tif_files = [os.path.join(path, file) for file in os.listdir(path) if file.endswith(".tif")]
    for tif_file in tqdm(tif_files, desc="Processing"):
        y = os.path.basename(tif_file)
        calculate(y, workspace_path, tif_file, temperaturetif, raintif, NPPtif)

modification这个函数用于处理数组中的零值。它首先计算数组中非零值的平均数，然后将所有零值替换为平均数。做这一步是因为有的湿地面积太小了，没有被赋予数值，需要这一步把 0值的湿地赋予值。

calculation函数计算碳汇量，使用温度、降水量和NPP数据及公式来进行计算，返回计算后的结果。
process_and_calculate函数负责遍历指定路径下的所有TIF文件，并调用calculate函数进行处理。

并且使用tqdm库，提供进度条显示，方便实时监控处理状态。

def calculate(y, workspace_path, tif_file, temperaturetif, raintif, NPPtif):
    if y == 'river.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, "monthly_temperature.tif")
        Carbon_sink = total_area_result * ((0.6 + 0 + 3.8) * 3.6666 + 0.9 * 3.6666 * 0.1)
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '碳汇量']]
        excel_filename = "river.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))

    elif y == 'lake.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, "monthly_temperature.tif")
        Carbon_sink = total_area_result * ((0.6 + 0 + 3.8) * 3.6666 + 0.9 * 3.6666 * 0.1)
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '碳汇量']]
        excel_filename = "lake.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))

    elif y == 'canal.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, "monthly_temperature.tif")
        Carbon_sink = total_area_result * (0.6 + 0 + 0) * 3.6666
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '碳汇量']]
        excel_filename = "canal.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))
    elif y == 'marsh.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, temperaturetif)
        temperature = modification(
            (raster_processing.process_raster_data(workspace_path, tif_file, temperaturetif))) - 273.1
        rain = modification(raster_processing.process_raster_data(workspace_path, tif_file, raintif)) * 100
        NPP = modification(raster_processing.process_raster_data(workspace_path, tif_file, NPPtif)) * 1000 * 0.08333
        Carbon_sink = calculation(total_area_result, temperature, rain, NPP)
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '气温': temperature,
            '降水量': rain,
            'NPP': NPP,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '气温', '降水量', 'NPP', '碳汇量']]
        excel_filename = "marsh.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))
    elif y == 'swamp.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, temperaturetif)
        temperature = modification(
            (raster_processing.process_raster_data(workspace_path, tif_file, temperaturetif))) - 273.1
        rain = modification(raster_processing.process_raster_data(workspace_path, tif_file, raintif)) * 100
        NPP = modification(raster_processing.process_raster_data(workspace_path, tif_file, NPPtif)) * 1000 * 0.08333
        Carbon_sink = calculation(total_area_result, temperature, rain, NPP)
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '气温': temperature,
            '降水量': rain,
            'NPP': NPP,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '气温', '降水量', 'NPP', '碳汇量']]
        excel_filename = "swamp.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))
    elif y == 'field.tif':
        total_area_result = raster_processing.process_raster_data1(workspace_path, tif_file, temperaturetif)
        temperature = modification(
            (raster_processing.process_raster_data(workspace_path, tif_file, temperaturetif))) - 273.1
        rain = modification(raster_processing.process_raster_data(workspace_path, tif_file, raintif)) * 100
        NPP = modification(raster_processing.process_raster_data(workspace_path, tif_file, NPPtif)) * 1000 * 0.08333
        Carbon_sink = calculation(total_area_result, temperature, rain, NPP)
        serial_numbers = list(range(1, len(total_area_result) + 1))
        df = pd.DataFrame({
            '序号': serial_numbers,
            '面积': total_area_result,
            '气温': temperature,
            '降水量': rain,
            'NPP': NPP,
            '碳汇量': Carbon_sink
        })
        df = df[['序号', '面积', '气温', '降水量', 'NPP', '碳汇量']]
        excel_filename = "field.xlsx"
        df.to_excel(excel_filename, index=False)
        print("Data saved to {}".format(excel_filename))

calculate函数根据传入的文件名来决定执行哪个计算过程。对于特定的TIF文件（如river.tif、lake.tif、canal.tif等），调用raster_processing模块中的process_raster_data1函数处理数据。根据结果计算“碳汇量”（Carbon_sink），并将计算结果、面积等数据存入Pandas DataFrame。将DataFrame导出为Excel文件，文件名称设置为对应的地形。计算公式在之前的流程图里有显示。

具体的处理逻辑如下：
河流地形 (river.tif) 和湖泊地形 (lake.tif)，调用process_raster_data1函数获取总面积，然后计算碳汇量。创建DataFrame，包含序号、面积和碳汇量，并保存为Excel文件。
人工水渠(canal.tif)：
处理过程与河流和湖泊类似，碳汇量的计算公式略有不同。

乔木湿地 (marsh.tif)、草甸湿地 (swamp.tif) 和水田湿地 (field.tif)：
这些文件使用温度、降水量和NPP（净初级生产力）数据进行更复杂的计算。
使用modification和calculation函数来处理这些数据，计算相应的碳汇量，并将结果输出到Excel文件。

2.4脚本文件

workspace = "Run_time_folder"
current_time = datetime.datetime.now()
folder_name = current_time.strftime("%Y%m%d_%H%M%S")
folder_path = os.path.join(workspace, folder_name)
if not os.path.exists(folder_path):
    os.makedirs(folder_path)
workspace_path = folder_path
path = "Run_time_folder/tif"

temperature_tif = None
for filename in os.listdir("temperature_tif"):
    if filename.endswith(".tif"):
        temperaturetif = os.path.join(folder_path, filename)
        break
rain_tif = None
for filename in os.listdir("rain_tif"):
    if filename.endswith(".tif"):
        temperaturetif = os.path.join(folder_path, filename)
        break
NPP_tif = None
for filename in os.listdir("NPP_tif"):
    if filename.endswith(".tif"):
        temperaturetif = os.path.join(folder_path, filename)
        break

Calculate_tif.process_and_calculate(workspace_path, path, temperature_tif, rain_tif, NPP_tif)

最后写一个根据文件名找到相关文件并执行整个流程的脚本文件。这样就完成了整个步骤的串联，但是其实可以把它写成一个执行文件，比如exe之类的，用pyinstaller库，这样会更加便捷。当时时间比较紧张所以就没有继续，感兴趣的小伙伴可以尝试一下，或者有什么问题可以问我。项目完整文件我放在我github库里了https://github.com/Kukeer2/Using-the-Arcpy-Package-Included-with-ArcGIS-ArcMap-for-Simple-Python-Script-Automation.git