PostGIS空间分析实战

我一直是一个地图极客，可以追溯到 80 年代，当时我会拿一张道路地图集和一些描图纸，然后在我自己的道路网络中绘制。 我最喜欢的游戏之一是拿一张旧地图或地球仪，并尝试根据国家的名称和形状来确定它的制作时间。 在在线地图软件、大数据和数据科学时代，它变得更加有趣。 现在我可以下载地图的整个数据集并针对它编写程序了！

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 中操作方式和形状的函数。

示例：查找从波特兰先锋法院广场到时尚中心的方向和距离：

SELECT
  ST_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;

结果如下：

用人类的话说，这意味着 5650 英尺，航向为东北 42.7°。

再例如，我们可以使用函数  ST_Intersects() 来查找某个点所在的城市、县和州的名称。 请记住，市、县和州边界存储在数据库的“polygon”表中。 它们由“boundary”和“admin_level”字段标识，其中“admin_level”为 4 是州，6 是县，8 是市，10 是社区。

SELECT name
FROM planet_osm_polygon
WHERE
 boundary = 'administrative'
 AND ST_Intersects(
    ST_GeomFromText('POINT(-122.679 45.519)', 4326),
    ST_Transform(way, 4326)
)

结果如下：

这里需要解释一下，我们使用的数据集仅限于俄勒冈州，因此它不包含国家/地区等对象。 因此，你不会在结果中看到美国。 其他 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 azimuth
FROM roads
WHERE
  highway IN ('motorway')
ORDER BY distance
LIMIT 1

现在我们知道与我们的点相关的最近道路在哪里，我们只需要对 AI 进行编程。

对于这个项目，我们编写了几个非常简单的 AI，它们都遵循类似的计划：

• 查找到最近道路的距离和方向
• 将点沿相反方向移动一小段距离（例如 1000 英尺）
• 检查新点是否仍在俄勒冈州内
• 重新计算新点的距离和方向
• 重复这个过程几次迭代

同样，它不是太复杂。 如果你将起点放在高速公路环路内，例如环绕波特兰市中心的 I-5/I-405 环路，它甚至不会聪明到知道它位于封闭形状内。 因此，该算法也无法“逃脱”封闭形状以找到其外部的点。 如果离高速公路最远的点实际上在高速公路的另一边，则该算法将找不到它。 它所能做的就是找到局部最大距离处的点。 解决这个问题将成为未来AI研究的主题。

5、结束语

一旦我们完成了简单的 AI，我就有了一个更大的目标。 我想看看我是否能完成我最初的想法，即计算俄勒冈州最偏远的点。 而且，为了好玩，华盛顿也一样。

我的计划是在全州创建一个点网格，间隔大约 2 英里。 然后我会计算从每个点到任何类型的最近道路的距离。 几何数运算的计算量很大，因此需要几个小时才能计算出来。

一旦我确定了最偏远点的几个候选点，我就可以使用上述 AI 的修改版本。 我将其编程为远离数据库中的任何道路，甚至是未命名的四轮驱动轨道。

结果证明这是相当成功的。 我在俄勒冈州东北部（位于美丽的瓦洛厄山脉）确定了距离最近道路 7.5 英里的一个点，以及华盛顿奥林匹克国家公园距离最近道路 13.6 英里的一个点。 （我没有公布确切的位置。我不想冒险让这些原始地方被游客淹没。）

为了验证我的发现，我在 Project Remote 网站上找到了一张地图，上面用红点表示他们在每个州确定的远程点。 他们的图像非常缩小，故意不显示确切位置，但他们在俄勒冈州的东北角和华盛顿州的西北角确实有一个红点，在我确定的点附近。

我希望你和我一样觉得这很有趣。

原文链接：Spatial analysis with PostgreSQL and OpenStreetMap

BimAnt翻译整理，转载请标明出处