速通PostGIS地理数据分析

前言

因为一些缘故，接触了一下对空间地理数据的分析工作。忙完以后，想了想决定写下这篇文章，方便日后如有需要，可以快速查阅，也供后人参考。

背景条件

我们主要是针对给定的空间地理数据（如GeoJson，Shapefile）进行一些数据提取和分析，获得我们想要的分析结果。本文不讨论这些数据如何收集，抑或是如何编制，我们只讨论分析方面的内容，并且也不会涉及非常深，因为本文重点更在于，快速介绍对地理数据进行分析的大致方法，而具体的分析思路、分析流程，是各位灵感发挥的重要地方。

本文会涉及到对PostgreSQL数据库的使用，建议未接触过PostgreSQL的读者自行查阅相关的基础资料。不过，假如对基本的SQL，特别是MySQL/MariaDB有基本的了解，甚至是一定的使用经验，可以参考本文的前序文章在MySQL基础上速通PostgreSQL

因为我本身不是学地理信息专业的，因此以上知识与理解是临时抱佛脚速成学得，如果本文有什么错误，还请批评指正（请发邮件给我）！多谢各位包涵。

目标

我们的分析工作环境可能需要我们做到以下几点。

• 建立一个PostgreSQL数据库（下面简称psql），并启用PostGIS
• 将收集到的数据放到psql中
• 使用Python对这些数据进行进一步的分析计算

数据库为什么选择PostgreSQL呢？可能是由于PostgreSQL中的GIS插件PostGIS比较优秀，支持比较好。

而选用Python进行数据分析，是因为其提供了Numpy/Pandas/Geopandas/Scipy等十分通用的工具。如果你喜欢其他工具体系，完全可以替换，只要可以连接读取psql中的地理数据即可。这里我们以Geopandas作为主要工具。

环境配置

我觉得配环境向来是浪费口舌的，但又是最重要的一环，没有配好的环境，我们什么都做不了。

我们需要安装的软件有:

• PostgreSQL + PostGIS插件
• Python + Pandas/Numpy/Scipy/SQLAlchemy/Matplotlib模块
• Geopandas/Geoalchemy

一般来说，最上面的两项是比较容易配置的，安装就好了，有非常多的方法，PostgreSQL在Windows上有安装包，在Linux上有软件包管理器。而Python部分无论是使用Conda还是使用Pip，都非常方便。

但有个例外，就是Geopandas/Geoalchemy两项。不知道是怎么回事，Conda源中Geopandas版本号比较古老，而使用conda-forge渠道会把安装卡在环境依赖关系检测部分。我当时是尝试了几次都没弄出理想的效果。

如果你也遇到了这样的问题，或许可以试试下面的方法：容器/Windows Subsystem of Linux

独立的环境

实际上这个方法比较诡异，但至少是能非常顺利地使用的。

我当时是无意中，看到Ubuntu的APT源中，存在python3-geopandas和python3-geoalchemy2两个软件包，版本号也不过时(Debian源中也有，但Debian Bullseye使用的版本号比Ubuntu 22.04的老一些)，遂直接放弃使用Conda配环境的念头，直接安装这两个软件包，运行一切正常。

当然，这样的方法能不出乱子，也得益于这些软件包没有依赖特定的硬件配置或者组件版本号。如果是需要依赖某个特定版本CUDA Toolkit的软件包，或许这样会出现更大的问题。

不过，这样一来，环境隔离就变成了一个需要考虑的问题。我们希望将这些环境隔离起来，万一日后发生错误，可以整个删除，而不对骨干系统造成影响。而APT安装的软件包，是全局性的。

因此，我们可以通过将这部分环境，整个作为一个容器进行隔离。

• 如果你主力操作系统使用的是Windows，那么可以尝试在WSL中安装Ubuntu（目前WSL可以自己手动import发行版，只要空间足够，可以添加非常多个实例，作为各自隔离的子系统，即便它们是同一个发行版，这和应用商店安装的有所不同）。
• 如果你主力操作系统使用的是Linux，那么可以试试LXC/LXD解决方案，直接运行一个Ubuntu容器。Docker应该也是不错的选择。
• 如果你主力操作系统使用的是Mac OS X，那………我也不清楚有什么好的解决方案，毕竟我不是苹果用户，反正思路和上面的是一样的。实在不行还有虚拟机对吧。

创建容器后，就可以直接通过APT安装python3-geopandas和python3-geoalchemy2，以及其他需要的Python模块。

如果你喜欢使用Visual Studio Code 配合 Jupyter Notebook，那么只需要安装python3-ipykernel软件包，并将Visual Studio Code连接到你的WSL或者容器中，打开工作文件夹，创建一个.ipynb文件，选择内核为系统目录下的python，你喜欢的工作流程就呈现在你的面前了。

PostGIS速通

接下来，我们试着介绍以下，如何使用PostGIS对地理数据进行查询。假设我们需要在名为gis_test1的数据库上存储空间数据。

顺带一提，我好像并没有接触3D形状的查询分析，我接触的都是在平面上的地理数据分析。不过，方法几乎是相同的。

启用PostGIS

若想要在一个数据库(Database)中存储地理空间数据，则必须启用PostGIS扩展。

这里的数据库(Database)概念，指的是PostgreSQL服务器下,一个个相对独立的Database，并不是指整个PostgreSQL服务器实例。

