基于arcpy实现导出区域内网格中心坐标功能

在进行数据采集的时候经常会用到基于“周边检索”结果的采集：就是利用平台（网站或APP）提供的"附近"检索功能，搜索"某个位置"周边“X千米”范围内的某类信息（例如POI），然后采集搜索出来的结果。这里的“某个位置”就是搜索圆形区域的圆心（搜索中心点），“X千米”指的是搜索半径。

这个搜索中心点的选取至关重要。因为如果选的少了会因为区域覆盖不全导致数据遗漏，选的太多（密）了，会增加搜索次数，影响采集效率。所以如何合理的选择搜索中心点很重要。例如，我们在采集"北京市房山区"内POI信息的时候就曾遇到过这样的问题。刚开始我们使用了"行政区、商圈、加油站"这三类信息点的位置作为搜索中心点，采集下来发现有不少遗漏。仔细检查后发现，原因是由于房山区相对比较偏远，这三类信息点比较少，搜索中心点比较少导致有很多地区覆盖不到，从而造成数据缺失。

一种合理的搜索中心点选取方法：将待采集区域划分成面积想等的若干网格，每个网格的面积由搜索半径确定，然后取网格中心点的坐标作为搜索中心点，这样区域内每个地方都能被覆盖到。

下面介绍基于arcpy实现上面的思路。arcpy是ArcGIS里包含的一个Python地理数据分析库。在安装完ArcGIS之后就能使用该库了。需要注意的是ArcGIS安装的时候会自带安装一个32位的Python，我们需要使用它自带的这个Python，否则（例如，使用自己安装的64位Python）会找不到arcpy库，或者出现因为和64位版本Python不兼容导致的异常问题。

第一步，获取待采集行政区的边界坐标。关于行政区的边界坐标获取方法，网上介绍有各种途径，这里推荐一种最简单的方法，使用阿里云datav里的工具，链接是http://datav.aliyun.com/tools/atlas/。
第二步，根据待采集区域的边界坐标，画出该区域范围（多边形）。如下图所示是使用arcpy根据"北京房山区"的边界坐标，绘制出的多边形区域。

第三步，根据待采集区域边界上极限（最大最小）坐标，计算出每个网格（正方形）的顶点坐标，画出网格（渔网图）。这里的网格大小根据搜索半径确定，如果搜索半径为2KM，这里网格边长就选用2KM（近似等于0.009 * 2 经度）大小。如下图所示是在“房山区”上画出的网格后的效果。

第三步，根据待采集区域边界上极限（最大最小）坐标，计算出每个网格（正方形）的顶点坐标，画出网格（渔网图）。这里的网格大小根据搜索半径确定，如果搜索半径为2KM，这里网格边长就选用2KM（近似等于0.009 * 2 经度）大小。如下图所示是在“房山区”上画出的网格后的效果。

"北京市房山区"共被拆分为738个"2KM*2KM"的网格，最后导出的网格中心点坐标列表如下所示。

上述过程的完整代码如下：

# coding: utf-8
# create_boundary_fishnet_coords.py
# 导出行政区边界内渔网格中心点坐标


import sys
import os
import math
import arcpy
import re

km = 3
# 网格大小, 经度0.009度相当于1公里
GRID_WIDTH = 0.009 * km


