"网格分析在不少层面都能用到，近期的接触当中，发现有对区域地价作分析的需求。针对该状况，笔者找了相关的资料实现了一个范围的网格化。同时用生成的数据，调用高德的开放接口，尝试了蜂窝热力图。该基础功能，能为后续作数据增值服务作准备。"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)]
03
—
经纬度投影
地理坐标要先转化成平面投影，在平面投影里作几何运算。
下面代码包含了转极坐标与存储数据为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)))
04
—
获取正多边形
这里主要是选择不一样的多边形，正三角，正方形，正六边形。用到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]
05
—
正多边形填充
用正多边形填充范围，这里面用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))
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 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米
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-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>
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>蜂窝热力图（按米分箱）- 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源创计划”，欢迎正在阅读的你也加入，一块儿分享。