我们连接到目标数据库后，使用以下代码启用PostGIS。操作前请确认执行以下语句时，是否已经切换到对应的数据库。

CREATE EXTENSION postgis

GIS中的形状

小学数学课上学过的，点动成线，线动成面，面动成体，在我们三维空间中都遵循着这些道理。

而如果我们将目光投向真实世界平面地图，我们会发现，地图上的许多事物，可以概括为以下几种类型:

• 点(Point): 地图上的一个点，例如，家门口的小卖部，抑或是市内的一座医院，就是一个点
• 线(Linestring): 地图上的一条线，例如，一条没有分支的地铁线路，从头到尾就是一条线
• 多边形(Polygon): 地图上的一个多边形，例如，在中国，一个县级行政区划（市辖区、县级市、县、旗等）所管辖的区域，就是一个多边形

GIS中Polygon的示例 - 来自postgis.net

这几种类型，正是在GIS中，对二维地图上的空间事物进行描述的基本类型。参考9. Geometries — Introduction to PostGIS

在此之上，如果我们将多个同类型的事物，放到一起，就会得到几个集合类型(Collection)，其中几个都是以 Multi 开头的

• 点集(MultiPoint): 地图上的多个点，例如，市内同名连锁便利店的集合，就是一个点集
• 线集(MultiLineString): 地图上的多条线，例如，市内地铁网络的线路集合，就是一个线集
• 多边形集(MultiPolygon): 地图上的多个多边形，例如，在中国，一个地级市下属的所有县级行政区划（市辖区、县级市、县、旗等），就是多个多边形组成的集合
• 地理集合(GeometryCollection): 如你所猜想的，以上几种任意类型的集合，而且也可以包含集合类型

这时候一定会有人问，上文中为什么要将医院、便利店标记成点(Point)？许多地点都有一定的尺寸、占地面积，将其抽象成点是否合适

参照上面参考页面:

Points are used to represent objects when the exact details, such as shape and size, are not important at the target scale. For example, cities on a map of the world can be described as points, while a map of a single state might represent cities as polygons.

Polygons are used to represent objects whose size and shape are important. City limits, parks, building footprints or bodies of water are all commonly represented as polygons when the scale is sufficiently high to see their area. Roads and rivers can sometimes be represented as polygons.

省流版本: 这取决于考察的区域范围与这些地点大小的相对关系，以及这些地点的尺寸是否会对分析结果产生影响。

举个例子，我们如果想要在城市火车站一角，张贴一张地铁线网概览图，那么我们可以放心地将地铁线路抽象为一条条线，而且这些线在图中的粗细，除了影响游客的可读性外，在传达信息的准确性上，并不会带来什么负面影响。

但如果，我们需要在地铁线路附近新建一栋大楼，需要向下挖掘地基，那么我们就不能用上面例子中的概览图来确定施工范围边界！毕竟一条线本身并没有粗细，在考察施工边界这样精度的问题时，这样的一条条线不能告诉我们地铁隧道具体位于何处，也不能告诉我们应该从何处向下挖掘地基。我们必须将地下的铁路隧道视为具有具体粗细的线，也就是一个多边形，才能在地图中准确表现出地铁隧道的位置，从而让我们准确避开地铁隧道进行施工。

PostGIS的形状运算

我们定义出上述元素，自然是为了方便接下来的运算工作。

通常，我们存储到数据库中的数据，常常是地图上事物的原始数据，如城市中公交车站的位置，各建筑物的具体地理边界等。然而，对于我们的分析需求来说，从数据库中读取这些原始数据往往是不够的，我们通常需要进行更复杂的运算，例如查询某个区域内公交车站的数量，或是某个区域内道路的总长度。

PostGIS提供的函数可以帮助我们通过SQL语句，构造上述查询目标的SQL表达式，从而轻松从数据库中获得上述信息。

形状的关系

在中小学时期，我们就学习过，两个图形是否相交，如何相交，会构成不同的模式。在PostGIS中，这一问题进一步和图形的内部、边界、外部结合起来，构成了GIS系统中的九维相交模型(9-Intersection Model)。

具体的解释，可以参考https://postgis.net/docs/manual-3.2/using_postgis_query.html

简单来说，我们常常涉及到的运算，其实主要以计算相交(Intersect)和包含(Contain)关系为主。而PostGIS的函数，通常以ST_开头。上述链接也提供了相关的示例，可以看看。

以我编写的一个SQL查询为例。

SELECT 
    n.area_id,
    n.area_name,
    ca.crime_area,
    ST_Area(ssa2.geometry) AS nei_area
FROM
    neighbourhoods n
JOIN
(
    SELECT
        n.area_id,
        SUM(ST_Area(ST_Intersection(ssa.geometry, bae.geometry))) AS crime_area
    FROM
        neighbourhoods n
    JOIN sa2_statistical_areas ssa ON
        n.area_id = ssa."SA2_MAIN16"
    JOIN break_and_enter bae ON
        ST_Intersects(ssa.geometry,
        bae.geometry)
    GROUP BY
        n.area_id
) ca
ON
    ca.area_id = n.area_id
JOIN sa2_statistical_areas ssa2 
ON
    n.area_id = ssa2."SA2_MAIN16"

在上述查询中，我们用ST_Intersects函数，计算ssa表中的geometry列，与bae表中的geometry列的交集，如果有两条记录产生了交集，则JOIN条件为真，联表查询生成的结果中会多出一行。