用正六边形分割地图

给定地图上的区域(用多边形顶点的经纬度表示), 需要用正多边形(三角形/正方形/六边形)对地图上的区域进行填充. 在一些实际应用中, 这样做的是为了划分标准的作业单元, 例如定义共享自行车的骑行范围, 物流的配送范围, 外卖骑手的接单范围等.

效果如下(杭州市西湖区: 六边形半径1000米).

image

问题描述

给定地图区域边界顶点的集合. 用正多边形对定义的区域进行分割. 为了保证正多边形恰好填满平面, 可选的正多边形为三角形, 四边形, 六边形.

正多边形

给定正边形的中心点和初始角度(如下图所示).

image

用代表其顶点坐标, 我们有:

为了方便后面的计算, 先定义一个PolyGetter的类用来生成正多边形对象.

#poly_getter.py

import math

from shapely.geometry import Polygon, Point


class PolyGetter(object):
    """ 生成正多边形对象
    """

    def __init__(self, radius, k, theta=0):
        self.radius = radius
        self.k = k  # 正多边形的边数
        self.theta = theta  # 起始角度: degree

    def from_center(self, center):
        """ 输入中心点的坐标，返回对应的正多边形
        :param center: Point对象
        """

        def get_xy(i):
            x = center.x + self.radius * math.cos(2 * math.pi * (i / self.k + self.theta / 360))
            y = center.y + self.radius * math.sin(2 * math.pi * (i / self.k + self.theta / 360))
            return x, y

        return Polygon([Point(get_xy(i)) for i in range(self.k)])

    def from_vertex(self, vertex, i):
        """ 输入顶点的坐标，返回对应的多边形
        :param vertex: 顶点坐标，Point对象
        :param i: 顶点的编号（按极坐标顺序编号）
        """
        c_x = vertex.x - self.radius * math.cos(2 * math.pi * i / self.k + self.theta)
        c_y = vertex.y - self.radius * math.sin(2 * math.pi * i / self.k + self.theta)
        return self.from_center(Point(c_x, c_y))

    def neighbors_of(self, poly):
        """ 输入正多边形，返回它所有邻接的多边形
        :param poly: 多边形，Polygon对象
        """
        dist = self.radius * math.cos(math.pi / self.k)  # 计算中心到边的距离
        p = PolyGetter(2 * dist,
                       self.k,
                       self.theta + 180 / self.k)
        centers = list(p.from_center(poly.centroid).exterior.coords)
        return [Polygon(self.from_center(Point(c))) for c in centers]

坐标转换

考虑到地球是球面, 我们需要把经纬度坐标投影到平面, 然后在屏幕上按照指定的单位距离对地图区域进行分割. 下面两个函数可以分别把坐标点投影到平面以及映射到球面.

from pyproj import Proj


def project_to_plane(points):
    """ 把极坐标点投影到平面.
    """
    p = Proj(4508)
    return [p(point[0], point[1]) for point in points]


def project_to_polar(points):
    """ 把平面上的点转换成极坐标.
    """
    p = Proj(4508)

    def proj_and_round(point):
        q = p(point[0], point[1], inverse=True)
        return round(q[0], 6), round(q[1], 6)

    return [proj_and_round(point) for point in points]

填充

多边形填充的算法如下: 给定初始的多边形做宽度优先搜索(Breadth-First-Search)它的邻接多边形, 直到搜索的区域覆盖整个地图边界. 注意两点: 1. 已搜索的多边形需要标记; 2. 超出边界的多边形不要保存在结果中.

from shapely.geometry import Polygon, MultiPolygon

from poly_getter import PolyGetter


class PolyFill(object):
    """ 给定一个多边形区域, 用正多边形填充它.
    """

    def __init__(self, boundary, center):
        """
        :param boundary: 被分割的多边形, Polygon对象
        :param center: 锚点，以center为出发点进行分割, Point对象
        """
        self._boundary = boundary
        self._center = center
        self._radius = None
        self._k = None
        self._theta = None
        self._poly_getter = None
        self._result = []
        # 用来保存已经搜索过的正多边形（用中心点来表示）
        self._searched_polys = set({})

    def set_params(self, radius, k=6, theta=0):
        """ 参数设置.
        :param radius: 正多边形外接圆的半径
        :param k: 正k多边形.  k = 3, 4, 6
        :param theta: 正多边形的起始角度（度数）
        """
        self._radius = radius
        self._k = k
        self._theta = theta
        assert radius > 0 and k in {3, 4, 6}, \
            ValueError('radius > 0 and k in (3, 4, 6)')
        self._radius = radius
        self._k = k
        self._theta = theta
        self._poly_getter = PolyGetter(self._radius, self._k, self._theta)

        return self

    def run(self):
        """
        :return:
        [[(x11, y11), (x12, y12) ...]  # 多边形1
         [(x21, y21), (x22, y22) ...]  # 多边形2
         ...]                          # ...
        """
        assert self._radius, ValueError("set parameters first!")

        start_poly = self._poly_getter.from_center(self._center)
        # 以start_poly为起点执行BFS填充boundary
        self._fill(start_poly)
        return [list(poly.exterior.coords) for poly in self._result]

    def _fill(self, start_poly):
        """ 给定初始的填充多边形, 按照BFS的方式填充周围的区域.
        以k多边形为例(k=3,4,6), 有 360/k 个多边形与它相邻.
        """
        self._mark_as_searched(start_poly)
        q = [start_poly]
        while len(q):
            poly = q.pop(0)
            self._append_to_result(poly)
            # 把有效的多边形加入队列. 有效的定义:
            # 1. 与poly邻接;
            # 2. 未被搜索过;
            # 3. 在边界内（与boundary定义的区域有交集)
            q += self._get_feasible_neighbors(poly)

    def _mark_as_searched(self, poly):
        """ 把多边形标记为'已搜索'
        """
        self._searched_polys.add(self._get_poly_id(poly))

    @staticmethod
    def _get_poly_id(poly):
        """ 用多边形的中心点的位置判断两个多边形是否相同.
        注意浮点数精度问题.
        """
        c = poly.centroid
        return '%.6f,%.6f' % (c.x, c.y)

    def _is_searched(self, poly):
        """ 判断多边形是否存在
        """
        return self._get_poly_id(poly) in self._searched_polys

    def _append_to_result(self, poly):
        """ 把正多边形poly保存到结果集.
        """
        # poly与boundary取交集, 然后保存结果
        s = self._boundary.intersection(poly)
        if s.is_empty:
            return
        # Polygon对象则直接保存
        if isinstance(s, Polygon):
            self._result.append(s)
        # MultiPolygon对象则依次把它包含的Polygon对象保存
        elif isinstance(s, MultiPolygon):
            for p in s:
                self._result.append(p)

    def _get_feasible_neighbors(self, poly):
        """ 计算与poly邻接的有效的正多边形, 然后标记为'已搜索'.
        """
        def mark_searched(p):
            self._mark_as_searched(p)
            return p

        def is_feasible(p):
            if self._is_searched(p) or self._boundary.intersection(p).is_empty:
                return False
            return True

        # 1. 仅包含'未被搜索'和'不在界外'的正多边形
        # 2. 把poly所有的feasible多边形标记为'已搜索'
        return [mark_searched(p)
                for p in self._poly_getter.neighbors_of(poly)
                if is_feasible(p)]

可视化

• 地图可视化可以使用高德开放平台 - 数据可视化 JS API
• 完整代码

• 正方形分割(半径2000米)

  
  image

• 三角形分割(半径4000米)

  
  image