简介这份PDF文档面向地理信息、大数据算法方向的研究者与工程技术人员聚焦“天地图·常州”平台地理数据资源匮乏的现实问题探讨如何借助大数据算法对多源在线地理数据进行解析与聚合。内容从基础测绘数据电子地图、数字正射影像、地面高程模型的局限切入重点分析团购、房产、公交、公众点评等公共服务类数据的结构属性提出面向在线服务数据集与在线文本数据的获取、解析和聚合技术方案并给出国家、省、市级“天地图”节点建设框架下的原型系统设计最终在常州平台上实现多类数据的服务聚合。资源包为1个PDF文件大小约3.36MB适合作为地理信息公共服务、数据融合方向的学习与参考材料。目前已有74人学习可帮助读者理解多源异构地理数据的解析思路、聚合框架与落地路径为相关课题研究或平台建设提供方法借鉴。1. 天地图常州数据解析与聚合从瓦片到业务图层的那条链路拿到「面向天地图常州的地理数据解析与聚合方法研究」这个题目时很多人第一反应是去翻论文但真正落到工程里它讲的其实是一件很具体的事天地图常州节点对外提供的是标准 WMTS 瓦片服务和若干矢量接口而业务侧要的往往是「常州范围内所有学校周边 500 米有哪些餐饮 POI」这类聚合结果。中间隔着坐标解析、瓦片拼接、属性抽取、空间聚合四道工序任何一道没处理好最后出来的图层要么偏移几百米要么聚合结果对不上。这篇笔记面向两类人一类是手里有天地图常州底图、想把业务数据叠上去做空间分析的 GIS 开发另一类是做大数据方向毕业设计、需要一套能跑通的地理数据解析与聚合链路的同学。核心不是讲天地图怎么申请 key而是讲清楚从瓦片 URL 到可聚合数据集之间参数怎么设、坐标系怎么对齐、聚合算法怎么选。坐标系和瓦片层级这两件事是后面所有步骤的地基地基歪了聚合做得再花哨也是白搭。2. 天地图常州瓦片服务的解析URL 结构、层级与坐标对齐2.1 瓦片 URL 的构成与常州节点的层级选择天地图的标准 WMTS 瓦片 URL 通常形如http://t{0-7}.tianditu.gov.cn/DataServer?Tvec_wx{x}y{y}l{z}tk你的key其中vec_w是矢量底图经纬度投影img_w是影像底图cva_w是注记层。常州位于东经 119°08′ 到 120°12′、北纬 31°09′ 到 32°04′ 之间按 Web Mercator 切分在 z10 时常州大约占 2×2 块瓦片z14 时能看清街道z16 以上才有建筑轮廓。实际做解析时第一个要定的是层级。层级不是越高越好z18 的单块瓦片覆盖范围约 150 米常州全域需要几十万块瓦片下载和拼接成本极高而 z12 虽然只有几十块但聚合到街道级别时精度不够。我一般会按业务粒度反推——如果聚合单元是社区z14 足够如果是建筑物级别才上到 z17。# 计算常州某点在某层级下的瓦片行列号Web Mercator import math def lonlat_to_tile(lon, lat, zoom): lat_rad math.radians(lat) n 2.0 ** zoom x_tile int((lon 180.0) / 360.0 * n) y_tile int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n) return x_tile, y_tile # 常州市中心近似坐标 print(lonlat_to_tile(119.974, 31.811, 14)) # 输出该层级下的瓦片行列号这段代码的关键在math.asinh(math.tan(lat_rad))它是 Web Mercator 纬度投影的标准公式不能用简单的线性换算替代否则在高纬度会明显偏移。常州纬度不高误差相对小但如果你的研究区往北延伸到 40° 以上这个公式必须用对。zoom参数直接决定瓦片数量每加一级瓦片总数翻四倍这是后面做聚合时要算清楚的成本账。2.2 坐标系对齐WGS84、GCJ02 与天地图的实际输出天地图官方服务输出的是 CGCS2000 坐标系日常可以近似当作 WGS84 处理。但问题在于很多业务数据来自不同渠道手机采集的 GPS 是 WGS84从某些地图 API 拿到的 POI 可能是 GCJ02 偏移后的坐标。如果直接把 GCJ02 的点和天地图瓦片叠在一起常州地区会有 300 到 500 米的系统性偏移这个偏移在 z14 下肉眼可见聚合结果会整体错位。常见做法是统一转到 WGS84 再做聚合。GCJ02 转 WGS84 没有精确的解析公式工程上用迭代逼近或查表法。下面是一个简化的迭代反解精度在米级够用# GCJ02 转 WGS84 的迭代逼近简化版米级精度 def gcj02_to_wgs84(lon, lat): dlon, dlat lon - 119.974, lat - 31.811 # 以常州中心为参考 # 实际工程中应使用完整的偏移量计算此处示意迭代思路 for _ in range(5): # 正向计算 GCJ02 偏移量并修正 lon - dlon * 0.000001 lat - dlat * 0.000001 return lon, lat提示上面是示意迭代结构真实项目里建议直接用成熟的开源转换库不要自己手写偏移公式否则边界情况容易翻车。坐标对齐做完后建议做一次可视化验证把转换后的点和天地图 z14 瓦片叠在一起看常州火车站、文化宫这些地标是否吻合。如果整体偏移一个固定方向说明坐标系没转对如果局部扭曲说明数据本身有问题。2.3 瓦片下载与拼接的工程参数下载瓦片时并发数和请求间隔是两个必须调的参数。天地图对未认证的 key 有频率限制并发太高会被临时封禁。我一般设并发 4 到 8每次请求间隔 100 到 200 毫秒单机一天能稳定拉几万块瓦片。拼接时用 Pillow 或 GDAL 都行GDAL 更适合大范围拼接因为它支持分块读写不会一次性把几十 GB 的图全读进内存。# 用 gdal_translate 把拼接后的瓦片转成 GeoTIFF带上地理参考 gdal_translate -of GTiff -a_srs EPSG:3857 \ -a_ullr 119.0 32.1 120.2 31.0 \ merged_tiles.png changzhou_base.tif-a_ullr四个参数分别是左上经度、左上纬度、右下经度、右下纬度必须和你下载的瓦片范围严格对应差一点整个图层就偏了。-a_srs EPSG:3857指定 Web Mercator 投影如果后续要转成 WGS84 做空间分析再用gdalwarp转一次。3. 地理数据解析从瓦片到结构化属性的抽取方法3.1 矢量数据解析GeoJSON 与 Shapefile 的字段映射天地图常州除了瓦片还可能提供行政区划、道路、POI 等矢量接口返回格式常见的是 GeoJSON。解析 GeoJSON 的核心是把properties里的字段映射到业务表结构。比如 POI 数据里name、address、type三个字段是聚合时最常用的type字段的编码体系要和你的业务分类对齐否则聚合出来的类别会乱。import json def parse_geojson_features(filepath): with open(filepath, r, encodingutf-8) as f: data json.load(f) records [] for feat in data[features]: props feat[properties] geom feat[geometry] records.append({ name: props.get(name, ), address: props.get(address, ), category: props.get(type, unknown), lon: geom[coordinates][0], lat: geom[coordinates][1] }) return records这段代码里geom[coordinates]对 Point 类型是[lon, lat]但如果是 Polygon结构会嵌套一层直接取[0]和[1]会报错。解析前要先判断geometry.typePoint、LineString、Polygon 三种类型的坐标结构完全不同。这是新手最容易踩的坑之一。3.2 属性清洗缺失值、重复点与异常坐标的处理从瓦片或接口拿到的原始数据缺失值和重复点很常见。缺失值分两种一种是字段为空可以用默认值填充另一种是坐标缺失这种记录必须丢弃因为没法参与空间聚合。重复点通常是同一 POI 被多次采集按name address去重即可但要注意同名不同址的情况不能只按 name 去重。异常坐标的判定有个简单规则常州范围外的点直接标记为异常。经度不在 119.0 到 120.2、纬度不在 31.0 到 32.1 之间的要么是坐标系没转对要么是数据本身错误。这一步放在聚合之前做能省掉后面很多排查时间。def clean_records(records): seen set() cleaned [] for r in records: # 坐标范围过滤 if not (119.0 r[lon] 120.2 and 31.0 r[lat] 32.1): continue # 去重键 key (r[name], r[address]) if key in seen: continue seen.add(key) cleaned.append(r) return cleanedseen集合用元组做键比字符串拼接更可靠避免name或address里含分隔符导致的误判。清洗后的数据量通常会减少 10% 到 30%具体取决于数据来源的规范程度。3.3 空间索引的建立为聚合做准备数据清洗完下一步是建空间索引。聚合操作本质上是「给定一个中心点找半径 R 内的所有点」如果每次都全表扫描几万条数据还能忍上百万条就不可接受了。常见做法是用 R 树或网格索引。Python 里shapely配合rtree可以快速建索引或者用geopandas的sjoin做空间连接。from shapely.geometry import Point import geopandas as gpd gdf gpd.GeoDataFrame( cleaned, geometry[Point(r[lon], r[lat]) for r in cleaned], crsEPSG:4326 ) # 建立空间索引 gdf.sindexcrsEPSG:4326指定 WGS84 经纬度如果后面要做距离计算建议先投影到 UTM 或 Web Mercator因为经纬度下的距离计算单位是度不是米。常州位于 UTM 50N 带投影后距离计算更直观。4. 数据聚合网格聚合、聚类与业务指标的计算4.1 网格聚合把连续空间切成可统计的单元网格聚合是最常用也最容易解释的方法。把常州全域按固定边长切成网格每个网格统计落入其中的 POI 数量、类别分布等指标。网格边长怎么定如果聚合单元是社区级500 米到 1 公里比较合适如果是商圈分析200 米到 500 米更细。边长越小网格数量越多计算量越大但空间分辨率越高。import numpy as np def grid_aggregate(gdf, cell_size0.005): # cell_size 单位为度0.005 度约 500 米 gdf gdf.copy() gdf[grid_x] (gdf.geometry.x / cell_size).astype(int) gdf[grid_y] (gdf.geometry.y / cell_size).astype(int) result gdf.groupby([grid_x, grid_y]).agg( count(name, size), categories(category, lambda x: list(x)) ).reset_index() return resultcell_size0.005在常州纬度下约等于 500 米这个换算关系是1 度纬度约 111 公里1 度经度在 31.8° 纬度下约 94.5 公里。所以 0.005 度经度约 472 米0.005 度纬度约 555 米网格不是正方形但工程上够用。如果要严格正方形先投影到米制坐标系再切。4.2 聚类聚合DBSCAN 在 POI 热点识别中的参数调优网格聚合的缺点是边界效应——两个相邻网格的点可能属于同一个商圈却被切开了。DBSCAN 能解决这个问题它按点的密度聚类不需要预先指定簇数量。两个核心参数eps是邻域半径min_samples是核心点的最小邻居数。常州 POI 数据下eps设 0.002 度约 200 米、min_samples设 5 到 10能识别出大多数商圈。from sklearn.cluster import DBSCAN coords np.array([[r[lon], r[lat]] for r in cleaned]) db DBSCAN(eps0.002, min_samples5).fit(coords) labels db.labels_ # labels 为 -1 的是噪声点不归属任何簇eps调大簇会合并可能把两个商圈连成一个调小簇会碎一个商圈被拆成多个。min_samples调大噪声点增多小簇被过滤。这两个参数没有万能值建议先用网格聚合的结果做参照看聚类结果和网格热点是否吻合再微调。4.3 聚合指标的计算密度、多样性与可达性聚合不只是数点数。业务上常用的指标有三类密度单位面积内的 POI 数量、多样性类别熵衡量一个区域内业态丰富程度、可达性到最近设施的距离。密度计算简单多样性用香农熵可达性需要做最近邻查询。from scipy.spatial import cKDTree tree cKDTree(coords) # 查询每个点到最近邻的距离 distances, indices tree.query(coords, k2) nearest_dist distances[:, 1] # 排除自身k2是因为查询结果第一个是点自身距离为 0第二个才是最近邻。nearest_dist的单位是度要转成米需要乘以 111000 再按纬度修正。这一步在评估公共服务设施覆盖时特别有用比如算每个小区到最近学校的距离。5. 避坑与排查天地图常州数据解析聚合中的五个血泪教训5.1 瓦片拼接后整体偏移几百米现象拼接后的底图和业务点叠在一起所有点统一往一个方向偏。原因瓦片下载时用的投影和拼接时指定的-a_ullr不匹配或者坐标系没转。解决先确认瓦片是 Web Mercator 还是经纬度直投-a_srs和-a_ullr必须和瓦片实际投影一致再用gdalinfo检查输出文件的角点坐标。5.2 聚合结果里出现大量空网格现象网格聚合后很多网格的 count 为 0地图上大片空白。原因网格范围设得比数据范围大太多或者cell_size太小导致点分散到过多网格。解决先统计数据的经纬度范围网格范围紧贴数据边界cell_size根据数据密度调整一般保证平均每个网格有 3 到 5 个点。5.3 DBSCAN 聚类把整个常州聚成一簇现象eps设得过大所有点都被归为一个簇。原因eps的单位是度0.01 度在常州约 1 公里如果 POI 密集这个半径会把相邻商圈全连起来。解决eps从 0.001 开始试逐步增大观察簇的数量变化找到数量骤降的拐点。5.4 坐标转换后点落在常州范围外现象GCJ02 转 WGS84 后部分点跑到经度 118 或 121 以外。原因迭代反解不收敛或者原始数据本身就是错误坐标。解决加范围过滤转换后不在常州范围内的点直接丢弃同时检查原始数据的采集来源。5.5 空间索引建立后查询仍然很慢现象建了 R 树索引但半径查询还是慢。原因查询时没有真正用上索引比如用gdf[gdf.distance(point) r]这种写法会退化成全表扫描。解决用gdf.sindex.query()或geopandas.sjoin()确保查询走索引路径。6. 进阶技巧用聚合结果反推天地图瓦片缓存策略聚合做完之后还有一个容易被忽略的优化点瓦片缓存。如果你的应用需要频繁加载常州某几个区域的底图可以把这些区域的瓦片预下载并缓存到本地减少实时请求。缓存策略的核心是确定「热点区域」——聚合结果里 POI 密度最高的网格就是最值得缓存的区域。具体做法取网格聚合结果中 count 排名前 20% 的网格把这些网格覆盖的瓦片层级和行列号算出来批量下载到本地目录用z/x/y.png的结构存储。加载时先查本地缓存没有再请求天地图。这样能把常用区域的加载速度提升一个数量级。def get_hot_tiles(grid_result, zoom14, top_ratio0.2): threshold grid_result[count].quantile(1 - top_ratio) hot grid_result[grid_result[count] threshold] tiles set() for _, row in hot.iterrows(): # 网格中心经纬度转瓦片行列号 lon row[grid_x] * 0.005 0.0025 lat row[grid_y] * 0.005 0.0025 x, y lonlat_to_tile(lon, lat, zoom) tiles.add((zoom, x, y)) return tilestop_ratio0.2表示取前 20% 的热点网格这个比例可以根据缓存空间调整。0.0025是网格半边长用来从网格索引反算中心点经纬度。算出来的(zoom, x, y)三元组就是缓存键下载时按这个结构存文件即可。验证缓存效果的方法很简单在浏览器开发者工具里看瓦片请求的响应时间缓存命中时应该在 10 毫秒以内未命中时取决于网络通常 100 毫秒以上。如果命中率低于 60%说明热点区域选得不对或者缓存层级和实际加载层级不一致。我自己做这类项目时习惯在聚合完成后先画一张热点图肉眼确认热点区域是否符合预期再决定缓存哪些瓦片。这一步多花十分钟后面能省掉很多「为什么加载还是慢」的排查时间。希望帮到你。本文还有配套的精品资源点击获取