对点图层数据添加网格编号，实现分组……

craybb

于 2024-08-19 09:20:24 发布

阅读量311

点赞数 10

分类专栏：解决方案文章标签：分组蜂巢多边形网格

本文链接：https://blog.csdn.net/craybb/article/details/141310921

版权

解决方案专栏收录该内容

8 篇文章 0 订阅

订阅专栏

对点数据添加网格编号

一、背景需求

将点数据按规则格网分组的原因很多，如制图中对点位进行抽稀，并行处理时对点要素集合进行分组等。或许你说，我直接使用点要素的ID对电脑CPU核数取余，也能实现分组，但当这些点位之间需要考虑临近关系时，就得按一定的规则对点位进行分组。常用的方法，就是对点集合进行网格编号。

二、添加网格编号

规则形状格网仅可由等边三角形、正方形或六边形组成，因为只有这三种多边形形状可以进行镶嵌（使用相同的形状重复地拼接，以无缝或无重叠覆盖一块区域）以创建均匀的格网。

等边三角形网格编号，使用较少，本文主要介绍正方形和正六边形两种规则网格编号的实现。

三、为何选择六边形

正方形（渔网）格网是 GIS 分析和专题制图中主要使用的形状类型。在制图点位标注放置位置中，可以看出，标注最佳位置放置与正方形邻域所占的格网是对应的。

实际上，在点位制图标注的时候，我们很少将八个方位的区域设置位置权重都设置为非0。主要原因是八个方位都允许放置点位标注信息，最后的制图效果并不佳，图面会呈现过多的标注信息而略显杂乱。而若只选择正东、正南、正北、正西四个方位进行标注放置，标注符号与标注之间的关系都是沿着X轴，Y轴放置，很美观，但是标注内容略显单薄，图面承载内容偏少。

六边形在某些情况下更适合进行分析：

（1）六边形在所有六个方向上到质心的距离相同，如果您要使用距离范围查找邻域，如果使用与渔网格网相对的六边形格网，您会在各要素的计算中得到更多的邻域；

（2）使用六边形格网查找邻域更加直接。由于每个边的接触角度和边长都相同，每个邻域的质心都相等。但是，使用渔网格网时，（上/下/右/左）邻域的质心都距离有 N 单位远，而对角线的质心邻域距离更远（正好等于 2 的平方根乘以 N 单位）；

蜂巢网格：

四、生成蜂巢多边形

生成蜂巢多边形比较简单，重点需考虑奇数行和偶数行正六边形顶点的X轴坐标值计算方法不同。蜂巢多边形的编号字段“row_col”属性值，从1_1，自西南角开始，至东北方向逐行逐列加1。

根据输入要素类图层范围和指定的正六边形边长，生成蜂巢多边形要素类图层。代码实现如下：

from pystd.utils.ReadWrite.geoDB import *
from shapely.geometry import Polygon
import pandas as pd

def _change_prj_type(layer, _org_crs):
    prj_type = _org_crs.type_name
    if prj_type.split(' ')[0] == 'Geographic':
        # 如果是地理坐标系，先转换成Web Mercator
        layer.to_crs(epsg='102100', inplace=True)


def _get_extent(layer):
    x_min, y_min, x_max, y_max = layer.total_bounds
    return x_min, y_min, x_max, y_max


def create_hexagon_polygon(in_feature, honeycomb_size, save_path):
    obj = GeoFC()
    layer = obj.set_oid(in_feature, "geo_oid")
    _org_crs = layer.crs
    _change_prj_type(layer, _org_crs)
    minX, minY, maxX, maxY = _get_extent(layer)

    x, y = (minX, minY)
    geom_array = []
    row_num = []
    i = 1
    row = 0
    col = 0
    while y <= maxY:
        col += 1
        while x <= maxX:
            row += 1
            geom = Polygon([(x, y - honeycomb_size),
                            (x + pow(3, 0.5) / 2 * honeycomb_size, y - 0.5 * honeycomb_size),
                            (x + pow(3, 0.5) / 2 * honeycomb_size, y + 0.5 * honeycomb_size),
                            (x, y + honeycomb_size),
                            (x - pow(3, 0.5) / 2 * honeycomb_size, y + 0.5 * honeycomb_size),
                            (x - pow(3, 0.5) / 2 * honeycomb_size, y - 0.5 * honeycomb_size),
                            (x, y - honeycomb_size)])

            x += pow(3, 0.5) * honeycomb_size
            geom_array.append(geom)
            row_num.append('{}_{}'.format(col, row))

        row = 0
        i += 1
        y += 1.5 * honeycomb_size

        if i % 2 == 0:
            x = minX + pow(3, 0.5) / 2 * honeycomb_size
        else:
            x = minX

    df = pd.DataFrame({"row_col": row_num})

    gdf = gpd.GeoDataFrame(df, geometry=geom_array, crs='epsg:102100')

    gdf.to_crs(_org_crs, inplace=True)

    out_path = os.path.join(save_path, "hexagon_grid_num" + ".shp")
    gdf = gdf.loc[:, ['geometry', 'row_col']]
    gdf.to_file(out_path,
                driver='ESRI Shapefile',
                encoding='utf-8')

if __name__ == "__main__":
    in_feature = r"D:\3.测试数据\test.shp"
    save_path = r"D:\6.scratch_ws"
    cell = 10000

    # 生成正六边形
    create_hexagon_polygon(in_feature, cell, save_path)

五、添加网格编号

网格编号规则与蜂巢多边形编号的规则相同。都是从西南开始（minX, minY），至东北（maxX, maxY）结束。本文仅展示了正六边形编号的效果图，正方形编号比较简单，可在下文给出的“添加网格编号”代码中的函数add_grid_number中，将参数type的值改为“square”即可，如add_grid_number(in_feature, cell, type="square")。

