0
点赞
收藏
分享

微信扫一扫

回忆2012-入坑OpenStreetMap10年纪

醉东枫 2022-01-13 阅读 74
OpenStreetMapPostGISLinux

DOOMED

小贝：如果有一天灾难降临地球，人类文明崩塌，活到后天的幸存者最需要准备什么？
小瓦：指南针，低功耗太阳能 LCD 手持计算机，灾难前最后一周的 OpenStreetMap 镜像，短波电台，USRP B210，我想想，还有摩尔斯代码表。
小贝：难道我们不需要去弄些罐头、卫生纸、盘尼西林吗？
小瓦：当然，找吃的这样光荣的任务就交给你了！

2012年，电影《2012》上映，引发灾难电影的热潮。《后天》、《天地浩劫》风靡全球，玛雅神话，气候变迁，都成为茶余饭后的话题。某天，笔者突然意识到OpenStreetMap应该是应对灾后世界的必备法宝之一。于是，我的第一个OpenStreetMap渲染服务器虚拟机镜像搭建完成，同时完成了部分Qt/MFC客户端代码。
2012-2022，几乎每年都会推出一版虚拟机镜像，并在网盘开源。
2012-2016，用Qt开发客户端，完成了从PC向低功耗树莓派的迁移。
2018-2021，陆续加入数字高程DEM图层、海洋深度图层。
2022年1月，基于Manjaro的最新VMBox服务器虚拟机制作完成。
在此十年之际，聊聊玩OpenStreetMap的开心事，以此纪念逝去的时光。虚拟机镜像可以从我的网盘下载。

1.OSM安全与法律问题提醒

为了大家安全开心的玩OSM，这里重要的话写在前头：

** 重要 ** ：在很多国家包括我国，由于地理信息数据关系到冠名、划界等问题，对公众发布需要严格的审批流程，并承担严苛的法律责任。没有许可证，不能向公众开放地理信息服务！

** 重要 ** ：下载OSM数据是合法的，但请不要随便采集、上传国内地标。上传地标可构成非法测绘，若上传敏感地标可能官司缠身。

** 重要 ** ：OSM的数据精度、地标的准确性是没有严格保证的。使用这些数据进行导航、户外穿越需要额外注意。

举例，测试网站www.goldenhawking.site:8088 ，若不备案，需要按照相关要求，需要封堵80端口，仅在非标准大端口提供到OSM官网的链接。不能使用动态网页，不能使用含有用户注册、留言功能。也就是只能发行非80端口的静态的网页。

2.OSM知识点走马观花

除了作为末日危机、全球断网前必须收集的重要数据财富，OpenStreetMap也有很多值得书写的知识点、亮点。每年更新一版全球镜像，可玩性很高。

  • PostgreSQL灵魂应用：PostGIS，PostgreSQL，hstore， pgRouting
  • Linux下二进制数据操作：changeset、osmsis osm2pgsql
  • 渲染服务器搭建：apache2 mod_tile renderd fcgi
  • Web前端：CesiumJS、OpenLayers
  • 后端与桌面开发：Qt、QGIS、C#
  • 数据处理、矢量服务 QGIS MapBox GeoServer
  • 运维：虚拟化、装机、SSD阵列、手持全球服务器
  • Linux：Debain、Fedora，CentOS，Ubuntu、Mint、ArchLinux、Manjaro

10年来，不停的寻找最香的发行版，购买SDD，配置强大的8TB SSD盘阵，索引表空间物理分离，都可以拿出来写。当然，也要感谢家人的支持，一会玩OSM，一会SDR，陆续花了好几万，二手服务器、工作站，还有居高不下的电费。这还不算树莓派这些零碎。教研室微薄的奖金全用来折腾了。

这些知识点太零散，用业余时间是搞不深的，蜻蜓点水可以看我的专栏了解其中的少数，或者直接从官网维基系统学习。今天说一些轻松的小知识点。

注意，下面蜻蜓点水的介绍，不足以搭建一个瓦片服务器。真实的流程比较复杂。

3 获取和导入数据

最新的PBF格式数据超过60GB，且下载速度不快（也可能和我的校园网和别人分享有关） 。建议找一个树莓派外挂移动硬盘不间断下载。

3.1 硬件条件

