SpatiaLite¶

SQLite 的 SpatiaLite 模块增加了处理地理和空间数据的功能。有关如何使用它的示例，请参阅教程《使用 SpatiaLite 构建位置到时区 API》。

要在 Datasette 中使用它，您需要安装 mod_spatialite 动态库。然后可以使用 --load-extension 命令行选项将其加载到 Datasette 中。

如果您这样运行 Datasette，它可以在常见的安装位置查找 SpatiaLite

datasette --load-extension=spatialite --setting default_allow_sql off

如果 SpatiaLite 位于其他位置，请改用扩展的完整路径

datasette --setting default_allow_sql off \
  --load-extension=/usr/local/lib/mod_spatialite.dylib

警告¶

警告

SpatiaLite 扩展添加了大量额外的 SQL 函数，其中一些对于不受信任的用户执行来说可能不安全：它们可能导致 Datasette 服务器崩溃。

除非采取额外的安全措施来防范潜在的有害 SQL 查询，否则不应将启用了 SpatiaLite 的 Datasette 实例暴露在公共互联网上。

建议采取以下步骤

禁止不受信任的用户执行任意 SQL 查询。请参阅控制执行任意 SQL 的能力了解如何执行此操作。最简单的方法是使用 datasette --setting default_allow_sql off 选项启动 Datasette。
使用您希望用户能够执行的、使用了 SpatiaLite 函数的 SQL 查询定义预设查询。

Datasette SpatiaLite 教程包含使用这些技术安全运行 SpatiaLite 的详细说明

安装¶

在 OS X 上安装 SpatiaLite¶

在 OS X 上安装 SpatiaLite 最简单的方法是使用 Homebrew。

brew update
brew install spatialite-tools

这将安装 spatialite 命令行工具和 mod_spatialite 动态库。

您现在可以像这样运行 Datasette

datasette --load-extension=spatialite

在 Linux 上安装 SpatiaLite¶

SpatiaLite 已为大多数 Linux 发行版打包。

apt install spatialite-bin libsqlite3-mod-spatialite

根据您的发行版，您应该可以像这样运行 Datasette

datasette --load-extension=/usr/lib/x86_64-linux-gnu/mod_spatialite.so

如果您不确定模块的位置，请尝试运行 locate mod_spatialite 并查看结果。

对经度/纬度列进行空间索引¶

这里有一个方法：取一个带有现有经度/纬度列的表，向该表添加一个 SpatiaLite POINT 几何列，填充新列，然后填充一个空间索引

import sqlite3

conn = sqlite3.connect("museums.db")
# Lead the spatialite extension:
conn.enable_load_extension(True)
conn.load_extension("/usr/local/lib/mod_spatialite.dylib")
# Initialize spatial metadata for this database:
conn.execute("select InitSpatialMetadata(1)")
# Add a geometry column called point_geom to our museums table:
conn.execute(
    "SELECT AddGeometryColumn('museums', 'point_geom', 4326, 'POINT', 2);"
)
# Now update that geometry column with the lat/lon points
conn.execute(
    """
    UPDATE museums SET
    point_geom = GeomFromText('POINT('||"longitude"||' '||"latitude"||')',4326);
"""
)
# Now add a spatial index to that column
conn.execute(
    'select CreateSpatialIndex("museums", "point_geom");'
)
# If you don't commit your changes will not be persisted:
conn.commit()
conn.close()

利用空间索引¶

SpatiaLite 空间索引是 R*Tree。它们允许您使用子查询高效地运行边界框查询，模式类似于使用自定义 SQL 进行搜索中使用的模式。

在上面的例子中，生成的索引将被称为 idx_museums_point_geom。这是一种 SQLite 虚拟表。您可以使用以下查询检查其内容

select * from idx_museums_point_geom limit 10;

这里有一个在线示例：timezones-api.datasette.io/timezones/idx_timezones_Geometry

pkid	xmin	xmax	ymin	ymax
1	-8.601725578308105	-2.4930307865142822	4.162120819091797	10.74019718170166
2	-3.2607860565185547	1.27329421043396	4.539252281188965	11.174856185913086
3	32.997581481933594	47.98238754272461	3.3974475860595703	14.894054412841797
4	-8.66890811920166	11.997337341308594	18.9681453704834	37.296207427978516
5	36.43336486816406	43.300174713134766	12.354820251464844	18.070993423461914

您现在可以构建利用索引的高效边界框查询，如下所示

select * from museums where museums.rowid in (
    SELECT pkid FROM idx_museums_point_geom
    -- left-hand-edge of point > left-hand-edge of bbox (minx)
    where xmin > :bbox_minx
    -- right-hand-edge of point < right-hand-edge of bbox (maxx)
    and xmax < :bbox_maxx
    -- bottom-edge of point > bottom-edge of bbox (miny)
    and ymin > :bbox_miny
    -- top-edge of point < top-edge of bbox (maxy)
    and ymax < :bbox_maxy
);

