使用ArcPy一键生成指定研究区范围的矢量格网数据

2025.01.02

格网作为研究区单元的一种形式，对于很多定量的空间分析有非常大的帮助。但是要制作指定研究区范围的格网数据却并不是很简便。以ArcGIS为例，一般来说，首先是根据自己的单元范围需要创建渔网，接着将创建的渔网与研究区数据以“INTERSECT”空间连接，然后再通过空间连接后的字段删除那些没有连接的字段。操作并不困难，只是有点繁琐。如果需要短时间做多个格网，难免会有点烦躁，所以做了个ArcPy代码，来一劳永逸这个问题。

事前准备

Python库

这个代码需要Arcpy。

工作目录

可以自主指定输入数据和输出数据的地址，数据格式均为shpfile。

代码

import arcpy
import os

# ========== 用户自定义参数 ==========
# 研究区路径
study_area_shp = r"C:\Address\studyarea.shp"

# 输出格网路径
output_grid_shp = r"C:\Address\output.shp"

# 输出格网点路径
output_grid_points_shp = r"C:\Address\output_points.shp"

# 格网大小（单位：米）
grid_size = 200

# 格网坐标系
output_spatial_reference = arcpy.SpatialReference(32650)  # WGS_1984_UTM_Zone_50N

# ========== 脚本逻辑 ==========

try:
    # 确保输出目录存在
    output_dir = os.path.dirname(output_grid_shp)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 设置环境
    arcpy.env.overwriteOutput = True

    # 获取研究区的范围
    study_area_extent = arcpy.Describe(study_area_shp).extent

    # 获取范围坐标
    x_min = study_area_extent.XMin
    y_min = study_area_extent.YMin
    x_max = study_area_extent.XMax
    y_max = study_area_extent.YMax

    # 向外扩展范围以确保格网覆盖完整的研究区
    x_min = int(x_min // grid_size) * grid_size
    y_min = int(y_min // grid_size) * grid_size
    x_max = (int(x_max // grid_size) + 1) * grid_size
    y_max = (int(y_max // grid_size) + 1) * grid_size

    # 创建临时格网路径
    temp_grid = os.path.join(output_dir, "temp_grid.shp")

    # 创建格网
    arcpy.management.CreateFishnet(
        out_feature_class=temp_grid,
        origin_coord=f"{x_min} {y_min}",
        y_axis_coord=f"{x_min} {y_min + 1}",
        cell_width=grid_size,
        cell_height=grid_size,
        number_rows="",
        number_columns="",
        corner_coord=f"{x_max} {y_max}",
        labels="LABELS",
        geometry_type="POLYGON"
    )

    # 创建格网点
    temp_points = os.path.join(output_dir, "temp_points.shp")
    arcpy.management.FeatureToPoint(temp_grid, temp_points, "INSIDE")

    # 使用空间连接工具保留完整正方形格网并移除不相交的单元
    intersected_grid = os.path.join(output_dir, "intersected_grid.shp")
    arcpy.analysis.SpatialJoin(
        temp_grid, 
        study_area_shp, 
        intersected_grid, 
        join_type="KEEP_COMMON", 
        match_option="INTERSECT"
    )

    # 设置坐标系并保存最终格网
    arcpy.management.DefineProjection(intersected_grid, output_spatial_reference)
    arcpy.management.CopyFeatures(intersected_grid, output_grid_shp)

    # 使用空间连接提取与最终格网相交的格网点
    intersected_points = os.path.join(output_dir, "intersected_points.shp")
    arcpy.analysis.SpatialJoin(
        temp_points, 
        intersected_grid, 
        intersected_points, 
        join_type="KEEP_COMMON", 
        match_option="INTERSECT"
    )

    # 设置坐标系并保存最终格网点
    arcpy.management.DefineProjection(intersected_points, output_spatial_reference)
    arcpy.management.CopyFeatures(intersected_points, output_grid_points_shp)

    # 删除临时文件
    arcpy.management.Delete(temp_grid)
    arcpy.management.Delete(temp_points)
    arcpy.management.Delete(intersected_grid)
    arcpy.management.Delete(intersected_points)
    arcpy.management.Delete(temp_grid_label)

    print(f"格网已成功生成并保存到: {output_grid_shp}")
    print(f"格网点已成功生成并保存到: {output_grid_points_shp}")

except Exception as e:
    print(f"发生错误: {e}")