目前，导入数据时的中间结果会超过1TB。因此，建议如下搭配：

  1. 物理主机直接导入，2TB SSD ，48GB内存
  2. 物理机+虚拟机： 2x2TB SSD （4TB），物理机64GB内存，虚拟机48G

CPU的关系不是很大，笔者使用了 i7 8700 CPU。

3.2 导入数据

在4TB SSD上，使用osm2pgsql导入数据耗费36小时。在国内网路条件下，使用changeset保持更新是不现实的。我们一般可采用每季度、每年更新一次的方法。因此，使用flatnode可以提高速率。

$ osm2pgsql -c -s -S"./openstreetmap-carto.style" -C28000 -dgis --drop --hstore --flat-nodes "./cache/flat_node"  ”./downloads/planet-latest.osm.pbf"
2022-01-01 14:29:00  osm2pgsql version 1.5.1
2022-01-01 14:29:00  Database version: 13.4
2022-01-01 14:29:00  PostGIS version: 3.0
2022-01-01 14:29:01  Setting up table 'planet_osm_point'
2022-01-01 14:29:01  Setting up table 'planet_osm_line'
2022-01-01 14:29:01  Setting up table 'planet_osm_polygon'
2022-01-01 14:29:01  Setting up table 'planet_osm_roads'
Processing: Node(7402483k 2028.1k/s) Way(825850k 28.61k/s) Relation(9527340 412.11/s)                                                                               
2022-01-02 05:56:17  Reading input files done in 55636s (15h 27m 16s).
2022-01-02 05:56:17    Processed 7402483451 nodes in 3650s (1h 0m 50s) - 2028k/s
2022-01-02 05:56:17    Processed 825850076 ways in 28869s (8h 1m 9s) - 29k/s
2022-01-02 05:56:17    Processed 9527913 relations in 23117s (6h 25m 17s) - 412/s
2022-01-02 05:56:19  Dropping table 'planet_osm_nodes'
2022-01-02 05:56:19  Done postprocessing on table 'planet_osm_nodes' in 0s
2022-01-02 05:56:19  Dropping table 'planet_osm_ways'
2022-01-02 05:56:20  Done postprocessing on table 'planet_osm_ways' in 1s
2022-01-02 05:56:20  Dropping table 'planet_osm_rels'
2022-01-02 05:56:20  Done postprocessing on table 'planet_osm_rels' in 0s
2022-01-02 05:56:20  Done postprocessing on table 'planet_osm_nodes' in 0s
2022-01-02 05:56:20  Done postprocessing on table 'planet_osm_ways' in 0s
2022-01-02 05:56:20  Done postprocessing on table 'planet_osm_rels' in 0s
2022-01-02 05:56:20  Clustering table 'planet_osm_roads' by geometry...
2022-01-02 05:56:20  Clustering table 'planet_osm_polygon' by geometry...
2022-01-02 05:56:20  Clustering table 'planet_osm_point' by geometry...
2022-01-02 05:56:20  Clustering table 'planet_osm_line' by geometry...
2022-01-02 07:22:58  Creating geometry index on table 'planet_osm_roads'...
2022-01-02 07:47:18  Analyzing table 'planet_osm_roads'...
2022-01-02 07:57:26  Creating geometry index on table 'planet_osm_point'...
2022-01-02 10:42:03  Analyzing table 'planet_osm_point'...
2022-01-02 10:43:41  All postprocessing on table 'planet_osm_point' done in 17241s (4h 47m 21s).
2022-01-02 15:59:45  Creating geometry index on table 'planet_osm_line'...
2022-01-02 17:32:06  Creating geometry index on table 'planet_osm_polygon'...
2022-01-02 20:08:14  Analyzing table 'planet_osm_line'...
2022-01-02 20:08:38  All postprocessing on table 'planet_osm_line' done in 51137s (14h 12m 17s).
2022-01-03 02:36:37  Analyzing table 'planet_osm_polygon'...
2022-01-03 02:37:31  All postprocessing on table 'planet_osm_polygon' done in 74471s (20h 41m 11s).
2022-01-03 02:37:31  All postprocessing on table 'planet_osm_roads' done in 6819s (1h 53m 39s).
2022-01-03 02:37:41  osm2pgsql took 130121s (36h 8m 41s) overall.