空间索引可以针对多边形列和点列创建，在这种情况下，它们将代表该多边形的最小边界矩形。这对于加速 within 查询非常有用，如 Timezones API 示例所示。

将 shapefile 导入 SpatiaLite¶

shapefile 格式是分发地理空间数据的一种常见格式。您可以使用 spatialite 命令行工具从 shapefile 创建新的数据库表。

现在尝试使用北卡罗来纳大学全球河流数据库项目提供的北美 shapefile。下载文件并解压缩（这将在当前目录中创建名为 narivs.dbf、narivs.prj、narivs.shp 和 narivs.shx 的文件），然后运行以下命令

$ spatialite rivers-database.db
SpatiaLite version ..: 4.3.0a       Supported Extensions:
...
spatialite> .loadshp narivs rivers CP1252 23032
========
Loading shapefile at 'narivs' into SQLite table 'rivers'
...
Inserted 467973 rows into 'rivers' from SHAPEFILE

这将把 narivs shapefile 中的数据加载到一个名为 rivers 的新数据库表中。

退出 spatialite（使用 Ctrl+D），然后像这样针对您的新数据库运行 Datasette

datasette rivers-database.db \
    --load-extension=/usr/local/lib/mod_spatialite.dylib

如果您浏览到 http://localhost:8001/rivers-database/rivers，您将看到新表... 但 Geometry 列将包含无法读取的二进制数据（SpatiaLite 使用基于 WKB 的自定义格式）。

将其转换为半可读数据最简单的方法是使用 SpatiaLite 的 AsGeoJSON 函数。尝试在 http://localhost:8001/rivers-database 的 SQL 查询界面中使用以下查询

select *, AsGeoJSON(Geometry) from rivers limit 10;

这将为您提供一个额外的 GeoJSON 列。您可以将此列中的 GeoJSON 复制并粘贴到 geojson.io 的调试工具中，以便在地图上将其可视化。

要查看一个更有趣的示例，请尝试首先按几何长度最长的记录进行排序。由于表中有 467,000 行，您首先需要增加 Datasette 设置的 SQL 时间限制

datasette rivers-database.db \
    --load-extension=/usr/local/lib/mod_spatialite.dylib \
    --setting sql_time_limit_ms 10000

现在尝试以下查询

select *, AsGeoJSON(Geometry) from rivers
order by length(Geometry) desc limit 10;

使用 Shapely 导入 GeoJSON 多边形¶

多边形数据的另一种常见形式是 GeoJSON 格式。可以直接将其导入 SpatiaLite，或使用Shapely Python 库。

Who's On First 是一个极好的开放许可 GeoJSON 多边形数据源。让我们导入威尔士的地理多边形数据。首先，我们可以使用 Who's On First Spelunker 工具查找威尔士的记录

spelunker.whosonfirst.org/id/404227475

该页面包含一个指向 GeoJSON 记录的链接，可在此处访问

data.whosonfirst.org/404/227/475/404227475.geojson

以下是用于创建 SQLite 数据库、启用 SpatiaLite、创建 places 表并为威尔士添加记录的 Python 代码

import sqlite3

conn = sqlite3.connect("places.db")
# Enable SpatialLite extension
conn.enable_load_extension(True)
conn.load_extension("/usr/local/lib/mod_spatialite.dylib")
# Create the masic countries table
conn.execute("select InitSpatialMetadata(1)")
conn.execute(
    "create table places (id integer primary key, name text);"
)
# Add a MULTIPOLYGON Geometry column
conn.execute(
    "SELECT AddGeometryColumn('places', 'geom', 4326, 'MULTIPOLYGON', 2);"
)
# Add a spatial index against the new column
conn.execute("SELECT CreateSpatialIndex('places', 'geom');")
# Now populate the table
from shapely.geometry.multipolygon import MultiPolygon
from shapely.geometry import shape
import requests

geojson = requests.get(
    "https://data.whosonfirst.org/404/227/475/404227475.geojson"
).json()
# Convert to "Well Known Text" format
wkt = shape(geojson["geometry"]).wkt
# Insert and commit the record
conn.execute(
    "INSERT INTO places (id, name, geom) VALUES(null, ?, GeomFromText(?, 4326))",
    ("Wales", wkt),
)
conn.commit()

使用 within() 查询多边形¶

within() SQL 函数可用于检查一个点是否在一个几何图形内

select
    name
from
    places
where
   within(GeomFromText('POINT(-3.1724366 51.4704448)'), places.geom);

GeomFromText() 函数接受一个 well-known text 字符串。请注意，此处使用的顺序是 longitude 然后 latitude。

要以利用空间索引的方式运行相同的 within() 查询，请使用以下命令

select
    name
from
    places
where
    within(GeomFromText('POINT(-3.1724366 51.4704448)'), places.geom)
    and rowid in (
        SELECT pkid FROM idx_places_geom
        where xmin < -3.1724366
        and xmax > -3.1724366
        and ymin < 51.4704448
        and ymax > 51.4704448
    );