oracle多边形经纬范围筛选_将地理范围网格化以及高德Loca API使用

"网格分析在很多层面都能用到，近期的接触当中，发现有对区域地价做分析的需求。针对该情况，笔者找了相关的资料实现了一个范围的网格化。同时用生成的数据，调用高德的开放接口，尝试了蜂窝热力图。该基础功能，能为后续做数据增值服务做准备。"

视频：

01

—

范围线准备

范围线可以稍微大一点，可以用shp，json,kml格式的数据。如果找不到合适的，最近发现一个下载到县区一级的网站，可以下载一个范围线。

http://datav.aliyun.com/tools/atlas/#&lat=31.769817845138945&lng=104.29901249999999&zoom=4

02

—

转化为列表

用Python读取js数据，把经纬度提取出来转化为列表格式。

范围文件：redline.json(部分)

{"type": "FeatureCollection", "features": [{"type": "Feature", "geometry": {"type": "LineString", "coordinates": [[113.27443651579965, 27.874858011710145], [113.27443937112358, 27.87483870018461], [113.27444062579825, 27.87483021439111], [113.27444426387372, 27.874805608836787], [113.2745222628382, 27.874665403047846]]}, "properties": {"FID_1": 0}}]}

 转化代码：

import jsonwith open('redline.json','r',encoding='utf8')as fp:    json_data = json.load(fp)listhh=json_data['features'][0]['geometry']['coordinates']tupleaa=[]for i in listhh:    tupleaa.append(tuple(i))print(tupleaa)

转化后的列表：

[(113.27443651579965, 27.874858011710145), (113.27443937112358, 27.87483870018461), (113.27444062579825, 27.87483021439111), (113.27444426387372, 27.874805608836787)]

03

—

经纬度投影

地理坐标要先转化成平面投影，在平面投影里做几何运算。

下面代码包含了转极坐标与存储数据为js格式：util.py

import jsonfrom pyproj import Projdef 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]def to_js(points, path, var_name):    """ 保存成js文件(前端展示).    """    data_js = {        'coordinates': [list(p) for p in points]    }    with open(path, 'w', encoding='utf-8') as f:        f.write('%s = [%s];' % (var_name, json.dumps(data_js)))

04

—

获取正多边形

这里主要是选择不同的多边形，正三角，正方形，正六边形。用到shaply模块。

poly_getter.py:

import mathfrom shapely.geometry import Polygon, Pointclass 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]

05

—

正多边形填充

用正多边形填充范围，这里面用BFS(广度优先搜索)进行填充，以范围的几何中心往外围填充。

poly_fill:

from shapely.geometry import Polygon, MultiPolygonfrom yl.poly_getter import PolyGetterclass 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)]

06

—

生成数据成果

引入上述的几个python文件，最终生成js文件，方便前端引入使用。

main.py:

from pathlib import Pathfrom shapely.geometry import Polygonfrom data import boundary_district_330106from poly_fill import PolyFillfrom util import project_to_plane, project_to_polar, to_jsdef map_seg(radius, k, theta=0):    # 把经纬度投影到二维平面    boundary_plane = Polygon(project_to_plane(boundary_district_330106))    # 用正多边形填充    result = PolyFill(boundary_plane, boundary_plane.centroid).set_params(radius, k, theta).run()    # 把结果转换成极坐标    result = [project_to_polar(poly) for poly in result]    # 保存结果    d = Path('web')    to_js(boundary_district_330106, d / 'data_boundaries.js', 'MS.data.blockBoundaries')    to_js(result, d / 'data_bricks.js', 'MS.data.bricks ')if __name__ == '__main__':    map_seg(100, 6)  # 六边形分割: 半径1000米    # map_seg(2000, 4)  # 四边形分割: 半径2000米    # map_seg(4000, 3)  # 三角形分割: 半径4000米

07

—

前端代码

这里前端代码，还是高德的api,还需要自己写部分js,拆分的很开!

index.html:

<html lang="zh-CN"><head>    <meta charset="utf-8">    <meta http-equiv="X-UA-Compatible" content="IE=edge">    <meta name="viewport" content="width=device-width, initial-scale=1">    <title>map-segtitle>    <script src="https://a.amap.com/Loca/static/dist/jquery.min.js">script>    <script src="https://webapi.amap.com/maps?v=1.4.15&key=c9ef357d99db1c97911d0105a24abb6f">script>    <script src="https://webapi.amap.com/loca?v=1.3.0&key=c9ef357d99db1c97911d0105a24abb6f">script>    <style>        html,        body,        #container {            margin: 0;            padding: 0;            width: 100%;            height: 100%;        }        #layerCtrl {            margin: 1%;            padding: 5px;            width: 6%;            background: #484848;            opacity: 0.5;            color: white;        }style>head><body>    <div id="container" class="container">        <script src="index.js">script>        <div id="layerCtrl" class="layerCtrl">            <p>                <input type="checkbox" id="cboxBricks" checked="checked" onclick="ctrlLayer(layerBricks.layer, this)">                <label for="cboxBricks">显示砖块label>            p>        div>    div>body>html>

08

—

高德蜂窝热力图

高德地图可以直接从网址中复制代码，输入自己的密钥。网址如下：

https://lbs.amap.com/api/loca-api/demos/hexagonlayer/loca_heatmap_hexagon_3d

要将风格一行进行调整，换成黑色风格底图，数据是tsv格式的，提取蜂窝的中心地理坐标。数据的值统一用1000.

gaode.html:

<html lang="zh-CN"><head>    <meta charset="utf-8">    <meta http-equiv="X-UA-Compatible" content="IE=edge">    <meta name="viewport" content="width=device-width, initial-scale=1">    <title>蜂窝热力图(按米分箱)- 3Dtitle>    <style>        html,        body,        #container {            margin: 0;            padding: 0;            width: 100%;            height: 100%;        }style>head><body>    <div id="container" class="container">div>    <script src="//webapi.amap.com/maps?v=1.4.15&key=6ebc532a86456a04773258fdf9921696">script>    <script src="//webapi.amap.com/loca?v=1.3.0&key=6ebc532a86456a04773258fdf9921696">script>    <script src="//a.amap.com/Loca/static/dist/jquery.min.js">script>    <script>        var map = new AMap.Map('container', {            viewMode: '2D',            resizeEnable: true,            pitch: 0,            mapStyle: 'amap://styles/2f4e6f376313300be862ab7cdc8340bd',            features: ['bg', 'road']        });        $.get('test.tsv', function (heatmapData) {            var layer = new Loca.HexagonLayer({                map: map,                fitView: true            });            layer.setData(heatmapData, {                lnglat: function (obj) {                    var val = obj.value;                    return [val['lng'], val['lat']]                },                value: 'count',                type: 'tsv'            });            layer.setOptions({                mode: 'count',                unit: 'meter',                style: {                    color: ['#0868AC', '#43A2CA', '#43A2CA', '#7BCCC4', '#BAE4BC', '#F0F9E8', '#F0F9E8'],                    radius: 100,                    opacity: 0.9,                    gap: 20,                    height: [0, 100000]                }            });            layer.render();        });script>body>html>

代码与提取码奉上：

链接：https://pan.baidu.com/s/1iG85fLkNRme_JMzC7GtaLQ 

提取码：o4sd 

后台回复关键词：教育，可查看我的个人在线课程，欢迎大家来沟通、交流、共同学习！

---