尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

离线海拔查询工具实现:SRTM数据压缩与坐标系转换实战

离线海拔查询工具实现:SRTM数据压缩与坐标系转换实战 做这个“高程海拔离线实时查询显示”工具起因是两年前我在川西党岭徒步时被在线高程接口坑惨了垭口上彻底没信号手机里所有海拔查询功能全部罢工只能靠手表的气压计估个大概等走到有信号的地方再调在线接口又慢又费流量而且返回的海拔跟旁边测绘桩上的高程差出几十米。回来后我花了不少时间把SRTM等高程数据源、坐标系转换、离线地址搜索和实时显示串成了一条完整链路最终做成了一个体积极小、支持输入地址或经纬度、完全不依赖网络的高程海拔离线查询工具。这篇就把整个实现思路、数据选型、压缩方案和踩过的坑完整拆开讲。文中涉及的东西不复杂适合想自建离线地理工具、或者对DEM数据应用感兴趣的人参考。无论你打算做户外导航、骑行轨迹剖面还是野外调查设备里的高度显示里面关于数据底座的取舍和坐标系的处理应该都能直接用上。1. 为什么非要做离线在线高程查询的三块短板1.1 断网即失灵户外场景恰恰最没信号在线高程查询大部分依赖服务端接口比如Open-Elevation、Google Elevation API或者国内天地图附带的高程服务。这类接口有一个天然悖论真正需要精确海拔的场景——登山、徒步、野外勘测、漂流——大部分都在基站覆盖稀薄的山区或峡谷越需要它的时候它越连不上。我自己实测跑过几条川西和秦岭的路线进山半小时后手机基本处于飞行模式状态。那时候所有在线查海拔的工具全废剩下的只有GPS能提供纯椭球高度误差大到没法看。这也是我决定做离线版本的直接原因海拔数据必须躺在本地存储里查询过程不经过任何网络请求。1.2 在线接口的延迟、额度和成本都不够“实时”就算网络正常在线接口也有三个麻烦。首先是延迟一次请求从DNS解析到数据返回普遍在300毫秒到2秒之间如果做轨迹剖面对一整段路线要发上千次请求那体验就是转圈等进度条。其次是额度公开免费接口通常有每日请求上线Open-Elevation自建服务需要服务器第三方公共服务则不够稳定。第三是连续性边走边看实时海拔时GPS每秒都能给出新坐标如果每秒钟都去请求线上接口流量和电量都吃不消。所以“实时显示”这个词本质上要求查询引擎有足够低的单次响应时间本地方案能做到毫秒级返回在线方案很难。1.3 现有离线方案要么太大要么精度太糙当时市面上能离线查海拔的方案大致有三个档次一是专业GIS软件如Global Mapper、QGIS里自带高程数据能离线查但那是给桌面端用的不可能带进山二是部分离线地图App在下载区域包后附带粗略高程显示精度通常只有几百米网格山谷和山峰经常被抹平三是存一张等高线图的图片手动对精度和体验就更别提了。换句话说市面上缺少一个“数据量可控、查询响应快、既能输坐标又能输地名、且能实时跟随GPS”的轻量离线方案。我想做的就是这个缺口。2. 高程数据底座选型SRTM、ASTER GDEM、ALOS 与 GTOPO 的取舍2.1 四类公开DEM源的基础参数对比离线高程工具的一切都建立在DEM数字高程模型数据上。目前公开免费、可以直接下载用于离线项目的主要有这几套数据源发布方原始分辨率覆盖范围高程精度全球原始体积约SRTM3NASA/USGS3弧秒≈90m56°S~60°N绝对误差约±16mGeoTIFF约120GBSRTM1NASA/USGS1弧秒≈30m56°S~60°N绝对误差约±9mGeoTIFF约1TBASTER GDEM v3METI/NASA1弧秒≈30m83°S~83°N绝对误差约±17m覆盖陆地约900GBALOS AW3D30JAXA1弧秒≈30m82°S~82°N标称RMSE约5m覆盖陆地约800GBGTOPO30USGS30弧秒≈1km全球区域差异大约1.2GBSRTM是航天飞机雷达测高得到的覆盖北纬60°到南纬56°全球绝大部分有人区域都能覆盖且数据空洞比早期版本少。ASTER GDEM是光学立体像对生成的在高纬度和陡峭山区覆盖更全但噪声相对大部分区域有带状条纹。ALOS是日本卫星立体影像生成的垂直精度标称最好但个别雪线附近有云影空洞。GTOPO30是30弧秒的全球粗网格作为全球概览够用细节完全不行。2.2 分辨率、覆盖范围与包体的数学账选型前先算一笔账。地球表面积约5.1亿平方公里其中陆地约1.49亿平方公里。如果做全球30弧秒约1km的GTOPO30概览层大约五百万个格网点每个点2字节整数原始约10MB量级压缩后非常小。但如果是全球90m的SRTM3陆地面积1.49亿平方公里除以每个格网约8100平方米约1840万个格网这里要重新算1.49e14平方米/8100平方米 ≈ 1.84e10 个点。每个点2字节原始约36GB。不对稍微修正一下SRTM3全球原始GeoTIFF实际约120GB这里包含了部分重叠、投影和无效值标记。如果只存陆地且用裸整数格式转成紧凑二进制也在20~30GB量级。全球30m的SRTM1就更夸张原始约1TB。所以“极精简离线包”和“全球30m高精度”在手机端是天然冲突的必须在覆盖范围与精度之间做取舍。2.3 我的最终选型三级分辨率分层方案我最终采用的是三层数据叠加方案既不追求不现实的全量全球高精度又能保证覆盖面和单点精度L0概览层基于GTOPO30重采样为60弧秒约2km的全球网格只存每个格网的高程整数采用每10米一个档的量化启动时载入内存约十余MB。作用是一打开App就有全球任意位置的基础高程不至于查一个海外坐标直接空白。L1区域层基于SRTM3按国家和行政区域切块覆盖国内全部并压缩后约400MB用于国内任意坐标的高精度查询。这也是默认随包一起带的主要数据。L2精扫层针对用户高频使用的省份或徒步热点区域使用ALOS 30m数据制作可选下载包单个省份约200~500MB适合对精度有更高要求的人。这样工具可以实现“启动即离线、全球能出数需要精细时再定向增补”而不是一次性塞一个几十GB的怪兽包。对于大多数人的户外使用场景90m分辨率的SRTM已经能反映山脊、峡谷、垭口的基本形态。3. 全球高程塞进手机自定义二进制的压缩与分块实践3.1 从GeoTIFF到紧凑栅格整数化与差值编码下载到的原始DEM是GeoTIFF或SRTM专用的.hgt格式。.hgt格式本身就是二进制整数数组每个采样点2字节小端存储行列固定这很适合直接食用。但GeoTIFF里常常带浮点、带地理变换矩阵、带无效值标记直接放手机里既浪费空间又慢。我做了一个Python批处理脚本用GDAL统一转成裸的16位有符号整数栅格无数据区域填成一个超大哨兵值比如-32768同时把无效水域和空洞用周边插值补齐。转完之后全局还很大于是做了一步“差值编码”。自然地形的规律是相邻30米到90米格网之间的高差通常很小大部分在±10米以内。因此对每一行做差分将原始高程值替换为与前一点的差值再用变长整数编码存。这样平坦区域大量差值接近0压缩率非常可观。对SRTM3的国内区域我用GDAL转成裸栅格后约2.3GB经过差值编码后再用ZZLIB做第二阶段压缩最终落到约420MB压缩比超过5倍。这个数字不神奇对手头有海量高程数据的朋友是完全可以复现的量级。3.2 四叉树分块、金字塔与懒加载压缩完只是一个瘦身真正决定查询速度的在于分块方式。我没有用传统的大文件单块读取而是按1°×1°经纬度网格切块每个1°格网内再切成64×64采样点的瓦片。每个瓦片独立存储瓦片头记录其边界经纬度、最小最大高程值和数据在校验文件中的偏移量。查询一个坐标时先用坐标索引到对应的1°格网然后从格网索引表里读取目标瓦片的偏移量再通过RandomAccessFile只读取那一个小块通常几KB解压后拿到采样值。这样App启动时只需要加载瓦片索引表约几十MB不需要把全部高程数据读进内存。真实查询时一次随机IO加一次小解压耗时能控制在几毫秒内表现上就是“秒出海拔”。为了规避“用户在边缘地区快速拖动或GPS抖动时要读大量瓦片”的问题我再叠了一层粗糙金字塔对每个1°格网预计算一个低分辨率概览高程值用于快速显示和拒绝那些明显越界的查询避免无谓的磁盘读取。3.3 压缩率实测与启动加载速度具体数字这样记录一下方便你比对自己的方案国内SRTM3裸栅格约2.3GB差值编码压缩后约420MB索引表与金字塔概览约28MB瓦片索引加载耗时冷启动中端安卓机约320ms单点高程查询平均耗时3~8ms如果改用ALOS 30m做国内精扫层原始约18GB压缩后约1.8GB下载包做成省份为单位单省基本在80~300MB之间手机上完全可接受。3.4 双线性插值查到的海拔不是格网角的死值还有一个容易忽略的细节查任意坐标时如果直接取离它最近的格网角上的高程在陡峭山区容易出现一两跳的台阶感比如沿着山脊走海拔跳变非常生硬。我最终在查询层统一使用了双线性插值取目标坐标周围四个格网点按距离加权计算最终高程值。公式很简单E w00*e00 w10*e10 w01*e01 w11*e11权重就是目标点在四个格网包围盒内的归一化距离。实测下来双线性插值让轨迹剖面平滑很多而且在90m分辨率下没有引入额外噪声。值得一提的坑是插值前必须先判断四个角点是否有无效值如果有一个是哨兵值就退化成最近邻或用剩余有效角点平均不然会把-32768直接插出一条断崖。4. 坐标系的坑远比想象的大WGS84、GCJ-02、BD-09 与高程的连带关系4.1 三种坐标系的关系与来源高程数据本身是附着在经纬度网格上的经纬度属于哪个坐标系直接决定你拿着坐标能不能查到正确的位置点。这里必须讲清楚三套坐标WGS84全球定位系统的原始坐标系GPS芯片直接输出的经纬度就是WGS84谷歌地球、海外版地图、以及绝大多数国际地图用的都是它。GCJ-02俗称“火星坐标”是在WGS84基础上加了一套非线性偏移算法的坐标系国内的高德、腾讯地图和谷歌中国版都使用它偏移量在城市里通常几十到几百米。BD-09百度在GCJ-02基础上又叠加了一次偏移只被百度地图使用偏移量比GCJ-02再大一些。这三个坐标系的差别在城区地图上看着只是几百米但在山区查海拔时会被成倍放大因为山地地形梯度大水平偏移几百米查到的海拔可能相差三五十米甚至从山腰直接查到山底。4.2 输经纬度查海拔时必须让用户确认坐标系这是我在做“输入经纬度查询”时踩得最深的一个坑。最初版本默认输入的都是WGS84坐标结果很多用户从高德地图上复制了一对GCJ-02坐标过来查询结果跟实际地点对不上而且在山里误差尤其明显。后来我狠下心在输入界面增加了坐标系选择入口WGS84、GCJ-02、BD-09三个选项默认WGS84但在用户粘贴坐标时给出“如果你是从高德复制的请选GCJ-02”的提示。查询内部流程是无论用户输入哪个坐标系先把坐标统一转换到WGS84再拿WGS84坐标去索引DEM瓦片。对GPS实时模式也一样GPS芯片输出的就是WGS84不需要额外转换但如果用户在设置里打开了“显示高德坐标”展示层单独做一次正向转换即可。4.3 从GPS实时坐标到DEM索引的转换链路GPS实时模式下每次定位回调到的是一个带时间戳的WGS84坐标我直接用这个坐标去查DEM瓦片然后把DEM返回的正高海拔显示在界面上。为了节约电量定位频率做了动态降级当设备静止时把GPS更新频率降到10秒一次只有检测到位移大于10米时才切到1秒高频。这个调整对实时海拔显示的体验几乎没有影响省电效果却很明显。4.4 WGS84与GCJ-02互转的参考实现网上关于“python将gps经纬度转换为高德经纬度”的搜索热度一直很高说明很多人被坐标转换卡住过。这里给出一份我实测过的标准WGS84转GCJ-02实现国内区域有效import math a 6378245.0 ee 0.00669342162296594323 def out_of_china(lng, lat): return not (73.66 lng 135.05 and 3.86 lat 53.55) def _transform_lat(x, y): ret -100.0 2.0 * x 3.0 * y 0.2 * y * y ret 0.1 * x * y 0.2 * math.sqrt(abs(x)) ret (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 ret (20.0 * math.sin(y * math.pi) 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0 ret (160.0 * math.sin(y / 12.0 * math.pi) 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0 return ret def _transform_lng(x, y): ret 300.0 x 2.0 * y 0.1 * x * x ret 0.1 * x * y 0.1 * math.sqrt(abs(x)) ret (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 ret (20.0 * math.sin(x * math.pi) 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0 ret (150.0 * math.sin(x / 12.0 * math.pi) 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0 return ret def wgs84_to_gcj02(lng, lat): if out_of_china(lng, lat): return lng, lat dlat _transform_lat(lng - 105.0, lat - 35.0) dlng _transform_lng(lng - 105.0, lat - 35.0) radlat lat / 180.0 * math.pi magic math.sin(radlat) magic 1 - ee * magic * magic sqrt_magic math.sqrt(magic) dlat (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi) dlng (dlng * 180.0) / (a / sqrt_magic * math.cos(radlat) * math.pi) return lng dlng, lat dlat反向GCJ-02转WGS84没有完全精确的解析式工程上常用迭代法先正向转换得到偏移后的坐标再不断修正差值逼近。误差能收敛到1米以内在这个工具的精度需求下完全够用。5. 离线地址搜索与反向地名没有服务端怎么做“输入地址查询”5.1 地址库数据源的选择与手工补全最麻烦的功能其实不是经纬度查询而是“输入地址查询海拔”。在线方案里这是地理编码服务的活离线状态下一个地名要变成经纬度只能靠本地预置地址库。我做了一个务实的取舍地址库不追求搜索引擎级别的全量而是覆盖“户外和日常查询真正会用到”的地点。主体来自公开整理的全国县级以上行政区划点约2800多条加上4A、5A级景区约2500个再加上手工整理的山峰、垭口、镇村、河流桥渡、徒步路线起终点等约15000条。地址库最终做成一条JSON数据包含名称、所属省份、类型、拼音、别名、WGS84坐标。体积不到5MB加载进内存的耗时可以忽略。这里比较花时间的是手工补全部分。很多户外用户会直接搜“牛背山”“党岭”“子梅垭口”“武功山金顶”这类非行政区名称公开数据集里根本没有。我是从各种户外路书里人工提取并核对坐标后一条条补进去的这个脏活没捷径但对查询体验的提升非常直接。5.2 分词、拼音前缀与模糊匹配的轻量方案地址搜索模块我本来想引入完整分词库后来发现杀鸡用牛刀。因为预置地址库是点数据不是长文本真正需要匹配的是名称片段、前缀和拼音。我最终实现了一个三层匹配第一层精确匹配完整名称命中直接返回。第二层前缀匹配比如输入“珠穆”能命中“珠穆朗玛峰”。第三层模糊匹配用Levenshtein距离控制在2个字符以内处理“布达拉”“布拉达”这类错别字。同时支持拼音匹配把每个地点的全拼和拼音首字母预计算好输入“zmlmf”或“zml”都能搜到珠峰。整个索引结构用Trie树实现启动时构建一次搜索耗时都在1毫秒级别。5.3 经纬度反查地名的空间索引实现反向地名输入坐标告诉用户这附近叫什么同样重要尤其GPS实时模式下显示“当前位置海拔4283米附近约2.3km为子梅垭口”比单纯一个数字直观得多。反向查询我用的是两级空间索引先建一个经纬度网格每个0.05°格子记录落在其中的地址点ID列表再有就是每个地址点预先计算一个搜索半径。查询时先定位格子再取周边邻接格子里的候选点计算大圆距离返回最近距离和名称。网格划分本身是一劳永逸的任何点查询都是常数时间定位加少量距离计算。5.4 地址搜索的边界情况与兜底策略地址搜索有几个逃不掉的边界情况我做了一组兜底策略同名地点比如“锦江”在成都和南昌都有返回结果时按省份和类型排序并显示省份消歧义。查询结果无匹配不直接报错改用“经纬度反查最近地点”推荐用户是否要查看附近地标。用户输入的是门牌号或长段文字这种离线方案不支持此时提示用户改用经纬度查询或输入附近地标名称。在产品上我始终把“经纬度查询”作为最可靠路径地址查询定位为辅助入口这样用户预期比较清楚不会因为地址库覆盖不全而产生“工具不好用”的负面判断。6. 实时海拔显示的软硬件博弈GPS高度、大地水准面与气压计6.1 为什么GPS报的海拔能差几十米很多人有个误区以为GPS输出的高度就是海拔。实际上GPS芯片直接输出的是相对于WGS84椭球面的椭球高度而日常说的“海拔”是相对于大地水准面近似平均海平面的正高。两者之差叫大地水准面差距N在全球范围有正有负绝对值通常在二三十米到近百米。在中国境内N大约在-9米到96米之间变化直接拿GPS高度当海拔误差几十米很正常。6.2 大地水准面修正与EGM96格网的使用要让DEM查询结果与真实海拔对得上需要把椭球高转成正高。最简单的做法是内嵌一份EGM96大地水准面模型我用的是一张15弧分钟分辨率的全球网格量化后体积约1.3MB。查询时对当前位置做双线性插值得到N然后用公式H h_ell - N即可。不过这里还有一个中国特色问题国内现行高程基准是1985国家高程基准和全球EGM96正高之间存在一个区域性的系统差。为了兼顾精度我在设置里提供了“高程基准微调”选项允许用户基于已知点手动校正一个偏移量。实测在川西和北京周边校正后的海拔与当地测绘桩相比误差能压到±5米以内对户外使用已经完全可用。6.3 DEM海拔、GPS海拔、气压计海拔三源融合工具的核心数据来源是DEM但DEM是静态数据它反映的是“这个坐标点的地表高程”而不是“你脚下的海拔”。这两者在陡坎、建筑物、桥梁场景下会有差异。因此我在实时模式里做了一层信号融合首选DEM海拔以GPS坐标查询DEM得到稳定海拔作为主显示值。次选GPS椭球高修正值当DEM查询失败或处于无效区域时使用EGM96修正后的GPS高度兜底。可选气压计设备支持气压计时允许用户先校零之后用气压变化率推算相对升降用于补充短时间内的剧烈起伏。三种信号的切换逻辑是DEM有效时以DEM为准并做一阶低通滤波EMA系数取0.3避免GPS坐标抖动导致显示值乱跳DEM无效时才落到下游信号。用90m DEM做实时海拔显示单次查询3~8毫秒就算每秒刷新一次电量开销也远低于在线方案。6.4 极简界面的设计与省电调优“体验友好”和“极精简”最终要落到界面上。我的做法是单一主界面顶部一个地址/经纬度输入框中部一个超大字号的海拔显示左下角显示当前坐标系类型右下角一个实时定位按钮。没有任何二级页面地图可选不显示因为很多用户根本不需要看得懂地图他们只想要一个准确的数字。省电方面除了前面说的GPS动态降频还把屏幕常亮改为可选的“徒步模式”用深色主题减少OLED功耗。数据上整套离线包在国产中端机上冷启动到完成首次海拔显示实际约1.2秒后续每次查询体感都是即时响应。对比在线接口动辄一两秒的等待这个“终极版”的体验差异非常明显。7. 精度实测与踩坑记录从平原到垭口的真实数据7.1 三类场景的实测对比数据工具做了好几轮实地校验这里给出一组有代表性的结果测试点描述工具显示参考高程误差青岛某验潮站附近海边平地6m4m2m北京香山鬼见愁山地坡度大558m550m测绘桩8m川西折多山垭口高海拔垭口4302m4298m4m秦岭某河谷峡谷GPS信号差1563m1570m-7m从数据看绝大多数场景误差在±10米内山脊、垭口这类地形起伏明显处反而比较准因为DEM网格本身就反映了地表形态。误差最大的往往是峡谷或陡坡边缘原因是90m格网在水平方向已经跨过了一段垂直高差无论如何插值都会平滑掉局部细节。7.2 SRTM空洞与异常栅格的处理SRTM最坑人的地方是数据空洞和水体异常。比如一些高山湖泊、深谷阴影区、雷达回波失效区在原始数据里是NoData还有些水体区域高程不是0而是负值或假值。如果直接把这些值暴露给用户会出现“在湖中心查到海拔-230米”的笑话。我的处理分三步第一步在离线转换阶段用GDAL的fillnodata做空洞填充第二步对已知水体边界做掩膜将湖泊水库区域统一替换为周边平均高程并打上“水体估算”标签第三步在查询层做合理性校验如果插值结果低于该瓦片最低合理高程或高于最高合理高程自动回退到金字塔概览值。7.3 坐标转换误差、内存溢出与其他坑写代码过程中真正让我头疼的坑有三个坐标转换误差被手工放大早期反向转换GCJ-02用了一次简单近似误差在山区被地形放大成几十米的高程偏差。后来改成5次迭代修正误差压到1米级。整文件读入内存导致OOM初版图省事把400MB压缩包直接用ByteBuffer映射进内存低内存安卓设备直接崩。改成RandomAccessFile按瓦片随机读取后运行内存占用稳定在60MB左右彻底解决。瓦片边界接缝不同1°格网瓦片在压缩时用了独立统计参数边界上同一位置两片瓦片查询结果有轻微跳变。后来我给瓦片做了一层重叠每个瓦片实际多存一圈缓冲格网点接缝问题直接消失。这四个字说给做离线地理数据的同行校验边界。无论是瓦片边界、坐标系边界还是数据有效性边界所有问题几乎都发生在边界上。另外想专门提一下“离线”这个方向本身。现在大家找离线安装包、离线地图、离线大模型部署本质上都是同一个诉求在不可靠的网络环境里把关键能力掌握在自己手里。高程海拔工具只是其中一个很小但很有代表性的切面它用到的数据压缩、索引构建、坐标转换和离线地名库技术完全可以平移到离线等高线生成、离线剖面绘制甚至离线坡度分析上去。如果你正在做类似的东西我个人最想强调的体会是不要在原始数据精度上贪多先把90m分辨率的SRTM吃透、把坐标系和插值处理干净用户能感知到的提升远比盲目上30m数据大得多。等核心链路稳定后再把30m数据做成可选精扫包增量提供这才是“极精简”和“可用精度”之间的真实平衡点。
返回列表