使用 SpatiaLite 构建位置到时区 API

SQLite 的 SpatiaLite 扩展增加了大量用于地理空间分析的函数,可以与 Datasette 结合使用来构建 GIS(地理信息系统)应用程序。

本教程将展示如何结合使用 SpatiaLite 和 Datasette 来创建一个 JSON API,该 API 可以返回地球上特定经度和纬度点的时区。

我们将构建什么

您可以在此处试用该 API。提供经度和纬度,它将返回相应的时区 ID:https://timezones.datasette.io/timezones/by_point

一些示例

设置开发环境

您需要为本教程准备好两件事

推荐:使用 GitHub Codespaces

GitHub Codespaces 可以为您提供一个免费的开发环境,通过您的网络浏览器即可访问,并且预装了您所需的所有工具。

这是完成本教程最简单的方法。

访问此链接 创建一个新的 Codespace,其中包含了本教程其余部分所需的一切。

使用 Mac

在 Mac 上,您可以使用 Homebrew 安装 SpatiaLite 和 Datasette,如下所示

brew install spatialite-tools datasette

如果您使用的是 Mac,可能会发现您的 Python 安装无法加载外部 SQLite 模块。您可以通过运行以下命令进行检查

datasette --load-extension spatialite

如果您收到关于 enable_load_extension 的错误消息,请查阅 此页面 获取解决问题的提示。

构建数据库

要构建此项目,我们首先需要全球所有时区的地理形状数据。

timezone-boundary-builder 是 Evan Siroky 的一个项目,它使用来自 OpenStreetMap 的数据创建详细的时区多边形,然后以 GeoJSON 和 Shapefile 格式发布。结果数据根据 开放数据共享开放数据库许可证 (ODbL) 提供。

首先从最新版本下载 timezones-with-oceans.shapefile.zip 文件

wget https://github.com/evansiroky/timezone-boundary-builder/releases/download/2022g/timezones-with-oceans.shapefile.zip

(如果您正在使用 Codespaces,应在 Codespaces 界面的“终端”选项卡中运行所有这些命令。)

一个 Shapefile 是一组描述地理要素集合的文件。无需解压此 zip 文件——我们的工具可以直接使用它。

我们可以使用 shapefile-to-sqlite 命令将文件加载到新的 SpatiaLite 数据库中。

如果您在 Codespaces 中运行本教程,该工具已安装。否则,您可以按如下方式安装它

pip install shapefile-to-sqlite

要将 Shapefile 加载到数据库中,请使用以下命令

shapefile-to-sqlite timezones.db \
    timezones-with-oceans.shapefile.zip \
    --table timezones \
    --spatial-index

这将创建一个名为 timezones.db 的新数据库文件,并将 Shapefile 中的形状加载到一个名为 timezones 的新表中。它还将在 timezones 表上设置一个空间索引,稍后会对此进行描述。

该命令将显示一个进度条,如下所示

