PostGIS + OpenStreeMap地理空间分析实战
我一直是一个地图极客，可以追溯到 80 年代，当时我会拿一张道路地图集和一些描图纸，然后在我自己的道路网络中绘制。 我最喜欢的游戏之一是拿一张旧地图或地球仪，并尝试根据国家的名称和形状来确定它的制作时间。 在在线地图软件、大数据和数据科学时代，它变得更加有趣。 现在我可以下载地图的整个数据集并针对它编写程序了！推荐：用 NSDT设计器 快速搭建可编程3D场景。
1、Project Remote
最近，我一直在关注 Project Remote，这是一个家庭冒险计算和记录美国每个州离道路最远的点。 在撰写本文时，他们已经完成了 33 个州，而俄勒冈州不是其中之一。 这让我开始思考。 Project Remote 没有公布他们的确切公式，但他们给出了一些提示。 基于这些，我可以计算出俄勒冈州最偏远的点吗？
好吧，这有点复杂。 俄勒冈州有成千上万条道路，我如何计算那么多数字？ 也许我会从小一点开始。
我将只使用高速公路而不是所有道路，因为高速公路并不多。 我将从波特兰开始，因为那是大多数高速公路所在的地方，也是人们熟悉的地方。
似乎有两种方法可以找到距离高速公路最远的点：
[*]使用蛮力方法，设置点网格并计算每个点与最近的高速公路的距离
[*]编写一个 AI 从任意点开始并远离高速公路
对于 Metal Toad 2017 年 12 月的数据科学黑客马拉松，我和我的同事 Kalina Wilson 认为 AI 选项会更有趣。
那么，怎么做呢？ 首先，我有以下问题：
[*]给定一组纬度/经度坐标，计算机能否确定该点是否在俄勒冈州内？
[*]计算机可以测量到最近的高速公路的距离吗？
[*]计算机能否将坐标移动到离高速公路较远的附近点？
事实证明，所有这些问题的答案都是“是”，而且比我预期的要容易。
2、数据和工具
首先，我们需要一些数据和一些工具来分析它。 我选择了以下工具，因为它们都是免费和开源的，而且它们都可以很好地协同工作。
2.1 地图数据 - OpenStreetMap
OpenStreetMap 是开源地图软件。 其地图数据完全由社区构建（即众包），可免费下载。 数据集包含州、县、市的形状信息； 道路和高速公路； 江河湖海； 国家公园和森林； 以及企业、邮局和山峰等兴趣点。
全球 OpenStreetMap 数据集压缩为 40 GB。 这太啰嗦了，所以各种组织都提取了特定的区域，比如一个国家或州。 对于这个项目，我们使用了俄勒冈州的数据，即 134 MB。 你可以在 OpenStreetMap wiki 上阅读所有不同的下载选项。
2.2 GIS数据库 - PostGIS
PostGIS 是 PostgreSQL 的扩展，用于处理几何和地理数据。 它提供了一个 Geometry 列类型和几十个与几何相关的函数。 可用函数可以执行如下计算：
[*]找出两个形状之间的距离
[*]找到从一个形状或点到另一个的方向
[*]判断一个点是否在多边形内
[*]单位转换，例如，将纬度/经度转换为米/英尺
OpenStreetMap 数据采用称为 PBF 的特定格式，我们可以使用称为 osm2pgsql 的开源工具将其加载到 PostgreSQL 中。
2.3 地理信息系统 - QGIS
QGIS 是一个用于查看地理信息系统数据的桌面应用程序。它有一个可以直接连接到 PostgreSQL 的数据库管理器。 你可以使用 PostGIS 函数和常规 SQL“where”子句编写 SQL 查询，并将结果作为图层包含在地图上。
2.4 地理编码处理库 - Geopy
Geopy 是一个 Python 库，用于处理地理编码数据并执行几何和地理计算。 它可以通过调用各种网络服务来查找位置和地址。 它还可以执行计算以测量两点之间的距离。 使用一个几乎没有记录的功能，我了解到它还可以做诸如“从（北纬 45.51°，西经 122.67°），以罗盘航向东北 45° 行驶 500 米”之类的事情。
2.5 前端可视化 - Leaflet|MapBox|Turf
为了在浏览器中可视化，我们使用了：
[*]Leaflet JS ， 这是一个流行的 JavaScript 库，用于创建交互式地图。
[*]MapBox GL JS ，最初基于 Leaflet JS，但使用矢量切片，这增加了重要的客户端功能（缩放、元数据查询、样式）。
[*]Turf JS ， 添加了空间分析和统计工具，可以进行测量、计算和选择。
3、处理数据
将俄勒冈数据集加载到数据库后，我可以使用数据库前端工具查看数据。 我发现我现在有了三个有趣的表：点、线和多边形。 每个代表一种几何形状并包含具有该几何形状的每个对象。 “点”包含商户、山峰等点源要素。 “线”包括道路、河流、湖泊和其他“开放”的几何形状。 “多边形”表示闭合形状，如州、县、城市、国家森林和湖泊。
每个表都包含许多包含分类信息（道路类型、水域类型、边界类型、海拔等）的文本列和另一列称为“way”的列。 “way”是 OpenStreetMap 所说的地理形状。 它表示为 PostGIS“几何”数据类型。 在数据库中，Way 看起来像 0101000020E6100000FA7E6ABC74AB5EC0DF4F8D976E424740，它是地理形状“POINT(-122.679 46.519)”的编码表示。
PostGIS 提供了许多用于在 SQL 中操作方式和形状的函数。
示例：查找从波特兰先锋法院广场到时尚中心的方向和距离：
SELECTST_Distance( ST_Transform(ST_GeomFromText('POINT(-122.679 45.519)', 4326), 2269), ST_Transform(ST_GeomFromText('POINT(-122.667 45.532)', 4326), 2269)) AS distance,degrees(ST_Azimuth( ST_GeomFromText('POINT(-122.679 45.519)', 4326), ST_GeomFromText('POINT(-122.667 45.532)', 4326) )) AS direction;结果如下：
Distance
Direction
5650.219171193315
42.709389957366675
用人类的话说，这意味着 5650 英尺，航向为东北 42.7°。
再例如，我们可以使用函数 ST_Intersects() 来查找某个点所在的城市、县和州的名称。 请记住，市、县和州边界存储在数据库的“polygon”表中。 它们由“boundary”和“admin_level”字段标识，其中“admin_level”为 4 是州，6 是县，8 是市，10 是社区。
SELECT nameFROM planet_osm_polygonWHERE boundary = 'administrative' AND ST_Intersects( ST_GeomFromText('POINT(-122.679 45.519)', 4326), ST_Transform(way, 4326))结果如下：
name
admin_level
Oregon
4
Multnomah County
6
Portland
8
Portland Downtown
10
这里需要解释一下，我们使用的数据集仅限于俄勒冈州，因此它不包含国家/地区等对象。 因此，你不会在结果中看到美国。 其他 OSM 数据导出可能包含更多类型的边界。
请注意上面 ST_Transform() 函数的使用，以及一些神秘数字，如 4326 和 2269。这就是单位转换在 PostGIS 中的工作方式。 这些数字称为空间参考 ID 或 SRID。 SRID 4326 是纬度和经度。 SRID 2269 是一个特定于区域设置的坐标系，称为“NAD83 / Oregon North (ft)”。 由于地球的形状有些不规则，而且当你接近两极时经线会汇聚，因此每个地方都有自己的坐标系。 为了找到距离，我们需要经常将形状从一个坐标系转换到另一个坐标系。
4、寻找距离道路最远的点
现在我们已经掌握了数据的外观，可以开始工作了。 ST_Distance()、ST_ClosestPoint() 和 ST_Azimuth() 函数适用于任何类型的形状，而不仅仅是点。 因此，例如，我们可以计算直线上任意一点到最近点的距离。 这将是我们“最远一条路”计算的核心。 要找到最近的高速公路及其方向，我们可以运行如下查询：
SELECT osm_id, name, highway, REF, ST_AsText(ST_Transform(ST_ClosestPoint(way, 'SRID=900913;POINT(-122.679 45.519)'::geometry), 4326)) AS closest_point,ST_Distance( way_local, ST_Transform(ST_GeomFromText('POINT(-122.679 45.519)',4326), 2285) ) AS distance, degrees(ST_Azimuth(ST_GeomFromText('POINT(-122.679 45.519)', 4326),ST_ClosestPoint(ST_Transform(way, 4326), ST_GeomFromText('POINT(-122.679 45.519)', 4326)) )) AS azimuthFROM roadsWHEREhighway IN ('motorway')ORDER BY distanceLIMIT 1现在我们知道与我们的点相关的最近道路在哪里，我们只需要对 AI 进行编程。
对于这个项目，我们编写了几个非常简单的 AI，它们都遵循类似的计划：
[*]查找到最近道路的距离和方向
[*]将点沿相反方向移动一小段距离（例如 1000 英尺）
[*]检查新点是否仍在俄勒冈州内
[*]重新计算新点的距离和方向
[*]重复这个过程几次迭代
同样，它不是太复杂。 如果你将起点放在高速公路环路内，例如环绕波特兰市中心的 I-5/I-405 环路，它甚至不会聪明到知道它位于封闭形状内。 因此，该算法也无法“逃脱”封闭形状以找到其外部的点。 如果离高速公路最远的点实际上在高速公路的另一边，则该算法将找不到它。 它所能做的就是找到局部最大距离处的点。 解决这个问题将成为未来AI研究的主题。
5、结束语
一旦我们完成了简单的 AI，我就有了一个更大的目标。 我想看看我是否能完成我最初的想法，即计算俄勒冈州最偏远的点。 而且，为了好玩，华盛顿也一样。
我的计划是在全州创建一个点网格，间隔大约 2 英里。 然后我会计算从每个点到任何类型的最近道路的距离。 几何数运算的计算量很大，因此需要几个小时才能计算出来。
一旦我确定了最偏远点的几个候选点，我就可以使用上述 AI 的修改版本。 我将其编程为远离数据库中的任何道路，甚至是未命名的四轮驱动轨道。
结果证明这是相当成功的。 我在俄勒冈州东北部（位于美丽的瓦洛厄山脉）确定了距离最近道路 7.5 英里的一个点，以及华盛顿奥林匹克国家公园距离最近道路 13.6 英里的一个点。 （我没有公布确切的位置。我不想冒险让这些原始地方被游客淹没。）
为了验证我的发现，我在 Project Remote 网站上找到了一张地图，上面用红点表示他们在每个州确定的远程点。 他们的图像非常缩小，故意不显示确切位置，但他们在俄勒冈州的东北角和华盛顿州的西北角确实有一个红点，在我确定的点附近。
我希望你和我一样觉得这很有趣。
<hr>原文链接：http://www.bimant.com/blog/spatial-analysis-with-postgresql-n-openstreetmap/ 转发了 转发了 转发了 转发了 转发了 转发了 转发了 转发了 转发了 转发了 转发了 转发了 接了个二开的单，前端Vue，后端PHP。本想，这不手到擒来！看了看后端源码，有点懵，仔细研究加百度搜索，原来用了个“Niucloud-admin”框架，还真没接触过，又涨知识了！
页:
[1]