def create(input_boundary_file):
    # 输出目录
    # output_dir = os.path.splitext(os.path.basename(input_boundary_file))[0] + '-output'
    # output_dir =
    output_dir = re.match('.*\d\\\\',input_boundary_file).group()
    output_dir = output_dir + 'km' +str(km) +'\\'
    print(output_dir)
    if not os.path.exists(output_dir):
        os.mkdir(output_dir)

        # 加载边界原始数据
    bounday_file_data = ''
    with open(input_boundary_file, 'r') as f:
        bounday_file_data = f.read()

        # 根据边界点创建行政区多边形面
    bounday_array = arcpy.Array()
    xmin, ymin, xmax, ymax = None, None, None, None
    is_first = True
    for xy in bounday_file_data.split(';'):
        x, _, y = xy.partition(',')
        x = float(x.strip())
        y = float(y.strip())
        if is_first:
            xmin = xmax = x
            ymin = ymax = y
            is_first = False
        else:
            if x > xmax:
                xmax = x
            if x < xmin:
                xmin = x
            if y > ymax:
                ymax = y
            if y < ymin:
                ymin = y
        bounday_array.add(arcpy.Point(x, y))
        # https://pro.arcgis.com/zh-cn/pro-app/arcpy/classes/polygon.htm
    bounday_polygon = arcpy.Polygon(bounday_array)
    # 导出边界多边形的shp文件,用gis软件(e.g. OpenJUMP)查看
    shp_file = '{}/boundary.shp'.format(output_dir)
    arcpy.CopyFeatures_management(bounday_polygon, shp_file)
    print('Shapefile "{}" is ready.'.format(shp_file))

    # 画出渔网图
    # 根据边界坐标经纬度最大和最小值，依次计算出每个网格正方形四个顶点的坐标
    # 计算网格的行列数
    grid_rows_num = int(math.ceil((ymax - ymin) / float(GRID_WIDTH)))
    grid_columns_num = int(math.ceil((xmax - xmin) / float(GRID_WIDTH)))
    # 依次计算各网格(0, 0), (0, 1), (0, 2) ... (grid_rows_num - 1, grid_columns_num-1)四个顶点的坐标
    grids = []
    for r in range(grid_rows_num):
        for c in range(grid_columns_num):
            grid_4coords = arcpy.Array()
            # 左上角坐标
            x_lt = xmin + c * GRID_WIDTH
            y_lt = ymax - r * GRID_WIDTH
            # 右上角坐标
            x_rt = x_lt + GRID_WIDTH
            y_rt = y_lt
            # 左下角坐标
            x_lb = x_lt
            y_lb = y_lt - GRID_WIDTH
            # 右下角坐标
            x_rb = x_rt
            y_rb = y_lb
            # 按"左上->右上->右下->左下->左上"顺序画一个封闭四边形，注意顺序不能乱，否则画出来的图形不对（我第一次的时候就画成两个对三角了）
            grid_4coords.add(arcpy.Point(x_lt, y_lt))
            grid_4coords.add(arcpy.Point(x_rt, y_rt))
            grid_4coords.add(arcpy.Point(x_rb, y_rb))
            grid_4coords.add(arcpy.Point(x_lb, y_lb))
            grid_4coords.add(arcpy.Point(x_lt, y_lt))
            # 创建一个网格（四边形）
            grids.append(arcpy.Polygon(grid_4coords))
            # 导出网格的shp文件
    shp_file = '{}/grids.shp'.format(output_dir)

    arcpy.CopyFeatures_management(grids, shp_file)
    print('Shapefile "{}" is ready.'.format(shp_file))

    # 对比每个"网格四边形"和"区域多边形"，找到两者"相交"的网格，导出对应的网格中心点坐标
    # 如果 disjoint 返回 False，则两个几何相交，详见 https://pro.arcgis.com/zh-cn/pro-app/arcpy/classes/polygon.htm
    center_points = []
    for grid in grids:
        if bounday_polygon.disjoint(grid) == False:
            # 网格左上角表座
            grid_x_lt = grid.getPart(0)[0].X
            grid_y_lt = grid.getPart(0)[0].Y
            # 计算出网格中心点坐标
            grid_x_center = grid_x_lt + 0.5 * GRID_WIDTH
            grid_y_center = grid_y_lt - 0.5 * GRID_WIDTH
            print(grid_x_center, grid_y_center)
            center_points.append((grid_x_center, grid_y_center))
            # 导出中心点的shp文件
    shp_file = '{}/center_points.shp'.format(output_dir)
    # 导出中心点坐标到csv文件
    csv_file = '{}/center_points.csv'.format(output_dir)
    writer = open(csv_file, 'w')
    writer.write('"longitude","latitude"\n')
    # with open(csv_file,'w') as f:
    #     f.write('"longitude","latitude"\n')
    center_points_geo = []
    for p in center_points:
        center_points_geo.append(arcpy.PointGeometry(arcpy.Point(p[0], p[1])))
        writer.write('"{}","{}"\n'.format(p[0], p[1]))
    writer.close()
    print('CSV file "{}" for center points is ready.'.format(csv_file))
    # 参考 https://gis.stackexchange.com/questions/16122/creating-shapefile-from-lat-long-values-using-arcpy
    arcpy.CopyFeatures_management(center_points_geo, shp_file)
    print('Shapefile "{}" is ready.'.format(shp_file))


if __name__ == '__main__':
    create('H:\\Spider-Projects\\areaMeshing\\boundaryCoordinates\\510100\\coordinatesTo1629529673.txt')

文章转载自：http://www.site-digger.com/html/articles/20201012/817.html

附加：参考原文中的代码是python2.7的，本文中的代码支持python3.目前ArcGIS Pro已经支持3.6了，但是ArcGIS Pro必须激活才能使用（包括内置的python）。
下面放出ArcGIS Pro无限试用puding：
ArcGISPro.2.5 and all extension Patch.rar
效果如图：