zip://timezones-with-oceans.shapefile.zip
  [########----------------------------]   22%  00:00:38

这可能需要几分钟才能完成。

在 Datasette 中浏览数据库

我们现在可以针对该文件启动 Datasette 来浏览数据

datasette timezones.db --load-extension spatialite

您应该会看到以下输出

INFO:     Started server process [5385]
INFO:     Waiting for application startup.
INFO:     Application startup complete.
INFO:     Uvicorn running on http://127.0.0.1:8001 (Press CTRL+C to quit)

如果您使用的是自己的计算机,现在可以访问 http://127.0.0.1:8001 查看 Datasette 界面。

如果您在 Codespaces 中运行,该工具应该会提供一个“在浏览器中打开”按钮。

您现在可以浏览数据库了。应该有一个 timezones 表,包含三列:idtzidgeometry

不过目前没什么可看的!geometry 列只是一些大的二进制数据块。

在地图上查看时区形状

要在地图上查看时区几何图形,我们可以安装一个 Datasette 插件。

Chris Amico 开发的 datasette-geojson-map 插件增加了直接在地图上渲染 GeoJSON 和 SpatiaLite 几何图形的功能。

在终端中按 Ctrl+C 停止 Datasette 服务器,然后运行以下命令

datasette install datasette-geojson-map

现在再次启动 Datasette

datasette timezones.db \
  --load-extension spatialite \
  --setting default_page_size 10

我们在这里添加了一个额外的设置,将默认分页大小设置为 10。这是因为有些时区多边形非常大,默认分页大小 100 可能需要很长时间才能渲染。

再次访问 timezones 表,您应该会看到类似这样的内容

The timezones table page in Datasette shows a map at the top with a small number of time zones represented, above a table listing them by name.

查找点的时区

现在我们有了包含时区几何图形的完整表,我们可以构建一个 SQL 查询来显示地球上任何特定点的时区。

我们可以使用 SpatiaLite 的 within() 函数来实现这一点,该函数接受两个几何图形并检查一个是否包含在另一个之中。

首先,我们需要一个表示特定经纬度点的几何图形。

41.798, -87.696 是芝加哥市内的一个点 - 在 Google 地图上查看

我们可以使用 SpatiaLite 的 MakePoint(longitude, latitude) 函数(请注意这里的经度和纬度是反过来的)来创建一个几何图形,然后可以将其与 within() 函数一起使用

select
  tzid
from
  timezones
where
  within(
    MakePoint(-87.696, 41.798),
    timezones.Geometry
  ) = 1

如果第一个几何图形包含在第二个几何图形中,within(geom1, geom2) 返回 1,否则返回 0

果然,这个查询返回了 America/Chicago

在 SQL 查询中使用参数

Datasette 有一个功能,SQL 查询可以包含命名参数 :like_this,这些参数将被转换为表单字段并用于向查询提供新值。

尝试使用以下查询

select
  tzid
from
  timezones
where
  within(
    MakePoint(:longitude, :latitude),
    timezones.Geometry
  ) = 1

这看起来应该可行...但如果您使用先前的坐标进行尝试,会发现它没有返回任何结果。

这是因为 MakePoint() 需要接收浮点数值,但 Datasette 将所有参数作为字符串传递。

如果您传递了无效类型,MakePoint() 返回 null - 那么 within(null, geometry) 返回 -1 表示结果无效。

我们可以通过将字符串转换为浮点数来解决这个问题,如下所示

select
  tzid
from
  timezones
where
  within(
    MakePoint(cast(:longitude as float), cast(:latitude as float)),
    timezones.Geometry
  ) = 1

这个查询达到了预期的效果:给定地球上任何一个点的经纬度,它将返回正确的时区。

使用索引加速

这个查询还有一个遗留问题:它相对较慢。在我的测试中,查询运行时间在 150ms 到 750ms 之间,这是因为需要将该点与数据库中的所有 455 个多边形进行比较。

我们第一次导入数据时使用了 --spatial-index 选项。以下是如何利用这个空间索引来加速的方法

select tzid
from
  timezones
where
  within(
    MakePoint(cast(:longitude as float), cast(:latitude as float)),
    timezones.Geometry
  ) = 1
  and rowid in (
    select
      rowid
    from
      SpatialIndex
    where
      f_table_name = 'timezones'
      and search_frame = MakePoint(cast(:longitude as float), cast(:latitude as float))
  )

这里的技巧是额外的 and rowid in (...) 子查询。

SpatialIndex 是一个特殊的表 - 它是 SpatiaLite 创建的一个虚拟表。

您可以通过传递与其关联空间索引的另一个表的名称(在本例中是 timezones)来查询该表。然后,您将一个几何图形作为 search_frame 传入,该表将返回一个 rowid 值列表,代表与您传入的几何图形的粗略边界框重叠的任何多边形。

请注意,这并非精确比较:您获得的一些行 ID 可能不会与几何图形精确相交。

但这已经足够近似了。您可以将这些值与更精确的 within() 检查结合起来,这样就只需对多边形总集的一个子集运行完整计算。

在我的本地测试中,这将查询所需的时间从 150ms 降低到不到 8ms - 这是显著的加速!

添加国家多边形

哪些时区与特定国家相关?

目前我们的数据无法告诉我们这一点。我们可以只过滤以 America/ 开头的时区,但这将涵盖北美和南美的所有内容。如果我们能按国家浏览时区,那岂不是很棒?

我们可以尝试通过加载世界各国的多边形来回答这个问题。

datahub.io/core/geo-countries 包含我们所需的数据,该数据源自 Natural Earth 并根据 PDDL 许可证发布。

我们可以像这样下载一个包含国家数据的 GeoJSON 文件

wget https://datahub.io/core/geo-countries/r/countries.geojson

现在我们可以使用 geojson-to-sqlite 将其加载到 countries 数据库表中,并同样创建一个空间索引

pip install geojson-to-sqlite

geojson-to-sqlite timezones.db countries \
  countries.geojson --spatial-index 

Datasette 应该能识别新的 countries 表,并且 datasette-geojson-map 现在会显示这些国家的渲染轮廓。

要查找与特定国家相交的时区,我们可以使用以下查询

select
  tzid,
  geometry
from
  timezones
where intersects(
  timezones.geometry, (
    select geometry from countries where admin = 'United Kingdom'
  ) = 1
)

我们在这里使用了一些额外的技巧。

intersects() 函数是 SpatiaLite 的一个函数,用于检查两个几何图形是否以任何方式相互相交。

我们还使用了子查询来访问我们感兴趣的国家的几何图形

select geometry from countries where admin = 'United Kingdom'

countries 表中的 admin 列包含国家名称 - 这里我们获取的是英国的几何图形。

以下是此查询的结果

简化多边形

尝试将国家名称更改为“United States of America”,您可能会遇到一个问题:与美国相交的时区几何图形非常大,可能无法在地图上渲染!

我们可以使用另一个 SpatiaLite 函数来帮助解决这个问题:simplify()。此函数应用 Douglas–Peucker 算法 将多边形简化为较少数量的点,这使得渲染变得容易得多。

simplify() 接受一个几何图形和一个精度值。经过反复试验,我发现精度值设置为 0.05 对于这些时区多边形效果很好

select
  tzid,
  simplify(geometry, 0.05) as geometry
from
  timezones
where intersects(
  timezones.geometry, (
    select geometry from countries where admin = 'United States of America'
  )  = 1
)

这里必须使用 simplify(...) as geometry,因为 datasette-geojson-map 插件会寻找一个名为 geometry 的输出列以便在地图上渲染。

以下是此查询针对美国的结果

使用索引加速

我们可以使用与上面列出的时区查询类似的方式来使用空间索引

with country as (
  select
    geometry
  from
    countries
  where
    admin = :country
)
select
  timezones.tzid,
  simplify(timezones.geometry, 0.05) as geometry
from
  timezones,
  country
where
  timezones.id in (
    select
      SpatialIndex.rowid
    from
      SpatialIndex,
      country
    where
      f_table_name = 'timezones'
      and search_frame = country.geometry
  )
  and intersects(timezones.geometry, country.geometry) = 1

这使用了与之前相同的技巧 - 一个 where timezones.id in (...) 子查询,它从 SpatialIndex 虚拟表中返回一个可能的 rowid 值列表。

如果一个表有一个整数主键 - 例如我们这里的 timezones 表 - SQLite 会将 rowid 设置为相同的值,这就是为什么我们可以在这个查询中比较 timezones.idrowid

我们在这里使用了另一个技巧。为了避免两次运行查询来选择国家几何图形,我们转而将其捆绑到一个 CTE(公用表表达式)中,放在查询的开头。

with country as (...) 部分使得该查询的结果在 SQL 查询期间可以作为名为 country 的临时表使用。

使用预设查询定义元数据

现在我们已经弄清楚了驱动 API 所需的所有查询,是时候将它们捆绑到一个可以部署到互联网的配置中了。

Datasette 的元数据系统可用于提供有关数据库的额外信息,也可用于配置预设查询 - 带有名称的 SQL 查询,例如本教程开头展示的 by_point 示例。

这是一个 metadata.yml 文件,其中定义了 by_point 查询以及一个按国家代码查找时区的查询。

title: Time zones API
description: |
  An API for looking up time zones by latitude/longitude
about: simonw/timezones-api
about_url: https://github.com/simonw/timezones-api
license: ODbL
license_url: http://opendatacommons.org/licenses/odbl/
source: timezone-boundary-builder
source_url: https://github.com/evansiroky/timezone-boundary-builder
allow_sql: false
databases:
  timezones:
    tables:
      countries:
        source: Natural Earth
        source_url: https://www.naturalearthdata.com/
        license: Open Data Commons Public Domain Dedication and License (PDDL) v1.0
        license_url: https://opendatacommons.org/licenses/pddl/1-0/
        about: geo-countries
        about_url: https://datahub.io/core/geo-countries
    queries:
      by_point:
        title: Find time zone by lat/lon
        sql: |
          select tzid
          from
            timezones
          where
            within(
              MakePoint(cast(:longitude as float), cast(:latitude as float)),
              timezones.Geometry
            ) = 1
            and rowid in (
              select
                rowid
              from
                SpatialIndex
              where
                f_table_name = 'timezones'
                and search_frame = MakePoint(cast(:longitude as float), cast(:latitude as float))
            )
      by_country:
        title: Find time zones that intersect a country
        sql: |
          with country as (
            select
              geometry
            from
              countries
            where
              admin = :country
          )
          select
            timezones.tzid,
            simplify(timezones.geometry, 0.05) as geometry
          from
            timezones,
            country
          where
            timezones.id in (
              select
                SpatialIndex.rowid
              from
                SpatialIndex,
                country
              where
                f_table_name = 'timezones'
                and search_frame = country.geometry
            )
            and intersects(timezones.geometry, country.geometry) = 1

除了提供名为 by_nameby_country 的预设查询外,此文件还包含显示我们用于数据库的数据来源的元数据。

它还设置了另一个重要选项

allow_sql: false

此选项阻止用户对我们发布的数据库执行自己的自定义 SQL 查询。只有我们定义的预设查询可用。

我们在这里使用该选项是因为 SpatiaLite 有大量函数,其中一些可能会导致底层的 Datasette 实例崩溃。

我们可以通过像这样启动 Datasette 来测试新的 metadata.yml 文件

datasette timezones.db --load-extension spatialite -m metadata.yml

将应用程序部署到 Fly

datasette publish 命令可用于将 Datasette 实例部署到各种不同的托管提供商。

Fly 是托管此应用程序的绝佳选择,因为该 API 可能会吸引大量流量。

Google Cloud Run 根据实例持续接收的流量量收费,这对该应用程序来说可能变得昂贵。

Fly 对实例收取固定的月费,外加带宽的额外费用。

此处可查看当前的 Fly 定价。截至撰写本文时,一个拥有 256MB RAM 的实例 - 足以舒适地托管此 API - 每月费用为 1.94 美元。

您需要安装 Fly CLI 工具

curl -L https://fly.io/install.sh | sh

然后运行以下命令向 Fly 进行身份验证

flyctl auth login

您还应该安装 datasette-publish-fly 插件

pip install datasette-publish-fly

准备好所有这些后,您可以像这样部署应用程序

datasette publish fly timezones.db \
  --app timezones-api \
  --setting default_page_size 10 \
  --install datasette-geojson-map \
  --metadata metadata.yml \
  --spatialite

您需要一个唯一的 --app 名称(我已为此演示占用了 timezones-api)。

--spatialite 标志确保为部署的应用程序配置了 SpatiaLite。

部署命令可能需要几分钟才能完成。完成后,您可以访问 https://your-app-name.fly.dev/ 查看完成的应用程序,它已在线上。

以下是 timezones.datasette.io,就是使用此命令部署的。

源代码

本教程中所有内容的源代码都可以在 GitHub 上的 simonw/timezones-api 仓库中找到。