对点数据添加网格编号
一、背景需求
将点数据按规则格网分组的原因很多,如制图中对点位进行抽稀,并行处理时对点要素集合进行分组等。或许你说,我直接使用点要素的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")