根据输入的点要素类图层，输出为点位在指定边长的正方形或正六边形中所在的编号。函数“add_grid_number”用在了制图工具之---POI显示抽稀中。添加网格编号实现代码如下：

import math
from shapely.geometry import Point, Polygon
from pystd.utils.ReadWrite.geoDB import *


def even_swatch_dict(row, col, key):
    swatch_dict = {
        0: f"{row - 1}_{col}",
        1: f"{row}_{col - 1}",
        2: f"{row + 1}_{col}",
        3: f"{row + 1}_{col + 1}",
        4: f"{row}_{col + 1}",
        5: f"{row - 1}_{col + 1}"
    }
    return swatch_dict.get(key)


def odd_swatch_dict(row, col, key):
    swatch_dict = {
        0: f"{row - 1}_{col - 1}",
        1: f"{row}_{col - 1}",
        2: f"{row + 1}_{col - 1}",
        3: f"{row + 1}_{col}",
        4: f"{row}_{col + 1}",
        5: f"{row - 1}_{col}"
    }
    return swatch_dict.get(key)


def _change_prj_type(layer, _org_crs):
    prj_type = _org_crs.type_name
    if prj_type.split(' ')[0] == 'Geographic':
        # 如果是地理坐标系，先转换成Web Mercator
        layer.to_crs(epsg='102100', inplace=True)


def calculate_hexagon_grid_number(p, honeycomb_size, minX, minY):
    dx = p.x - minX
    dy = p.y - minY
    # 初略计算点所在正六边形的行和列
    if dx == 0:
        col = 1
    else:
        col = math.ceil(dx / (pow(3, 0.5) * honeycomb_size))

    if dy == 0:
        row = 1
    else:
        row = math.ceil(dy / (1.5 * honeycomb_size))

    # 根据行列号，计算中心点坐标
    center_y = (row - 1) * 1.5 * honeycomb_size + minY
    if row % 2 == 0:
        center_x = pow(3, 0.5) * honeycomb_size * (col - 0.5) + minX
    else:
        center_x = pow(3, 0.5) * honeycomb_size * (col - 1) + minX

    geom = [(center_x, center_y - honeycomb_size),
            (center_x - pow(3, 0.5) / 2 * honeycomb_size, center_y - 0.5 * honeycomb_size),
            (center_x - pow(3, 0.5) / 2 * honeycomb_size, center_y + 0.5 * honeycomb_size),
            (center_x, center_y + honeycomb_size),
            (center_x + pow(3, 0.5) / 2 * honeycomb_size, center_y + 0.5 * honeycomb_size),
            (center_x + pow(3, 0.5) / 2 * honeycomb_size, center_y - 0.5 * honeycomb_size),
            (center_x, center_y - honeycomb_size)]

    # 判断当前点位与正六边形的关系
    hexagon = Polygon(geom)
    # 判断点P是否在当前正六边形内
    if hexagon.contains(p):
        grid_num = "{}_{}".format(row, col)
    else:
        distance = []
        for i in range(len(geom) - 1):
            center_seg = (geom[i][0] + geom[i + 1][0]) / 2, (geom[i][1] + geom[i + 1][1]) / 2
            center_p = Point(center_seg)
            dis = center_p.distance(p)
            distance.append(dis)
        min_index = distance.index(min(distance))
        if row % 2 == 0:
            grid_num = even_swatch_dict(row, col, min_index)
        else:
            grid_num = odd_swatch_dict(row, col, min_index)

    return grid_num


def _get_extent(layer):
    x_min, y_min, x_max, y_max = layer.total_bounds
    return x_min, y_min, x_max, y_max


def calculate_square_grid_number(p: Point, minX, minY, cell) -> str:
    x, y = p.x, p.y
    grid_no_x = math.ceil((x - minX) / cell)
    grid_no_y = math.ceil((y - minY) / cell)
    return "{}_{}".format(grid_no_x, grid_no_y)


def add_grid_number(in_feature, cell, type='square'):
    obj = GeoFC()
    layer = obj.set_oid(in_feature, "geo_oid")
    _org_crs = layer.crs
    _change_prj_type(layer, _org_crs)
    x_min, y_min, *_ = _get_extent(layer)
    if type == "square":
        layer.loc[:, 'gridNo'] = layer.geometry.apply(lambda row: calculate_square_grid_number(row, cell, x_min, y_min))
    elif type == "hexagon":
        layer.loc[:, 'gridNo'] = layer.geometry.apply(lambda row: calculate_hexagon_grid_number(row, cell, x_min, y_min))

    layer.to_crs(_org_crs, inplace=True)

    out_path = os.path.join(save_path, "{}_grid_num".format(type) + ".shp")
    layer.to_file(out_path,
                  driver='ESRI Shapefile',
                  encoding='cp936')

if __name__ == "__main__":
    in_feature = r"D:\3.测试数据\test.shp"
    save_path = r"D:\6.scratch_ws"
    cell = 10000

    # 添加正六边形编号
    add_grid_number(in_feature, cell, type="square")

craybb

关注

10
点赞
踩
12

收藏

觉得还不错? 一键收藏
打赏
0
评论
对点图层数据添加网格编号，实现分组……

将点数据按规则格网分组的原因很多，如制图中对点位进行抽稀，并行处理时对点要素集合进行分组等。或许你说，我直接使用点要素的ID对电脑CPU核数取余，也能实现分组，但当这些点位之间需要考虑临近关系时，就得按一定的规则对点位进行分组。常用的方法，就是对点集合进行网格编号。
复制链接

扫一扫