在完成了必要的汉化、数据纠错和处理后，更新整体数据时间戳

sudo touch  /var/lib/mod_tile/planet-import-complete

4 渲染瓦片

渲染瓦片是非常耗时的工作。我们需要一个专用的工作站进行最初的数据渲染工作。一旦完成了常见区域10级以上的渲染，后续按需绘制会快得多。
OSMServer

4.1 建立额外的索引解决renderd速度过慢的问题

额外的索引让renderd更快速, 在OSM的openstreetmap-carto 工程里，有一个脚本叫做index.py，运行后会生成SQL：

-- These are indexes for rendering performance with OpenStreetMap Carto.
-- This file is generated with scripts/indexes.py
CREATE INDEX planet_osm_line_ferry  ON planet_osm_line USING GIST (way)  WHERE route = 'ferry' AND osm_id > 0;
CREATE INDEX planet_osm_line_label  ON planet_osm_line USING GIST (way)  WHERE name IS NOT NULL OR ref IS NOT NULL;
CREATE INDEX planet_osm_line_river  ON planet_osm_line USING GIST (way)  WHERE waterway = 'river';
CREATE INDEX planet_osm_line_waterway  ON planet_osm_line USING GIST (way)  WHERE waterway IN ('river', 'canal', 'stream', 'drain', 'ditch');
CREATE INDEX planet_osm_point_place  ON planet_osm_point USING GIST (way)  WHERE place IS NOT NULL AND name IS NOT NULL;
CREATE INDEX planet_osm_polygon_admin  ON planet_osm_polygon USING GIST (ST_PointOnSurface(way))  WHERE name IS NOT NULL AND boundary = 'administrative' AND admin_level IN ('0', '1', '2', '3', '4');
CREATE INDEX planet_osm_polygon_military  ON planet_osm_polygon USING GIST (way)  WHERE (landuse = 'military' OR military = 'danger_area') AND building IS NULL;
CREATE INDEX planet_osm_polygon_name  ON planet_osm_polygon USING GIST (ST_PointOnSurface(way))  WHERE name IS NOT NULL;
CREATE INDEX planet_osm_polygon_nobuilding  ON planet_osm_polygon USING GIST (way)  WHERE building IS NULL;
CREATE INDEX planet_osm_polygon_water  ON planet_osm_polygon USING GIST (way)  WHERE waterway IN ('dock', 'riverbank', 'canal')    OR landuse IN ('reservoir', 'basin')    OR "natural" IN ('water', 'glacier');CREATE INDEX planet_osm_polygon_way_area_z10  ON planet_osm_polygon USING GIST (way)  WHERE way_area > 23300;
CREATE INDEX planet_osm_polygon_way_area_z6  ON planet_osm_polygon USING GIST (way)  WHERE way_area > 5980000;
CREATE INDEX planet_osm_roads_admin  ON planet_osm_roads USING GIST (way)  WHERE boundary = 'administrative';
CREATE INDEX planet_osm_roads_admin_low  ON planet_osm_roads USING GIST (way)  WHERE boundary = 'administrative' AND admin_level IN ('0', '1', '2', '3', '4');
CREATE INDEX planet_osm_roads_roads_ref  ON planet_osm_roads USING GIST (way)  WHERE highway IS NOT NULL AND ref IS NOT NULL;

** 如果没有上述索引，renderd在13级别比例尺一下（14，15…）会灰常灰常慢！建立索引能够提高渲染速度1000倍以上。 **

4.2 新服务器-预先渲染前10层

如果是崭新的服务器，没有瓦片，则建议预先渲染前10层，注意，在不同环境中，花销不同。内存、磁盘、CPU有一个瘸腿，就可能花费超过72小时。

render_list -n 4 -a -z0 -Z9 -t /var/lib/mod_tile

预先渲染10层后，再按需渲染更深的层次即可。

4.3 旧服务器-更新存在的旧瓦片

如果是旧的服务器，瓦片已经存在，则执行

render_old -n 4 -t /var/lib/mod_tile

可以更新所有瓦片

4.4 技巧：用Morton编码显著降低磁盘读写

