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;

结果如下:

DistanceDirection
5650.21917119331542.709389957366675

用人类的话说,这意味着 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)
)

结果如下:

nameadmin_level
Oregon4
Multnomah County6
Portland8
Portland Downtown10

这里需要解释一下,我们使用的数据集仅限于俄勒冈州,因此它不包含国家/地区等对象。 因此,你不会在结果中看到美国。 其他 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翻译整理,转载请标明出处