將地理範圍網格化以及高德Loca API使用

"網格分析在不少層面都能用到，近期的接觸當中，發現有對區域地價作分析的需求。針對該狀況，筆者找了相關的資料實現了一個範圍的網格化。同時用生成的數據，調用高德的開放接口，嘗試了蜂窩熱力圖。該基礎功能，能爲後續作數據增值服務作準備。"html

視頻：

01前端

—python

範圍線準備jquery

範圍線能夠稍微大一點，能夠用shp，json,kml格式的數據。若是找不到合適的，最近發現一個下載到縣區一級的網站，能夠下載一個範圍線。web

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

02api

—微信

轉化爲列表markdown

用Python讀取js數據，把經緯度提取出來轉化爲列表格式。app

範圍文件：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)]

—

經緯度投影

地理座標要先轉化成平面投影，在平面投影裏作幾何運算。

下面代碼包含了轉極座標與存儲數據爲js格式：util.py

import json
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]

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)))

—

獲取正多邊形

這裏主要是選擇不一樣的多邊形，正三角，正方形，正六邊形。用到shaply模塊。

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]

—

正多邊形填充

用正多邊形填充範圍，這裏面用BFS（廣度優先搜索）進行填充，以範圍的幾何中心往外圍填充。

poly_fill:

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

—

生成數據成果

引入上述的幾個python文件，最終生成js文件，方便前端引入使用。

main.py:

from pathlib import Path
from shapely.geometry import Polygon

from data import boundary_district_330106from poly_fill import PolyFillfrom util import project_to_plane, project_to_polar, to_js

def 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米

—

前端代碼

這裏前端代碼，仍是高德的api,還須要本身寫部分js,拆分的很開!

index.html:

<!DOCTYPE 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-seg</title> <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>

—

高德蜂窩熱力圖

高德地圖能夠直接從網址中複製代碼，輸入本身的密鑰。網址以下：

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

要將風格一行進行調整，換成黑色風格底圖，數據是tsv格式的，提取蜂窩的中心地理座標。數據的值統一用1000.

gaode.html:

<!DOCTYPE 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>蜂窩熱力圖（按米分箱）- 3D</title> <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

後臺回覆關鍵詞：教育，可查看個人我的在線課程，歡迎你們來溝通、交流、共同窗習！

本文分享自微信公衆號 - 玩大數據的規劃師（paitashuju）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。