如果需要从特定的txt文件中渲染更深层次的历史瓦片记录，则需要准备一个txt文件，每行一个瓦片，坐标 x y z。历史记录是来自老的服务器上的文件夹记录、数据库记录。

36 88 9
...
0 0 0
...
23 24 8

但这些记录很可能是随机的、零散的。渲染这些零散的瓦片，会反复读取磁盘，效率很差。我们使用一段程序，读取零散记录，用Morton编码排序后输出：

#include <cstdio>
#include <map>
#include <array>
#include <cassert>
#include <list>
int main()
{
	/*tilehis 存储历史记录，每行是空格分割的 x y z*/
	FILE * fpIn = fopen ("tilehis.txt","r");
	if (!fpIn)
		return -1;
	char buf[1024];
	std::map<unsigned long long, std::list<std::array<unsigned long long,3> > > mx;
	printf("Reading...\n");
	while (fgets(buf,1024,fpIn))
	{
		unsigned long long x,y,z;
		sscanf(buf,"%llu%llu%llu",&x,&y,&z);
		assert (z>=0 && z<=22);	
		unsigned long long xy = z<<56;
		for (unsigned int n = 0; n<32;++n)
		{
			unsigned long long bx = (x >> n) & 0x01 ;
			unsigned long long by = (y >> n) & 0x01 ;
			if (bx)
				xy |= (1llu<<(2u*n+1u));
			if (by)
				xy |= (1llu<<(2u*n));
		}
		std::array<unsigned long long,3> a;
		a[0] = x;
		a[1] = y;
		a[2] = z;
		mx[xy].push_back(a);
	}
	fclose(fpIn);	
	for (auto p = mx.begin();p!=mx.end();++p)
	{
		for (auto q = p->second.begin();q != p->second.end();++q)
			printf ("%llu %llu %llu\n",(*q)[0],(*q)[1],(*q)[2]);	
	}
	printf ("Finished.\n");
	return 0;
}

输出后的坐标，相邻二维地图的区域，在一维的文件记录里也相邻

...
25778 13410 15
25800 13418 15
25800 13419 15
25801 13418 15
25801 13419 15
25802 13418 15
25802 13419 15
...
25803 13422 15
25804 13418 15
25804 13419 15
25805 13418 15
25805 13419 15
25806 13418 15
25806 13419 15
...

如此处理，64GB内存能够缓存相当多的临近数据，缓存的命中率超级高。外加索引加持，很多时间硬盘都不闪，只有CPU 100%

注意不要换算不同z层比例尺的地图到一块，那样反而过犹不及。因为不同层次需要的要素是不同的，混在一起，9级比例尺需要河流，桥梁，17级需要商店，结果导致系统缓存装不下这莫多的数据。

5 OpenStreetMap的意义

目前，GIS的现状是一线大厂分天下，二线跟着吃周边。OpenStreetMap也是具备导航模块的，搭建
pgRoutings就可以按照地理距离或者时间来导航。然而，正如大部分开源工具一样，使用OSM也总是给人 “貌似能行，可还是差了一点” 的感觉。作为商业应用，国内的用某度、某德，跨国的用BING都是很好的选择。OpenStreetMap的工具链还是过于零碎，数据的质量也参差不齐。

OSM的意义，是允许独立的个人拥有较为完整的全球地理信息数据，而不是少量缓存在电脑里的瓦片图片。

  • 拥有全部的数据，就可以对很多元素进行统计，得到有意义的一些信息。比如我的专栏中统计城市发达程度。
  • 在脱离互联网的情况下，进行GIS分析。对于学校上机机房这种不许学生上互联网的地方，使用瓦片服务获得全球体验。
  • 迷途求生中需要准备的神器。想象一下如果互联网不再连接世界，我们重新回到信息的孤岛，手中有一份全球地图的意义就非常大了。

OpenStreepMap

举报

相关推荐

带你入逆向坑，怎样在win10上安装并使用Z3库

v8v9方程组HarmonyOS后端开发

测试人生 | 毕业2年，拒绝独角兽入职名企大厂涨薪10万+，这个95后小姐姐好飒

测试测试开发功能测试自动化测试虚拟化云计算

为什么都说测试岗是巨坑，趁早跳出去？10年测试人告诉你千万别上当了...

网络安全运维
0 条评论