PostGIS 笔记

空间数据库

定义：存储、操作空间数据的数据库。

空间数据特点

• 空间数据类型是指点、线、多边形等形状；
• 提供用于空间操作的高效处理的多维空间索引；
• 提供用于查询空间属性和关系额空间函数。

发展历程

• 第1代
  • 在传统的 GIS 实施中，所有空间数据都存储在平面文件中，并且需要特殊的 GIS 软件来解释和操作数据；
  • 这些第一代管理系统旨在满足所有所需数据都在用户组织域内的用户的需求；
  • 它们是专为处理空间数据而构建的专有、独立的系统。
• 第2代
  • 空间系统将一些数据存储在关系数据库中（通常是属性或非空间部分），但仍然缺乏直接集成所提供的灵活性。
• 第3代
  • 当人们开始将空间特征视为主流数据库对象时，真正的空间数据库诞生了；
  • 空间数据库将空间数据与关系数据库完全集成；
  • 系统定位从以 GIS 为中心转变为以数据库为中心。

空间数据库管理系统 可用于地理世界以外的应用，用于管理与人体解剖学、大规模集成电路、分子结构和电磁场等相关的数据。

空间数据类型

传统数据库有字符串、数字和日期等类型；空间数据库添加了额外的空间类型来表示地理特征，这些空间数据类型抽象和封装了边界和维度等空间结构。在许多方面，空间数据类型可以简单地理解为形状。

空间数据类型以类型层次结构进行组织，每个子类型都继承其超类型的结构（属性）和行为（方法或函数）。

空间索引和边界框

普通数据库提供索引以允许对数据子集进行快速和随机访问，标准类型（数字、字符串、日期）的索引通常使用B树索引完成。

• B树使用自然排序顺序对数据进行分区，以将数据放入分层树中；
• 数字、字符串和日期的自然排序顺序很容易确定——每个值都小于、大于或等于所有其他值。

但是由于多边形可以重叠，可以相互包含，并且可以排列在二维（或更多）空间中，因此不能使用B树来有效地索引它们。真实的空间数据库提供了一个“空间索引”来回答“哪些对象在这个特定的边界框内？”这个问题。

边界框是最小的矩形(平行于坐标轴)能够包含给定的特征。

使用边界框是因为回答A 包含 B 吗？的问题。对于多边形来说计算量很大，但在矩形的情况下速度很快。即使是最复杂的多边形和线串也可以用一个简单的边界框来表示。

索引必须快速执行才能有用，空间索引不像B树那样提供精确的结果，而是提供近似的结果。问题这个多边形里面有什么线？而是由空间索引解释为哪些线的边界框包含在该多边形的边界框内？

各种数据库实现空间索引差别很大

• 最常见的实现是R-Tree和Quadtree（用于 PostGIS）；
• 在其他空间数据库中实现的基于网格的索引和GeoHash索引。

空间函数

普通数据库提供了一些查询数据功能，如连接字符串、对字符串执行哈希运算、对数字进行数学运算以及从日期中提取信息等。

空间数据库提供了一整套用于分析几何组件、确定空间关系、操纵几何的功能；这些空间功能是任何空间项目的基石。

可进一步划分为

• 转换：在几何和外部数据格式之间转换的函数；
• 管理：管理有关空间表和 PostGIS 管理的信息的功能；
• 检索：检索几何的属性和测量值的函数；
• 比较：比较两个几何形状的空间关系的函数；
• 生成：从其他函数生成新几何图形的函数。

可能的函数列表非常多，PostGIS 实现的函数（以及其他有用的函数）由OGC、SFSQL定义。

PostGIS

PostGIS 向 PostgreSQL 添加对空间类型、空间索引、空间函数功能，将其转变为空间数据库，并继承了PostGIS企业特性以及实施的开放标准。

PostgreSQL 是什么

PostgreSQL 是一个功能强大的关系数据库管理系统 (RDBMS)，是在 BSD 风格的许可证下发布的，因此是免费的开源软件。与许多其他开源程序一样，PostgreSQL 不受任何一家公司控制，而是拥有一个由开发人员和公司组成的全球社区来维护。

PostgreSQL 在设计之初就考虑了类型扩展（在运行时添加新数据类型、函数和索引的能力）。正因为如此，PostGIS 扩展可以由一个单独的开发团队开发，但仍然非常紧密地集成到核心 PostgreSQL 数据库中。

为什么选择 PostgreSQL

熟悉开源数据库的人的一个常见问题是，为什么 PostGIS 不是基于 MySQL 构建的？。

PostgreSQL 有：

• 默认情况下经过验证的可靠性和事务完整性 (ACID)；

• 对 SQL 标准的仔细支持（完整的 SQL92）；

• 可插拔类型扩展和功能扩展；

• 面向社区的发展模式；

• 对列大小没有限制（TOAST元组）以支持大型 GIS 对象；

• 允许 R-Tree 索引的通用索引结构 (GiST)；

• 轻松添加自定义功能；

结合起来，PostgreSQL 提供了一个非常简单的开发路径来添加新的空间类型。在专有领域，只有Illustra（现为 Informix Universal Server）允许如此简单的扩展。这不是巧合，Illustra 是对 1980 年代原始 PostgreSQL 代码库的专有改造。

因为向 PostgreSQL 添加类型的开发路径非常简单，所以从那里开始是有意义的。当 MySQL 在 4.1 版中发布基本空间类型时，PostGIS 团队查看了他们的代码，加强了最初选择 PostgreSQL 的决定。

因为 MySQL 空间对象作为一种特殊情况必须在字符串类型之上被破解，所以 MySQL 代码分布在整个代码库中。PostGIS 0.1 的开发耗时不到一个月。做一个MyGIS0.1 会花费更长的时间，因此，可能永远不会看到曙光。

为什么不是文件

自从首次编写 GIS 软件以来， Shapefile（以及 Esri File Geodatabase 和GeoPackage等其他格式）一直是存储空间数据并与空间数据交互的标准方式。

但是，这些平面文件具有以下缺点：

• 文件需要特殊的软件才能读取和写入：SQL 是随机数据访问和分析的抽象。如果没有这种抽象，您将需要自己编写所有访问和分析代码。
• 并发用户可能会导致损坏和速度下降：虽然可以编写额外的代码来确保对同一个文件的多次写入不会损坏数据，但当您解决了问题并解决了相关的性能问题时，您将编写数据库系统的大部分内容。为什么不只使用标准数据库？
• 复杂的问题需要复杂的软件来解决：在针对文件进行编程时，可以用数据库中的一行 SQL 表达的复杂而有趣的问题（空间连接、聚合等）需要数百行专门的代码才能回答。

PostGIS 的大多数用户都在设置需要多个应用程序访问数据的系统，因此使用标准 SQL 访问方法可以简化部署和开发。一些用户正在处理大型数据集；对于文件，它们可能被分割成多个文件，但在数据库中它们可以存储为单个大表。

总之，对多用户的支持、复杂的即席查询、大型数据集的性能是空间数据库与基于文件的系统不同的地方。

PostGIS 简史

2001 年 5 月，Refractions Research 发布了 PostGIS 的第一个版本。PostGIS 0.1 有对象、索引和一些函数。结果是一个适合存储和检索但不适合分析的数据库。

随着功能数量的增加，对组织原则的需求变得清晰起来。来自开放地理空间联盟的 Simple Features for SQL（SFSQL）规范为这种结构提供了函数命名和要求的指南。

借助 PostGIS 对简单分析和空间连接的支持，Mapserver成为第一个在数据库中提供数据可视化的外部应用程序。

在接下来的几年中，PostGIS 功能的数量不断增加，但其功能仍然有限。许多最有趣的函数（例如，ST_Intersects()、ST_Buffer()、ST_Union()）都很难编码。从头开始编写它们承诺了多年的工作。

幸运的是，第二个项目“几何引擎，开源”或GEOS出现了。GEOS 库为实现SFSQL规范提供了必要的算法。通过在 GEOS 中的链接，PostGIS 在0.8 版本中提供了对SFSQL的完整支持。

随着 PostGIS 数据容量的增长，另一个问题浮出水面：用于存储几何图形的表示证明效率相对较低。对于像点和短线这样的小对象，表示中的元数据有高达 300% 的开销。出于性能原因，有必要让代表节食。通过缩小元数据标头和所需的维度，大大减少了开销。在 PostGIS 1.0 中，这种新的、更快的、轻量级的表示成为默认值。

PostGIS 的最新版本继续添加功能和性能改进，以及对 PostgreSQL 核心系统中的新功能的支持。

谁使用 PostGIS

有关案例研究的完整列表，请参阅PostGIS 案例研究页面。

法国国家地理研究所

IGN 是法国的国家测绘机构，使用 PostGIS 存储该国的高分辨率地形图BDUni。BDUni 拥有超过 1 亿个特征，由 100 多名现场工作人员维护，他们每天验证观察结果并向数据库添加新映射。IGN 安装使用数据库事务系统来确保更新过程中的一致性，并使用热备用系统在系统发生故障时保持正常运行时间。

RedFin

RedFin是一家房地产代理公司，提供基于网络的服务，用于探索房产和估价。他们的系统最初是基于 MySQL 构建的，但他们发现迁移到 PostgreSQL 和 PostGIS在性能和可靠性方面提供了巨大的好处。

哪些应用程序支持 PostGIS

PostGIS 已成为一种广泛使用的空间数据库，支持使用它存储和检索数据的第三方程序的数量也在增加。支持 PostGIS 的程序包括服务器和桌面系统上的开源软件和专有软件。

下表列出了一些利用 PostGIS 的软件

开放/免费	封闭/专有
加载/提取 Shp2Pgsql ogr2ogr Dxf2PostGIS 基于网络 Mapserver GeoServer (Java-based WFS / WMS -server ) SharpMap SDK - for ASP.NET 2.0 MapGuide Open Source (using FDO) 桌面 uDig QGIS mezoGIS OpenJUMP OpenEV SharpMap SDK for Microsoft.NET 2.0 ZigGIS for ArcGIS/ArcObjects.NET GvSIG GRASS	加载/提取 安全 FME 桌面 翻译器/转换器 基于网络 Ionic Red Spider (now ERDAS) Cadcorp GeognoSIS Iwan Mapserver MapDotNet Server MapGuide Enterprise (using FDO) ESRI ArcGIS Server 桌面 Cadcorp SIS Microimages TNTmips GIS ESRI ArcGIS Manifold GeoConcept MapInfo (v10) AutoCAD Map 3D (using FDO)

PostGIS 安装

探索 PostgreSQL/PostGIS 数据库，并了解如何在 SQL 中编写空间查询，我们需要一些软件，可以在本地安装或在云端远程使用。

• 下面有关于如何访问 PostgreSQL 以在 Windows 或 MacOS 上安装的说明。PostgreSQL for Windows 和 MacOS 要么包括 PostGIS，要么有一种简单的方法来添加它。
• 下面有关于如何安装PgAdmin的说明。PgAdmin 是一个图形数据库浏览器和 SQL 编辑器，它为全世界的数据库引擎提供面向用户的界面。

有关安装 PostgreSQL 的最新说明，请转到PostgreSQL 下载页面并选择您正在使用的操作系统。

Windows

对于 Windows 安装：

1. 转到Windows PostgreSQL 下载页面。

2. 选择最新版本的 PostgreSQL 并将安装程序保存到磁盘。

3. 运行安装程序并接受默认值。

4. 查找并运行随数据库安装的StackBuilder程序。

5. 选择空间扩展部分并选择最新的PostGIS ..Bundle选项。
   

6. 接受默认值并安装。

MacOS

对于 MacOS 安装：

1. 转到Postgres.app站点，并下载最新版本。

2. 打开磁盘映像，然后将Postgres图标拖到Applications文件夹中。
   

3. 在Applications文件夹中，双击Postgres图标以启动服务器。

4. 单击Initialize按钮创建一个新的空白数据库实例。

   

5. 在Applications文件夹中，转到Utilities文件夹并打开Terminal。

6. 为方便起见，将命令行实用程序添加到您的PATH中。

sudo mkdir -p /etc/paths.d
echo /Applications/Postgres.app/Contents/Versions/latest/bin | sudo tee /etc/paths.d/postgresapp

客户端 DBeaver

DBeaver是一个面向开发人员和数据库管理员的免费开源通用数据库工具。

支持所有流行的数据库：MySQL、PostgreSQL、SQLite、Oracle、DB2、SQL Server、Sybase、MS Access、Teradata、Firebird、Apache Hive、Phoenix、Presto 等。

选择原因

支持空间数据查询结果在地图上落图！

安装

• 下载适合自己操作系统的软件版本；
• Windows 平台仅需选择安装路径，一直下一步就行，此处不再赘述；
• 其他平台请参考官方安装说明。

创建空间数据库

管理程序

1. 启动 DBeaver
   

2. 连接 PostgreSQL 数据库
   菜单选择数据库→新建数据库连接
   
   双击 PostgreSQL 图标，进入 PostgreSQL 连接配置界面，按实际情况填写即可。
   
   注意：若第一次使用缺少数据库JDBC连接驱动，该工具将引导你下载。
   

创建数据库

1. 右键选择新建数据库
   

2. 创建名为nyc的数据库
   

3. 选择新nyc数据库并将其打开以显示对象树，您将看到public架构。
   

4. 新建SQL编辑器，用于执行SQL。
   

5. 在加载 PostGIS 空间扩展插件

CREATE EXTENSION postgis;

8. 运行 PostGIS 函数确认 PostGIS插件安装成功

SELECT postgis_full_version();

函数列表

PostGIS_Full_Version()：报告完整的 PostGIS 版本和构建配置信息。

加载空间数据

PostGIS 通过各种库、应用程序的支持，提供了许多加载数据的选项。

首先从数据库备份文件中加载我们的工作数据，然后回顾一些使用常用工具加载不同 GIS 数据格式的标准方法。

示例数据下载

示例数据postgis-workshop-2020.zip为纽约市的四个 ShapeFile，以及一个社会人口变量属性表。我们已将 ShapeFile 加载为 PostGIS 表，稍后将在研讨中添加社会人口统计数据。

下面描述了我们每个数据集的记录数和表属性，这些属性值和关系是我们未来分析的基础。

要在 pgAdmin 中探索表的性质，请右键单击突出显示的表并选择Properties。您将找到表格属性的摘要，包括列选项卡中的表格属性列表。

nyc_census_blocks 人口普查区

人口普查区，是报告人口普查数据的最小地理区域，所有更高级别的人口普查地理（区块组、区域、都会区、县等）都可以从人口普查区块的联合构建。我们已将一些人口统计数据附加到我们的区块集合中。

记录数：36592

blkid	唯一标识
popn_total	总人数
popn_white	白人人数
popn_black	黑人人数
popn_nativ	美洲原住民人数
popn_asian	亚洲人人数
popn_other	其他类别人数
boroname	纽约自治市的名称（曼哈顿、布朗克斯、布鲁克林、史泰登岛、皇后区）
geom	空间数据

黑人人口占总人口的百分比

要将人口普查数据导入 GIS，您需要连接两条信息：实际数据（文本）和边界文件（空间）。获取数据有很多选择，包括从人口普查局的American FactFinder下载数据和边界。

nyc_neighborhoods 纽约社区

纽约拥有丰富的社区名称和范围历史，社区是不遵循政府规定的社会结构。例如，卡罗尔花园、红胡克和鹅卵石山的布鲁克林社区曾被统称为南布鲁克林。现在，根据您与哪个房地产经纪人交谈，该街区以前称为红钩的相同四个街区可以称为哥伦比亚高地、卡罗尔花园西或红钩！

记录数：129

name	小区名称
boroname	纽约自治市的名称（曼哈顿、布朗克斯、布鲁克林、史泰登岛、皇后区）
geom	空间数据

纽约市的街区

nyc_streets 纽约街道

街道中心线形成了城市的交通网络，这些街道标有类型，以便区分后巷、主干道、高速公路和较小的街道等通道。理想的居住区域可能位于住宅街道上，而不是靠近高速公路。

记录数：19091

name	街道名称
oneway	街道是单向的吗？
type	道路类型（主要、次要、住宅、高速公路）
geom	空间数据

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GqDOCXCF-1645930975448)(https://postgis.net/workshops/postgis-intro/_images/nyc_streets.png)]
纽约市的街道。主要道路为红色。

nyc_subway_stations 纽约地铁站

地铁站将人们居住的上层世界与地下无形的地铁网络连接起来，作为公共交通系统的门户，车站位置有助于确定不同的人进入地铁系统的难易程度。

记录数：491

name	车站名称
borough	纽约自治市的名称（曼哈顿、布朗克斯、布鲁克林、史泰登岛、皇后区）
routes	贯穿本站的地铁线路
transfers	您可以通过该站换乘的线路
express	是否为特快列车停靠的车站
geom	空间数据，车站点位

纽约市地铁站的点位置

nyc_census_sociodata 纽约市人口普查数据

在普查过程中收集了丰富的社会经济数据，但仅限于人口普查区域的更大地理层面。人口普查块组合形成人口普查区域（和块组）。我们在人口普查区收集了一些社会经济信息，以回答有关纽约市的一些更有趣的问题。

是nyc_census_sociodata一个数据表。在进行任何空间分析之前，我们需要将其连接到人口普查地理。

tractid	唯一标识
transit_total	工人数量
transit_private	使用私家车/摩托车的工人人数
transit_public	乘坐公共交通的工人人数
transit_walk	步行的工人数量
transit_other	使用其他形式（如步行/骑自行车）的工人数量
transit_none	在家工作的工人人数
transit_time_mins	所有工人在运输过程中花费的总分钟数（分钟）
family_count	家庭数量
family_income_median	家庭收入中位数（美元）
family_income_mean	家庭收入平均数（美元）
family_income_aggregate	所有家庭的总收入（美元）
edu_total	有教育经历的人数
edu_no_highschool_dipl	没有高中文凭的人数
edu_highschool_dipl	具有高中文凭且没有继续教育的人数
edu_college_dipl	具有大专学历且没有继续深造的人数
edu_graduate_dipl	拥有研究生学历的人数

导入示例数据

通过PostGIS 自带的PostGIS PostGIS Bundle 3 for PostgreSQL x64 13 Shapefile and DBF Loader Exporter软件导入 ShapeFile 示例数据。

1. 打开PostGIS PostGIS Bundle并登录nyc数据库
   
   
2. 批量添加要导入的Shape File文件，修改其SRID为26918，选择Import导入
   
3. 导入完成
   

ShapeFile 文件

ShapeFile 通常是指具公共前缀名称的 .shp、.shx、.dbf和其他扩展名的文件的集合。实际的 ShapeFile 专门与具有.shp扩展名的文件相关，但是.shp如果没有所需的支持文件，单独的文件对于分发是不完整的。

强制性文件

• .shp 形状格式；要素几何本身
• .shx 形状索引格式；特征几何的位置索引
• .dbf 属性格式；每个形状的柱状属性，在 dBase III 中
  可选文件包括
• .prj 投影格式；坐标系和投影信息，使用众所周知的文本格式描述投影的纯文本文件

通过shp2pgsql程序通过将 ShapeFile 文件从二进制数据 转换 为一系列 SQL 命令，然后加载到 PostGIS 数据库中运行。

SRID 是什么

空间参考标识符（SRID）定义了数据的地理坐标系和投影的所有参数。SRID 很方便，因为它将有关地图投影（可能非常复杂）的所有信息打包到一个数字中。

• 可以在在线数据库中查找地图投影的定义;
• 或直接在 PostGIS 内部查询spatial_ref_sys表。

SELECT srtext FROM spatial_ref_sys WHERE srid = 26918;

PostGISspatial_ref_sys表是一个OGC标准表，它定义了数据库已知的所有空间参考系统。PostGIS 附带的数据列出了 3000 多个已知的空间参考系统以及在它们之间转换/重新投影所需的详细信息。
在这两种情况下，都会看到26918空间参考系统的文本表示：

PROJCS["NAD83 / UTM zone 18N",
  GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
      SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
      AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4269"]],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["latitude_of_origin",0],
  PARAMETER["central_meridian",-75],
  PARAMETER["scale_factor",0.9996],
  PARAMETER["false_easting",500000],
  PARAMETER["false_northing",0],
  AUTHORITY["EPSG","26918"],
  AXIS["Easting",EAST],
  AXIS["Northing",NORTH]]

如果从数据目录打开nyc_neighborhoods.prj文件，将看到相同的投影定义。

从当地机构（例如纽约市）收到的数据通常会采用state plane或UTM标注的当地投影。我们的投影是Universal Transverse Mercator (UTM) Zone 18 North或 EPSG:26918。

SQL 基本查询

结构化查询语言（SQL: Structured Query Language）是一种向关系数据库查询、更新数据的方法。

通过以下四个基本的动词，就可以使用SQL语句来进行查询：

• SELECT：返回行以响应查询，详情；
• INSERT：向表中添加新行；
• UPDATE：更改表中的现有行;
• DELETE：从表中删除行。

现在通过 DBeaver 使用 SQL 来查询数据。

📌示例：纽约市所有街区的名称是什么？

SELECT name FROM nyc_neighborhoods;

选择查询

选择查询语法如下：

SELECT
	-- 列名
	some_columns
	-- 要查询的表
FROM some_data_source
-- 条件过滤器
WHERE some_condition;

📌示例：布鲁克林所有街区的名称是什么？

SELECT name
FROM nyc_neighborhoods
-- 过滤街区名为'Brooklyn'的记录
WHERE boroname = 'Brooklyn';

📌示例：布鲁克林所有街区的名称中有多少个字母？

SELECT
	name,
	-- 调用符串长度函数 char_length(string)
	char_length(name)
FROM nyc_neighborhoods
WHERE boroname = 'Brooklyn';

PostgreSQL 内置一系列的聚合函数 。

📌示例：布鲁克林所有街区名称中字母的平均数量和字母数量的标准差是多少？

SELECT
	-- 用于平均值的通用 avg()
	avg(char_length(name)),
	-- 用于标准偏差的 stddev()
	stddev(char_length(name))
FROM nyc_neighborhoods
WHERE boroname = 'Brooklyn';

如果希望对整个结果集进行分组汇总，聚合函数还需要搭配GROUP BY语句来使用。

📌示例：按区报告，纽约市所有街区名称中的平均字母数是多少？

SELECT 
	boroname,
	avg(char_length(name)),
	stddev(char_length(name))
FROM nyc_neighborhoods
-- 分组
GROUP BY boroname;

函数列表

• count(expression)：返回一组记录中的记录数；
• avg(expression)：返回数值列的平均值；
• sum(expression)：返回一组记录中记录的总和；
• char_length(string) : 返回字符串中的字符数；
• stddev(expression)：返回输入值的标准差。

练习

使用nyc_census_blocks表格联系SQL查询。

📌练习：nyc_streets 表中有多少条记录？

SELECT
	Count(*)
FROM nyc_streets;

📌练习：纽约市有多少条街道以B开头？

SELECT 
	Count(*)
FROM nyc_streets
WHERE name LIKE 'B%';

📌练习：纽约市的人口有多少？

SELECT
	Sum(popn_total) AS population
FROM nyc_census_blocks;

可以使用AS为表或列定义别名，使查询更易于编写和阅读。

📌练习：布朗克斯有多少人口？

SELECT
	Sum(popn_total) AS population
FROM nyc_census_blocks
WHERE boroname = 'The Bronx';

📌练习：每个行政区有多少个社区？

SELECT
	boroname,
	count(*)
FROM nyc_neighborhoods
GROUP BY boroname;

📌练习：对于每个行政区，白人占人口的百分比是多少？

SELECT
    boroname,
    (Sum(popn_white)/Sum(popn_total) * 100.0) AS white_pct
FROM nyc_census_blocks
GROUP BY boroname;

几何

示例数据

执行以下代码，构建示例数据。

CREATE TABLE geometries (name varchar, geom geometry);

INSERT INTO geometries VALUES
  ('Point', 'POINT(0 0)'),
  ('Linestring', 'LINESTRING(0 0, 1 1, 2 1, 2 2)'),
  ('Polygon', 'POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))'),
  ('PolygonWithHole', 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(1 1, 1 2, 2 2, 2 1, 1 1))'),
  ('Collection', 'GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0, 1 0, 1 1, 0 1, 0 0)))');

SELECT name, ST_AsText(geom) FROM geometries;

上面的示例创建一个表（geometries）然后插入五个几何：一个点、一条线、一个多边形、一个带孔的多边形和一个集合。

元数据表

根据SFSQL规范，PostGIS 提供了两个表（及视图）来跟踪、报告给定数据库中可用的几何类型。

• spatial_ref_sys：定义了数据库已知的所有空间参考系统；
• geometry_columns：提供了所有特征（定义为具有几何属性的对象）的列表。

让我们看一下geometry_columns我们数据库中的表。像以前一样将此命令粘贴到查询工具中：

SELECT * FROM geometry_columns;

• f_table_catalog：表示数据库；
• f_table_schema：表示数据库实例；
• f_table_name：表示表名；
• f_geometry_column：是包含几何列的列的名称，对于具有多个几何列的特征表，每个列将有一条记录；
• coord_dimension：表示存储空间数据的维度（2 维、3 维或 4 维）；
• srid：即SRID，spatial_ref_sys表引用该表的空间参考系统标识符。
• type：定义几何类型。

通过查询此表，GIS 客户端 / 库可以在执行投影、处理、渲染时确定检索数据的预期内容，无需检查每个几何图形。

若缺失srid，可以通过以下方式更新。

ALTER TABLE nyc_neighborhoods
ALTER COLUMN geom
TYPE Geometry(MultiPolygon, 26918)
USING ST_SetSRID(geom, 26918);

表示真实世界的对象

SFSQL 规范是 PostGIS 开发的原始指导标准，它定义了真实世界对象的表示方式。通过采用连续形状并以固定分辨率对其进行数字化，实现对象表示。SFSQL 只处理二维表示；PostGIS 已将其扩展到包括 3 维和 4 维表示；最近，SQL-Multimedia Part 3 ( SQL/MM) 规范正式定义了它们自己的表示。

示例表包含不同几何类型的混合，可以使用读取几何元数据的函数收集有关每个对象的一般信息。

• ST_GeometryType(geometry)返回几何的类型；
• ST_NDims(geometry)返回几何的维数；
• ST_SRID(geometry)返回几何的空间参考标识符号。

SELECT
	name,
	ST_GeometryType(geom),
	ST_NDims(geom),
	ST_SRID(geom)
FROM geometries;

点

一个空间点代表地球上的一个位置，由单个坐标表示。如，世界地图上的城市可以用点来描述，而单个州的地图可以用多边形来表示城市。

SELECT
	ST_AsText(geom)
FROM geometries
WHERE name = 'Point';

用于处理点的特定空间函数

• ST_X(geometry)：返回 X 纵坐标；
• ST_Y(geometry)：返回 Y 坐标。

因此，可以从这样的点读取坐标：

SELECT
	ST_X(geom),
	ST_Y(geom)
FROM geometries
WHERE name = 'Point';

纽约市地铁站 ( nyc_subway_stations) 表是一个以点表示的数据集。以下 SQL 查询将返回与一个点关联的几何图形（在ST_AsText列中）。

SELECT
	name,
	ST_AsText(geom)
FROM nyc_subway_stations
LIMIT 1;

线

线是位置之间的路径，采用两个或多个有序点表示，如道路和河流通常表示为线。

• 如果线在同一点开始和结束，则称它是闭合的。
• 如果它不与自身交叉或接触，则称它是简单的（如果它是封闭的，则在其端点处除外）。

纽约街道表 ( nyc_streets) 数据包含线名称和类型等详细信息。一条现实世界的街道可能由许多线串组成，每个线串代表一段具有不同属性的道路。

SELECT ST_AsText(geom)
FROM geometries
WHERE name = 'Linestring';

用于处理线的特定空间一些函数是

• ST_Length(geometry)：返回线串的长度；
• ST_StartPoint(geometry)：返回第一个坐标作为一个点；
• ST_EndPoint(geometry)：返回最后一个坐标作为一个点；
• ST_NPoints(geometry)：返回线串中的坐标数。

线串的长度

SELECT ST_Length(geom)
FROM geometries
WHERE name = 'Linestring';

面

面是一个区域的表示，外边界由一个环表示。

• 外圈（环）：是一个封闭、简单的线；
• 内圈（孔）：也是一个封闭、简单的线。

面可以用于表示空间对象的大小、形状。当比例尺足够高以查看其区域时，城市范围、公园、建筑足迹或水体通常都表示为多边形，如道路和河流有时可以表示为多边形。

SELECT ST_AsText(geom)
FROM geometries
WHERE name LIKE 'Polygon%';

第一个多边形只有一个环，第二个有一个内部洞。大多数图形系统都包含面的概念，但 GIS 系统在允许多边形显式具有孔方面相对独特。

用于处理多边形的特定空间一些函数是

• ST_Area(geometry)：返回多边形的面积；
• ST_NRings(geometry)：返回环的数量（通常为 1，如果有孔则更多）；
• ST_ExteriorRing(geometry)：将外环作为线串返回；
• ST_InteriorRingN(geometry,n)：将指定的内部环作为线串返回；
• ST_Perimeter(geometry)：返回所有环的长度。

可以使用 area 函数计算多边形的面积：

SELECT 
	name,
	ST_Area(geom) AS AREA
FROM geometries
WHERE name LIKE 'Polygon%';

注意：带孔的多边形的面积是外壳的面积（10x10 正方形）减去孔的面积（1x1 正方形）。

集合

有四种集合类型，将多个简单的几何图形组合成集合。

• MultiPoint , 点的集合；
• MultiLineString，线串的集合；
• MultiPolygon，多边形的集合；
• GeometryCollection，任何几何的异构集合（包括其他集合）；

集合是另一个在 GIS 软件中出现的概念，而不是在通用图形软件中。它们对于将现实世界对象直接建模为空间对象很有用。如，如何为用通行权分割地段建模？作为MultiPolygon，在路权的两侧都有一部分。

示例集合包含一个多边形和一个点：

SELECT 
	name,
	ST_AsText(geom)
FROM geometries
WHERE name = 'Collection';

用于处理集合的特定空间一些函数

• ST_NumGeometries(geometry)：返回集合中的零件数；
• ST_GeometryN(geometry,n)：返回指定部分；
• ST_Area(geometry)：返回所有多边形部分的总面积；
• ST_Length(geometry)：返回所有线性部分的总长度；

输入/输出

在数据库中，几何图形以仅由 PostGIS 程序使用的格式存储在磁盘上。为了让外部程序插入和检索有用的几何图形，需要将它们转换为其他应用程序可以理解的格式。

格式	输入	输出
WKT（Well-known text）	ST_GeomFromText(text, srid)	ST_AsText(geometry) ST_AsEWKT(geometry)
WKB（Well-known binary）	ST_GeomFromWKB(bytea)	ST_AsBinary(geometry) ST_AsEWKB(geometry)
GML（地理标记语言）	ST_GeomFromGML(text)	ST_AsGML(geometry)
KML（Keyhole Mark-up Language）	ST_GeomFromKML(text)	ST_AsKML(几何)
GeoJSON		ST_AsGeoJSON(geometry)
SVG		ST_AsSVG(geometry)

注意，除了具有几何表示的文本参数外，我们还有一个数字参数提供几何的SRID。

以下 SQL 查询显示了WKB表示的示例（需要调用encode()将二进制输出转换为 ASCII 格式以进行打印）：

SELECT encode(
	ST_AsBinary(
		ST_GeometryFromText('LINESTRING(0 0,1 0)')
	)
, 'hex');

出于本次研讨的目的，将继续使用 WKT 以确保您能够阅读和理解正在查看的几何图形。对于大多数实际过程，例如在 GIS 应用程序中查看数据、将数据传输到 Web 服务或远程处理数据，WKB 是首选格式。

由于 WKT 和 WKB 是在 SFSQL规范中定义，它们不处理 3 维或 4 维几何。对于这些情况，PostGIS 定义了EWKT和EWKB格式，提供了与 WKT 和 WKB 相同的格式化功能，并增加了维度。

以下是 WKT 中的 3D 线示例

SELECT ST_AsText(ST_GeometryFromText('LINESTRING(0 0 0,1 0 0,1 1 2)'));

注意，文本表示会发生变化！因为 PostGIS 的文本输入例程在消耗内容方面是自由的。

在输出端，ST_AsText函数是保守的，只发出 ISO 标准的WKT。

除了ST_GeometryFromText函数之外，还有许多其他方法可以从WKT或类似格式的输入创建几何图形：

-- Using ST_GeomFromText with the SRID parameter
SELECT ST_GeomFromText('POINT(2 2)',4326);

-- Using ST_GeomFromText without the SRID parameter
SELECT ST_SetSRID(ST_GeomFromText('POINT(2 2)'),4326);

-- Using a ST_Make* function
SELECT ST_SetSRID(ST_MakePoint(2, 2), 4326);

-- Using PostgreSQL casting syntax and ISO WKT
SELECT ST_SetSRID('POINT(2 2)'::geometry, 4326);

-- Using PostgreSQL casting syntax and extended WKT
SELECT 'SRID=4326;POINT(2 2)'::geometry;

一个使用 GML 并输出 JSON 的示例

SELECT ST_AsGeoJSON(ST_GeomFromGML('<gml:Point><gml:coordinates>1,1</gml:coordinates></gml:Point>'));

文本构建

到目前为止，看到的WKT字符串是文本类型的，我们一直在使用诸如**ST_GeomFromText()**之类的 PostGIS 函数将它们转换为几何类型。

PostgreSQL 包含一个允许数据从一种类型转换为另一种类型的短格式语法，即转换语法olddata::newtype。

例如，此 SQL 将双精度数转换为文本字符串

SELECT 0.9::text;

将WKT字符串转换为几何

SELECT 'POINT(0 0)'::geometry;

注意：强制转换时，除非指定 SRID，否则将获得具有未知 SRID 的几何图形。您可以使用EWKT形式指定 SRID，其中包括前面的 SRID 块。

SELECT 'SRID=4326;POINT(0 0)'::geometry;

函数列表

ST_Area：如果它是多边形或多多边形，则返回表面的面积。对于“几何”类型，区域以 SRID 为单位。“地理”面积以平方米为单位。
ST_AsText：返回没有 SRID 元数据的几何/地理的众所周知的文本 (WKT) 表示。
ST_AsBinary：返回没有 SRID 元数据的几何/地理的众所周知的二进制 (WKB) 表示。
ST_EndPoint：将 LINESTRING 几何的最后一个点作为 POINT 返回。
ST_AsEWKB：返回具有 SRID 元数据的几何图形的众所周知的二进制 (WKB) 表示。
ST_AsEWKT：返回具有 SRID 元数据的几何图形的众所周知的文本 (WKT) 表示。
ST_AsGeoJSON：将几何图形作为 GeoJSON 元素返回。
ST_AsGML：将几何图形作为 GML 版本 2 或 3 元素返回。
ST_AsKML：将几何图形作为 KML 元素返回。几个变种。默认版本=2，默认精度=15。
ST_AsSVG：在给定几何或地理对象的 SVG 路径数据中返回几何。
ST_ExteriorRing：返回表示 POLYGON 几何图形外环的线串。如果几何不是多边形，则返回 NULL。不适用于 MULTIPOLYGON
ST_GeometryN：如果几何是 GEOMETRYCOLLECTION、MULTIPOINT、MULTILINESTRING、MULTICURVE 或 MULTIPOLYGON，则返回从 1 开始的第 N 个几何。否则，返回 NULL。
ST_GeomFromGML：将几何的 GML 表示作为输入并输出 PostGIS 几何对象。
ST_GeomFromKML：将几何的 KML 表示作为输入并输出 PostGIS 几何对象
ST_GeomFromText：从众所周知的文本表示 (WKT) 返回指定的 ST_Geometry 值。
ST_GeomFromWKB：从众所周知的二进制几何表示 (WKB) 和可选的 SRID 创建几何实例。
ST_GeometryType：返回 ST_Geometry 值的几何类型。
ST_InteriorRingN：返回多边形几何的第 N 个内部线串环。如果几何不是多边形或给定的 N 超出范围，则返回 NULL。
ST_Length：如果它是线串或多线串，则返回几何图形的 2d 长度。几何以空间参考为单位，地理以米为单位（默认球体）
ST_NDims：以小整数形式返回几何的坐标维度。值为：2,3 或 4。
ST_NPoints：返回几何中的点（顶点）数。
ST_NRings：如果几何图形是多边形或多多边形，则返回环数。
ST_NumGeometries：如果几何是 GEOMETRYCOLLECTION（或 MULTI*），则返回几何数，否则返回 NULL。
ST_Perimeter：返回 ST_Surface 或 ST_MultiSurface 值边界的长度测量值。（多边形，多多边形）
ST_SRID：返回 ST_Geometry 的空间参考标识符，如 spatial_ref_sys 表中定义。
ST_StartPoint：将 LINESTRING 几何的第一个点作为 POINT 返回。
ST_X：返回点的 X 坐标，如果不可用，则返回 NULL。输入必须是一个点。
ST_Y：返回点的 Y 坐标，如果不可用，则返回 NULL。输入必须是一个点。

练习

📌练习：West Village 社区的面积是多少？

SELECT 
	ST_Area(geom)
FROM nyc_neighborhoods
WHERE name = 'West Village';

面积以平方米为单位。要获得公顷面积单位，请除以 10000。

📌练习：Pelham 街的几何类型及长度？

SELECT
	ST_GeometryType(geom),
	ST_Length(geom),
	ST_Transform(geom::geometry, 4326) as geom
FROM nyc_streets
WHERE name = 'Pelham St';

📌练习：Broad St 地铁站的 GeoJSON 表示是什么？

SELECT
	ST_AsGeoJSON(geom)
FROM nyc_subway_stations
WHERE name = 'Broad St';

📌练习：纽约市街道的总长度是多少公里？

SELECT
	Sum(ST_Length(geom)) / 1000
FROM nyc_streets;

📌练习：曼哈顿的面积是多少英亩？

SELECT 
	Sum(ST_Area(geom)) / 4047
FROM nyc_neighborhoods
WHERE boroname = 'Manhattan';

或是：

SELECT 
	Sum(ST_Area(geom)) / 4047
FROM nyc_census_blocks
WHERE boroname = 'Manhattan';

📌练习：最西边的地铁站在哪？

SELECT 
	name,
	ST_X(geom), 
	ST_Y(geom),
	ST_Transform(geom::geometry, 4326) as geom
FROM nyc_subway_stations
ORDER BY ST_X(geom)
LIMIT 1;

📌练习：Columbus Cir 有多长？

SELECT 
	ST_Length(geom)
FROM nyc_streets
WHERE name = 'Columbus Cir';

📌练习：纽约市街道的长度按类型分组？

SELECT 
	type, 
	Sum(ST_Length(geom)) AS length
FROM nyc_streets
GROUP BY type
ORDER BY length DESC;

该子句对结果按长度降序排序。

空间关系

空间数据库可以比较存储几何图形的空间关系。

相等

ST_Equals(geometry A, geometry B)：比较两个几何是否相等。

如果相同类型的两个几何图形具有相同的 x,y 坐标值，即如果第二个形状与第一个形状相等（相同）。

首先，从nyc_subway_stations表中筛选名为 Broad St 的点。

select
	name,
	encode(ST_AsBinary(geom), 'hex') as Binary
FROM nyc_subway_stations
WHERE name = 'Broad St';

点的表示不是很容易理解（0101000020266900000EEBD4CF27CF2141BC17D69516315141），但它是坐标值的精确表示。对于像相等这样的测试，使用精确的坐标是必要的。

然后，将几何表示重新插入ST_Equals测试

SELECT
	name
FROM nyc_subway_stations
WHERE ST_Equals(geom,'0101000020266900000EEBD4CF27CF2141BC17D69516315141');

相交

ST_Intersects(geometry A, geometry B)	ST_Disjoint(geometry A , geometry B)	ST_Crosses(geometry A, geometry B)	ST_Overlaps(geometry A, geometry B)
如果两个形状有任何公共空间，即它们的边界或内部相交。	如果两个几何图形不相交，则它们不相交。	如果相交产生的几何维度比两个原几何的最大维度小，并且相交集在两个原几何内部。	比较两个具有相同维度的几何，如果它们的交集几何不同但具有相同维度。
	事实上，测试not intersects通常比测试disjoint更有效，因为相交测试用到空间索引，而不相交测试则没有。

以 Broad Street 地铁站为例，使用ST_Intersects函数确定其邻域

SELECT 
	name,
	ST_AsText(geom)
FROM nyc_subway_stations
WHERE name = 'Broad St';

SELECT
	name,
	boroname
FROM nyc_neighborhoods
WHERE ST_Intersects(geom, ST_GeomFromText('POINT(583571 4506714)', 26918));

接触

ST_Touches(geometry A, geometry B)：比较两个几何是否内部不相交。且边界处接触。

包含

ST_Contains()（包含）、ST_Within()(被包含)比较一个几何图形是否完全包含在另一个几何图形内。

距离

ST_Distance(geometry A, geometry B)：计算两个几何之间的最短距离。

SELECT ST_Distance(
  ST_GeometryFromText('POINT(0 5)'),
  ST_GeometryFromText('LINESTRING(-2 2, 2 2)')
);

ST_DWithin(geometry A, geometry B, double precision)：测试两个对象是否在指定的距离内，不必计算实际缓冲区，只需测试距离关系。

找到 Broad Street 地铁站附近（10 米内）的街道

SELECT
	name
FROM nyc_streets
WHERE ST_DWithin(geom,ST_GeomFromText('POINT(583571 4506714)', 26918), 10);

落图验证

(
	select '-' as name, 
	ST_Transform(ST_GeomFromText('POINT(583571 4506714)', 26918)::geometry, 4326) as geom
)
union all
(
SELECT
	name,
	ST_Transform(geom::geometry, 4326) as geom
FROM nyc_streets
WHERE ST_DWithin(geom,ST_GeomFromText('POINT(583571 4506714)', 26918), 10)
);

函数列表

• ST_Contains(geometry A, geometry B)：当且仅当 B 的任何点都不位于 A 的外部，并且 B 的内部至少有一个点位于 A 的内部时，才返回 true。
• ST_Crosses(geometry A, geometry B)：如果提供的几何图形有一些但不是所有的公共内部点，则返回 TRUE。
• ST_Disjoint(geometry A , geometry B)：如果几何不“空间相交” - 如果它们不共享任何空间，则返回 TRUE。
• ST_Distance(geometry A, geometry B)：以投影单位返回两个几何之间的二维笛卡尔最小距离（基于空间参考）。
• ST_DWithin(geometry A, geometry B, radius)：如果几何在彼此之间的指定距离（半径）内，则返回 true。
• ST_Equals(geometry A, geometry B)：如果给定的几何表示相同的几何，则返回 true。方向性被忽略。
• ST_Intersects(geometry A, geometry B)：如果几何/地理“空间相交” - （共享空间的任何部分）返回 TRUE，否则返回 FALSE（它们是不相交的）。
• ST_Overlaps(geometry A, geometry B)：如果几何共享空间，具有相同的维度，但不完全相互包含，则返回 TRUE。
• ST_Touches(geometry A, geometry B)：如果几何至少有一个共同点，但它们的内部不相交，则返回 TRUE。
• ST_Within(geometry A , geometry B)：如果几何 A 完全在几何 B 内，则返回 true

练习

📌练习：获取 大西洋公地 街道的几何图形

SELECT
	ST_AsText(geom)
FROM nyc_streets
WHERE name = 'Atlantic Commons';

📌练习：大西洋公地 在哪个社区和行政区？

SELECT
	name,
	boroname
FROM nyc_neighborhoods
WHERE ST_Intersects(geom, ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918));

为什么从MULTILINESTRING变成LINESTRING？
在空间上，它们描述了相同的形状，因此从单项多几何体到单例体可以节省一些击键次数。
更重要的是，我们还对坐标进行了四舍五入以使它们更易于阅读，这实际上改变了结果：我们无法使用 ST_Touches() 谓词找出哪些道路连接大西洋公地，因为坐标不再完全相同。

落图验证

(
	select '-' as name, 
	'-' as boroname,
	ST_Transform(ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918)::geometry, 4326) as geom
)
union all
(
SELECT
	name,
	boroname,
	ST_Transform(geom::geometry, 4326) as geom
FROM nyc_neighborhoods
WHERE ST_Intersects(geom, ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918))
);

📌练习：大西洋公地与哪些街道相连？

SELECT 
	name
FROM nyc_streets
WHERE ST_DWithin(geom, ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918), 1);

落图验证

(
	select 'Atlantic Commons' as name, 
	ST_Transform(ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918)::geometry, 4326) as geom
)
union all
(
SELECT 
	name,
	ST_Transform(geom::geometry, 4326) as geom
FROM nyc_streets
WHERE ST_DWithin(geom, ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918), 1)
);

📌练习：大约有多少人住在大西洋公地（50 米范围内）？

SELECT
	Sum(popn_total)
FROM nyc_census_blocks
WHERE ST_DWithin(geom, ST_GeomFromText('LINESTRING(586782 4504202,586864 4504216)', 26918), 50);

空间连接

空间连接是将不同表的数据通过空间关系组合起来。

📌示例：查询有关地铁站及其所在社区的信息

SELECT
  subways.name AS subway_name,
  neighborhoods.name AS neighborhood_name,
  neighborhoods.boroname AS borough
FROM nyc_neighborhoods AS neighborhoods
JOIN nyc_subway_stations AS subways
ON ST_Contains(neighborhoods.geom, subways.geom)
WHERE subways.name = 'Broad St';

 subway_name | neighborhood_name  |  borough
-------------+--------------------+-----------
 Broad St    | Financial District | Manhattan

连接 和 聚合

使用 JOIN与GROUP BY` 组合分析。

📌示例：查询曼哈顿社区的人口构成

SELECT
  neighborhoods.name AS neighborhood_name,
  Sum(census.popn_total) AS population,
  100.0 * Sum(census.popn_white) / Sum(census.popn_total) AS white_pct,
  100.0 * Sum(census.popn_black) / Sum(census.popn_total) AS black_pct
FROM nyc_neighborhoods AS neighborhoods
JOIN nyc_census_blocks AS census
ON ST_Intersects(neighborhoods.geom, census.geom)
WHERE neighborhoods.boroname = 'Manhattan'
GROUP BY neighborhoods.name
ORDER BY white_pct DESC;

要统计的指定半径内的项目查询，可以使用距离测试测量函数作为连接键。

📌示例：使用距离查询来探索纽约的种族地理。

首先，了解一下这座城市的基本种族构成

SELECT
  100.0 * Sum(popn_white) / Sum(popn_total) AS white_pct,
  100.0 * Sum(popn_black) / Sum(popn_total) AS black_pct,
  Sum(popn_total) AS popn_total
FROM nyc_census_blocks;

在纽约的 800 万人口中，大约 44% 被记录为白人，26% 被记录为黑人。

哈莱姆区的非裔人口在曼哈顿遥遥领先（80.5%），杜克的A-train也是这样吗？

首先注意 nyc_subway_stations 表的 routes 字段内容就是我们感兴趣的找A-train。里面的值有点复杂。

SELECT
	-- 去除结果中重复行
	DISTINCT routes 
FROM nyc_subway_stations;

因此，要找到 A 列车，我们将需要routes中包含A的任何行。

SELECT
  DISTINCT routes
FROM nyc_subway_stations AS subways
WHERE strpos(subways.routes, 'A') > 0;

📌示例：统计 A-train 200米以内的种族构成

SELECT
  100.0 * Sum(popn_white) / Sum(popn_total) AS white_pct,
  100.0 * Sum(popn_black) / Sum(popn_total) AS black_pct,
  Sum(popn_total) AS popn_total
FROM nyc_census_blocks AS census
JOIN nyc_subway_stations AS subways
ON ST_DWithin(census.geom, subways.geom, 200)
WHERE strpos(subways.routes,'A') > 0;

A-train 沿线的种族构成与整个纽约市的构成并没有根本不同。

高级连接

📌示例：按线路统计人口构成

为了回答这个问题，将在查询中添加另一个连接，以便我们可以同时计算多条地铁线路的组成。

CREATE TABLE subway_lines ( route char(1) );
INSERT INTO subway_lines (route) VALUES
  ('A'),('B'),('C'),('D'),('E'),('F'),('G'),
  ('J'),('L'),('M'),('N'),('Q'),('R'),('S'),
  ('Z'),('1'),('2'),('3'),('4'),('5'),('6'),
  ('7');

将地铁线路表加入到我们的原始查询中。

SELECT
  lines.route,
  100.0 * Sum(popn_white) / Sum(popn_total) AS white_pct,
  100.0 * Sum(popn_black) / Sum(popn_total) AS black_pct,
  Sum(popn_total) AS popn_total
FROM nyc_census_blocks AS census
JOIN nyc_subway_stations AS subways
ON ST_DWithin(census.geom, subways.geom, 200)
JOIN subway_lines AS lines
ON strpos(subways.routes, lines.route) > 0
GROUP BY lines.route
ORDER BY black_pct DESC;

函数列表

• ST_Contains(geometry A, geometry B)：当且仅当 B 的任何点都不位于 A 的外部，并且 B 的内部至少有一个点位于 A 的内部时，才返回 true。
• ST_DWithin(geometry A, geometry B, radius)：如果几何在彼此之间的指定距离内，则返回 true。
• ST_Intersects(geometry A, geometry B)：如果几何/地理“空间相交” - （共享空间的任何部分）返回 TRUE，否则返回 FALSE（它们是不相交的）。
• round(v numeric, s integer)：四舍五入到 s 小数位的 PostgreSQL 数学函数
• strpos(string, substring)：PostgreSQL 字符串函数，返回指定子字符串的整数位置。
• sum(expression)：PostgreSQL 聚合函数，返回一组记录中记录的总和。

练习

📌练习：Little Italy 在哪个地铁站？在哪条地铁线路上？

SELECT s.name, s.routes
FROM nyc_subway_stations AS s
JOIN nyc_neighborhoods AS n
ON ST_Contains(n.geom, s.geom)
WHERE n.name = 'Little Italy';

📌练习：6号线服务的社区有哪些？（提示：表routes中的列nyc_subway_stations具有’B，D，6，V’和’C，6’之类的值）

SELECT DISTINCT n.name, n.boroname
FROM nyc_subway_stations AS s
JOIN nyc_neighborhoods AS n
ON ST_Contains(n.geom, s.geom)
WHERE strpos(s.routes,'6') > 0;

我们使用DISTINCT关键字从我们的结果集中删除重复值，其中附近有多个地铁站。

📌练习：911 之后，Battery Park 社区连续几天被禁止进入，有多少人需要疏散？

SELECT Sum(popn_total)
FROM nyc_neighborhoods AS n
JOIN nyc_census_blocks AS c
ON ST_Intersects(n.geom, c.geom)
WHERE n.name = 'Battery Park';

📌练习：哪个社区的人口密度最高（人/平方公里）？

SELECT
	n.name,
	Sum(c.popn_total) / (ST_Area(n.geom) / 1000000.0) AS popn_per_sqkm
FROM nyc_census_blocks AS c
JOIN nyc_neighborhoods AS n
ON ST_Intersects(c.geom, n.geom)
GROUP BY n.name, n.geom
ORDER BY 2 DESC;

空间索引

空间索引将数据组织到搜索树中来加速搜索，便于快速遍历来找到特定记录，是空间数据库的三个关键特征之一。

若没有索引，查询的开销将非常巨大，如连接查询两张包含 10,000 条记录的表，不用索引需比较 100,000,000 次，而使用索引仅需比较 20,000 次。

以nyc_census_blocks表举例说明，数据文件加载时已经建立了对应的空间索引，为了展示索引的功效，首先删除该表的索引。

1️⃣ 删除索引

DROP INDEX nyc_census_blocks_geom_idx;

2️⃣ 查询每个人口普查区块，并过滤以B开头的地铁站。

SELECT count(blocks.blkid)
FROM nyc_census_blocks blocks
JOIN nyc_subway_stations subways
ON ST_Contains(blocks.geom, subways.geom)
WHERE subways.name LIKE 'B%';

耗时：

3️⃣ 添加空间索引

CREATE INDEX nyc_census_blocks_geom_idx
ON nyc_census_blocks
USING GIST (geom);

4️⃣ 再次查询

SELECT count(blocks.blkid)
FROM nyc_census_blocks blocks
JOIN nyc_subway_stations subways
ON ST_Contains(blocks.geom, subways.geom)
WHERE subways.name LIKE 'B%';

耗时：

工作机制

空间索引不是通过被索引的列值来创建层次树，而是索引空间特征边框。

图中与黄星相交的线数为1，即红线;但与黄色框相交的特征边界框是两个，红色和蓝色。

需要计算哪些线与黄星相交？分以下步骤完成：

1. 计算哪些框与黄色框相交；
2. 精确计算哪些线与黄色框相交；

对于大型表，先评估近似索引，然后执行精确计算，从根本上减少了计算量。

PostGIS 和 Oracle Spatial 共享相同的 R-Tree 空间索引结构。R-Trees 将数据分解为矩形、子矩形、子子矩形等；是一种自调整的索引结构，可以自动处理可变数据密度、不同数量的对象重叠和对象大小。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-19AQIQf9-1646290507003)(https://postgis.net/workshops/postgis-intro/_images/index-01.png)]

索引函数

在空间索引可用的前提下，只有一部分函数会自动使用空间索引

• ST_Intersects
• ST_Contains
• ST_Within
• ST_DWithin
• ST_ContainsProperly
• ST_CoveredBy
• ST_Covers
• ST_Overlaps
• ST_Crosses
• ST_DFullyWithin
• ST_3D ST_3DIntersects
• ST_3DDWithin
• ST_3DDFullyWithin
• ST_LineCrossingDirection
• ST_OrderingEquals
• ST_Equals

仅索引查询

PostGIS 中的大多数常用函数（如：ST_Contains、ST_Intersects、ST_DWithin等）都自动包含索引过滤器。但是某些函数（如：ST_Relate）不包含索引过滤器。

可以通过使用&&运算符（表示边界框重叠或接触），对边界框搜索

📌示例：将 West Village 人口的仅索引查询与更精确的查询进行比较。

使用&&的索引查询方式：

SELECT
	Sum(popn_total)
FROM nyc_neighborhoods neighborhoods
JOIN nyc_census_blocks blocks
ON neighborhoods.geom && blocks.geom
WHERE neighborhoods.name = 'West Village';

49821

使用ST_Intersects的精确查询方式

SELECT
	Sum(popn_total)
FROM nyc_neighborhoods neighborhoods
JOIN nyc_census_blocks blocks
ON ST_Intersects(neighborhoods.geom, blocks.geom)
WHERE neighborhoods.name = 'West Village';

26718

注意：第一个查询统计了边界框与邻域边界框相交的每个块；第二个查询仅统计与邻域本身相交的那些块。

分析

PostgreSQL 查询计划器智能地选择何时使用或不使用索引来评估查询。索引搜索并不总是最快的，如果搜索要返回表中的每条记录，那么遍历索引树以获取每条记录实际上会比从头开始顺序读取整个表要慢。

明确查询矩形的大小并不能确定查询返回的数据量，如下图所示，红色方块很小，但会返回比蓝色方块多得多的记录。

PostgreSQL 数据库为了直到自己正在处理什么数据，它会定期收集、保留有关每个索引表列中数据分布的统计信息。

如果在短时间内大幅更改表里的内容，统计信息不一定能及时更新；建议在批量加载、修改、删除数据后，使用ANALYZE命令来同步表的统计信息。

ANALYZE nyc_census_blocks;

清理

维护索引不仅仅创建即可，当对表发出大量INSERT、UPDATE、DELETE时，配合VACUUM命令可回收表页留下的未使用空间，使得数据库能够高效运行。

PostgreSQL 默认提供了Autovacuum(自动清理)功能，可通过活动级别设定合理的间隔对表进行清理（恢复空间）和分析（更新统计信息）。虽然这对于高度事务性数据库至关重要，但不建议在添加索引、批量加载数据后等待 autovacuum 运行。每当执行大批量更新时，您应该手动运行VACUUM.

可以按需单独执行清理和分析数据库，通过VACUUM命令不会更新数据库统计信息；通过ANALYZE命令也不会恢复未使用的表行。
这两个命令都可以针对整个数据库、单个表或单个列运行。

VACUUM ANALYZE nyc_census_blocks;

函数列表

• geometry_a && geometry_b：如果 A 的边界框与 B 的边界框重叠，则返回 TRUE；
• geometry_a = geometry_b：如果 A 的边界框与 B 的边界框相同，则返回 TRUE；
• ST_Intersects(geometry_a, geometry_b)：如果几何/地理 空间相交 - （共享空间的任何部分）返回 TRUE。

投影数据

地球不是平的，没有简单的方法可以把它放在平面纸质地图（或电脑屏幕）上，所以人们想出了各种巧妙的解决方案，各有利弊。

• 一些投影保留区域，因此所有对象都具有彼此的相对大小；
• 其他投影保留角度，如墨卡托投影；
• 一些预测试图找到一个很好的中间混合，在几个参数上只有很小的失真。

所有投影的共同点是它们将（球形）世界转换为平面笛卡尔坐标系，选择哪种投影取决于您将如何使用数据。

PostGIS 可以通过ST_Transform(geometry, srid)函数来更改数据投影内置支持；通过ST_SRID( **geometry)、ST_SetSRID(geometry, srid)函数来管理空间参考标识符。

📌示例：使用ST_SRID函数来确认数据的SRID

SELECT ST_SRID(geom)
FROM nyc_streets LIMIT 1;

📌示例：SRID(空间参考标识符)的定义在spatial_ref_sys表中

SELECT * 
FROM spatial_ref_sys 
WHERE srid = 26918;

PostGIS 重投影引擎将尝试从spatial_ref_sys表中找到最佳投影

• 如果 proj 尝试在其内部目录中找到有效的authority name和authority srid；
• 如果 proj 尝试从srtext中解析；
• 最后 proj 尝试从proj4text中解析。

注意

• 这些冗余设计意味着在 PostGIS 中创建新投影所需的只是一个有效的srtext字符串或proj4text字符串；
• 默认情况下，所有常见的 authority name、authority srid 对都已加载到表中。
• 在创建自定义投影时，请填写srtext列，该列被GeoServer、QGIS和FME 等外部程序使用。

比较数据

PostGIS 函数在执行空间操作时，需要所操作空间数据属于同一个SRID。

坐标 和 SRID 定义了地球上的位置：

• 如果没有 SRID，坐标只是一个抽象概念;
• 笛卡尔坐标平面被定义为放置在地球表面的平面坐标系。

📌示例：比较不同 SRID 的几何图形，会抛出错误信息

SELECT ST_Equals(
	ST_GeomFromText('POINT(0 0)', 4490),
	ST_GeomFromText('POINT(0 0)', 26918)
);

注意：
空间索引是使用存储几何的 SRID 构建的；
在未指定 SRID 的情况下加载数据或创建新几何，则 SRID 值将为 0；
如果在不同的 SRID 中进行比较，则（通常）不使用空间索引；
最佳做法是为数据库中的所有表选择同一个SRID，只有在向外部应用程序读取或写入数据时才使用转换功能。

转换数据

📌示例：获取 SRID 为 26918 的 proj4 定义

SELECT srtext 
FROM spatial_ref_sys 
WHERE srid = 26918;

可以看到工作投影是 UTM zone 18N（通用横轴墨卡托），以米为测量单位。

PROJCS["NAD83 / UTM zone 18N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree|

📌示例：将投影中的一些数据转换为地理坐标——也称为 经度/纬度

地理坐标最常见的 SRID 是 4490，对应于 CGCS2000 椭球体上的 经度 / 纬度。
可以EPSG查看定义。

还可以从spatial_ref_sys表中提取定义：

SELECT srtext 
FROM spatial_ref_sys 
WHERE srid = 4490;

GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG|

让我们将Broad St地铁站的坐标转换为地理坐标

SELECT
	ST_AsText(ST_Transform(geom,4490))
FROM nyc_subway_stations
WHERE name = 'Broad St';

📌示例：查询geometry_columns表分配的SRID

SELECT
	srid,
	f_table_name AS name
FROM geometry_columns;

📌示例：通过ST_SetSRID设置坐标的 SRID

SELECT
    ST_AsText(
        ST_Transform(
            ST_SetSRID(geom, 26918),
        4490)
    )
FROM geometries;

函数列表

• ST_AsText：返回没有 SRID 元数据的几何/地理的众所周知的文本 (WKT) 表示。
• ST_SetSRID(geometry, srid)：将几何上的 SRID 设置为特定的整数值。
• T_SRID(geometry)：返回 ST_Geometry 的空间参考标识符，如 spatial_ref_sys 表中所定义。
• T_Transform(geometry, srid)：返回一个新几何，其坐标转换为整数参数引用的 SRID。

练习

在线资源：epsg.io

📌练习：用 UTM 18 测量纽约所有街道的长度是多少？

SELECT
	Sum(ST_Length(geom))
FROM nyc_streets;

📌练习：用 SRID 4490 测量纽约所有街道的长度是多少？

SELECT Sum(ST_Length(ST_Transform(geom,4490)))
FROM nyc_streets;

📌练习：查看 SRID 2831 的 WKT 定义

SELECT srtext
FROM spatial_ref_sys
WHERE SRID = 2831;

PROJCS["NAD83(HARN) / New York Long Island",GEOGCS["NAD83(HARN)",DATUM["NAD83_High_Accuracy_Reference_Network",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6152"]],PRIMEM["Greenwich",0,AUTHOR|

📌练习：查看有多少条街道穿过第 74 条子午线？

SELECT 
	Count(*)
FROM nyc_streets
WHERE ST_Intersects(ST_Transform(geom, 4326),'SRID=4326;LINESTRING(-74 20, -74 60)');

📌练习：
第 74 条子午线是一种奇特的说法，即地理上 X 值为 -74 的垂直线。可以构建这样一条线，然后将其与街道进行比较，也可以投影到地理上。

将线投影到 UTM 并在那里进行比较将返回略有不同的答案。要获得相同的答案，您需要对其进行分段，使其具有更多点，然后再进行转换。

SELECT Count(*)
FROM nyc_streets
WHERE ST_Intersects(geom, ST_Transform(ST_Segmentize('SRID=4326;LINESTRING(-74 20, -74 60)'::geometry,0.001), 26918));

地理坐标

地理（纬度/经度）坐标是参考子午线旋转角度（经度）和赤道的角度（纬度）的球坐标。

地理坐标与墨卡托、UTM 或 Stateplane 中的坐标不同

• 不表示与平面上绘制的原点的线性距离;
• 不是笛卡尔坐标。

注意
球坐标测量角距离，单位为度，测量的距离、长度和面积是没有意义的。
索引和检测（相交、包含）的近似结果可能是错误的，当接近两极或国际日期变更线等问题区域时，点之间的距离会变大。

示例坐标：洛杉矶（POINT(-118.4079 33.9434)）、巴黎（POINT(2.3490 48.8533)）。

📌示例：使用标准笛卡尔ST_Distance(geometry, geometry)计算洛杉矶和巴黎之间的距离

SELECT ST_Distance(
	-- Los Angeles (LAX)
	'SRID=4326;POINT(-118.4079 33.9434)'::geometry,
	-- Paris (CDG)
	'SRID=4326;POINT(2.5559 49.0083)'::geometry
);

SRID 4326 的单位是度，所以答案是122度。

在球体上，一个平方度的大小变化很大，随着远离赤道而变小，测量的距离没有任何意义。

为了计算出准确的距离，需要通过geography类型测量球体上真实的距离。

不同的空间数据库有不同的处理地理的方法

• Oracle 试图通过透明地进行地理计算来掩盖差异；
• SQL Server 使用两种空间类型，STGeometry用于笛卡尔数据，STGeography用于地理；
• Informix Spatial 是 Informix 的纯笛卡尔扩展，而 Informix Geodetic 是纯地理扩展；
• PostGIS 使用两种类型，几何和地理。

📌示例：使用geography类型来测量洛杉矶和巴黎之间的距离。

SELECT ST_Distance(
	-- Los Angeles (LAX)
    'SRID=4326;POINT(-118.4079 33.9434)'::geography,
    -- Paris (CDG)
    'SRID=4326;POINT(2.5559 49.0083)'::geography
);

计算的所有返回值geography都以米为单位，所以答案是 9125km。

旧版本的 PostGIS 支持使用ST_Distance_Spheroid(point, point, measure)函数对球体进行非常基本的计算。
但是ST_Distance_Spheroid功能仅适用于点，不支持跨两极或国际日期变更线的索引。

当遇到技术从洛杉矶-巴黎的航线与冰岛最近的距离时，就需要支持非点几何计算。

在笛卡尔平面上，使用地理坐标（紫色线）会产生错误的答案，使用大圆圈路线（红线）能得出正确答案；将 LAX-CDG 航线转换为线串，并计算到冰岛某个点的距离，将得到正确答案。

SELECT ST_Distance(
  ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'), -- LAX-CDG
  ST_GeographyFromText('POINT(-22.6056 63.9850)')                        -- Iceland (KEF)
);

在 LAX-CDG 航线上，离冰岛最近的距离是 502 公里。

跨越国际日期变更线的情况下，地理坐标的处理笛卡尔方法完全失效。从洛杉矶-东京的最短大圆路线穿越太平洋，最短的笛卡尔路线横跨大西洋和印度洋。

SELECT
	ST_Distance(ST_GeometryFromText('Point(-118.4079 33.9434)'), ST_GeometryFromText('Point(139.733 35.567)')) AS geometry_distance,
	ST_Distance(ST_GeographyFromText('Point(-118.4079 33.9434)'), ST_GeographyFromText('Point(139.733 35.567)')) AS geography_distance;

使用地理

为了将几何数据加载到地理表中，首先需要将几何投影到 EPSG:4326（经度/纬度），然后将其更改为地理。

转换地理坐标的方式

• ST_Transform(geometry, srid)函数；
• Geography(geometry)函数或::geography后缀。

CREATE TABLE nyc_subway_stations_geog AS
SELECT
  ST_Transform(geom,4326)::geography AS geog,
  name,
  routes
FROM nyc_subway_stations;

地理坐标建立空间索引（与几何完全相同）

CREATE INDEX nyc_subway_stations_geog_gix
ON nyc_subway_stations_geog USING GIST (geog);

地理索引将正确处理涵盖两极或国际日期变更线的查询，而几何索引则不会。

geography 类型只有少数原生函数

• **ST_AsText(geography)**返回text
• **ST_GeographyFromText(text)**返回geography
• **ST_AsBinary(geography)**返回bytea
• **ST_GeogFromWKB(bytea)**返回geography
• **ST_AsSVG(geography)**返回text
• **ST_AsGML(geography)**返回text
• **ST_AsKML(geography)**返回text
• **ST_AsGeoJson(geography)**返回text
• **ST_Distance(geography, geography)**返回double
• **ST_DWithin(geography, geography, float8)**返回boolean
• **ST_Area(geography)**返回double
• **ST_Length(geography)**返回double
• **ST_Covers(geography, geography)**返回boolean
• **ST_CoveredBy(geography, geography)**返回boolean
• **ST_Intersects(geography, geography)**返回boolean
• **ST_Buffer(geography, float8)**返回geography
• **ST_Intersection(geography, geography)**返回geography

创建地理表

创建具有地理坐标表 与 创建几何表非常相似；但是，地理坐标在创建表时需通过GEOGRAPHY(type)指定对象类型。

例如：

CREATE TABLE airports (
    code VARCHAR(3),
    geog GEOGRAPHY(Point)
);

INSERT INTO airports VALUES ('LAX', 'POINT(-118.4079 33.9434)');
INSERT INTO airports VALUES ('CDG', 'POINT(2.5559 49.0083)');
INSERT INTO airports VALUES ('KEF', 'POINT(-22.6056 63.9850)');

在表定义中，通过GEOGRAPHY(Point)将机场数据类型指定为点。

新添加的地理字段不会注册在geometry_columns视图中；它们注册在geography_columns中。

SELECT * FROM geography_columns;

转换为几何

PostgreSQL中用于转换的语法约定是：在需要转换的值添加::typename后缀。如：2::text将数字 2 转换为文本字符串"2"；'POINT(0 0)'::geometry将点的文本表示转换为几何点。

如何从地理坐标中读取 X 坐标呢？（ST_X(point)函数仅支持几何类型）

SELECT
	code,
	ST_X(geog::geometry) AS longitude 
FROM airports;

通过在 geography 值后添加::geometry后缀，将对象转换为 SRID 为 4326 的几何体。

为什么不直接使用地理坐标

每个人都理解纬度/经度，地理坐标是最常见的，而很少有人理解 UTM 坐标。为什么不直接使用地理坐标呢？

• 支持地理类型的函数还比较少，需要花费大量时间来解决地理类型限制；
• 在球体上的计算远比笛卡尔复杂。如，距离计算，笛卡尔公式(Pythagoras)涉及到对 sqrt() 的一次调用。距离计算，球形公式(havesine)包括两次 sqrt() 调用、一次 arctan()调用、四次 sin () 调用和两次 cos ()调用。三角函数计算非常复杂，而且球面计算涉过多。

密集数据（包含在州/市/县内），使用笛卡尔投影的几何(geometry)类型;
分散数据（覆盖世界大部分地区），使用地理(geography)类型。

函数列表

• ST_Distance(geometry, geometry)：对于几何类型返回投影单位中两个几何之间的二维笛卡尔最小距离（基于空间参考）。对于地理类型，默认返回两个地理区域之间的球体最小距离（以米为单位）。
• ST_GeographyFromText(text)：从 Well-Known Text 表示或扩展 (WKT) 返回指定的地理值。
• ST_Transform(geometry, srid)：返回一个新几何，其坐标转换为整数参数引用的 SRID。
• ST_X(point)：返回点的 X 坐标，如果不可用，则返回 NULL。输入必须是一个点。
• ST_Azimuth(geography_A, geography_B)：以弧度返回从 A 到 B 的方向。
• ST_DWithin(geography_A, geography_B, R)：如果 A 在 B 的 R 米内，则返回 true。

练习

📌练习：纽约（POINT(-74.0064 40.7142)）离西雅图（POINT(-122.3331 47.6097)）多远？

SELECT ST_Distance(
    'POINT(-74.0064 40.7142)'::geography,
    'POINT(-122.3331 47.6097)'::geography
);

📌练习：在球体上计算，纽约所有街道的总长度是多少？

SELECT 
    Sum(
        ST_Length(
            Geography(
                ST_Transform(geom, 4326)
            )
        )
    )
FROM nyc_streets;

📌练习：POINT(1 2.0001) 是否与地理或几何中的 POLYGON((0 0, 0 2, 2 2, 2 0, 0 0)) 相交？

SELECT
	ST_Intersects('POINT(1 2.0001)'::geography, 'POLYGON((0 0,0 2,2 2,2 0,0 0))'::geography) as geography_intersects,
	ST_Intersects('POINT(1 2.0001)'::geometry, 'POLYGON((0 0,0 2,2 2,2 0,0 0))'::geometry) as geometry_intersects;

• 正方形的上边缘在几何上是一条直线，从点的下方经过，正方形不包含点；
• 正方形的上边缘在地理上是一个大圆，在该点上方经过，该正方形确实包含该点。

几何构造函数

几何构造函数输入几何计算处理后输出几何。

ST_Centroid / ST_PointOnSurface

空间查询时经常需要用点表示多边形。

这对于空间连接很有用，因为在两个多边形层上使用ST_Intersects(geometry,geometry)通常会导致重复计算：边界上的多边形将与两侧的对象相交；用一个点替换它会迫使它位于一侧或另一侧，而不是两者。

• ST_Centroid(geometry)：返回几何质心点，注意返回的点不一定在本身上；
• ST_PointOnSurface(geometry)：返回内心点，作为计算空间连接的代理点使用。

📌示例：比较凹几何的质心和内心的位置

SELECT
	ST_Intersects(geom, ST_Centroid(geom)) AS centroid_inside,
	ST_Intersects(geom, ST_PointOnSurface(geom)) AS pos_inside
FROM (
	VALUES('POLYGON ((30 0, 30 10, 10 10, 10 40, 30 40, 30 50, 0 50, 0 0, 0 0, 30 0))'::geometry)
) AS t(geom);

ST_Buffer

ST_Buffer(geometry,distance)按给定的缓冲距离对指定的几何数据执行缓冲操作，将获得一个新的多边形。

📌示例：计算 nyc_census_blocks表中 自由岛 普查区块周围建立一个500 米的海上交通区

-- 用自由岛 500m 缓冲区制作新表
CREATE TABLE liberty_island_zone AS
SELECT ST_Buffer(geom,500)::geometry(Polygon,26918) AS geom
FROM nyc_census_blocks
WHERE blkid = '360610001001001';

ST_Buffer函数还接受负距离并在多边形输入中构建内接多边形。对于线和点，你只会得到一个空的回报。

ST_Intersection

通过计算两个多边形的叠加交集来创建新的覆盖范围，ST_Intersection(geometry A, geometry B)函数返回两个参数共有的空间区域（或线或点）,如果参数不相交，则该函数返回一个空几何。

SELECT ST_AsText(
    ST_Intersection(
        ST_Buffer('POINT(0 0)', 2),
        ST_Buffer('POINT(3 0)', 2)
    )
);

ST_Union

通过计算两个多边形的合并交集来创建新的覆盖范围，ST_Union(geometry, geometry) / ST_Union([geometry])函数接受输入并删除公共线路。

SELECT ST_AsText(
    ST_Union(
        ST_Buffer('POINT(0 0)', 2),
        ST_Buffer('POINT(3 0)', 2)
    )
);

作为ST_Union聚合的示例，请考虑我们的nyc_census_blocks表。

人口普查地理是精心构建的，因此可以从较小的地理中建立更大的地理。

📌示例：通过合并每个县块来创建县地图，可以通过合并所有具有相同前 5 位数字的几何图形来创建县地图blkid。

-- 通过合并人口普查块创建 nyc_census_counties 表
CREATE TABLE nyc_census_counties AS
SELECT
	SubStr(blkid,1,5) AS countyid,
	ST_Union(geom)::Geometry(MultiPolygon,26918) AS geom
FROM nyc_census_blocks
GROUP BY countyid;

📌验证：区域测试可以确认联合操作没有丢失任何几何图形。

首先，我们计算每个人口普查区块的面积，并将按人口普查县 id 分组的区域相加。

SELECT 
	SubStr(blkid,1,5) AS countyid,
	Sum(ST_Area(geom)) AS area
FROM nyc_census_blocks
GROUP BY countyid
ORDER BY countyid;

然后，从县表中计算每个新县多边形的面积：

SELECT
	countyid,
	ST_Area(geom) AS area
FROM nyc_census_counties
ORDER BY countyid;

一样的答案！已经从人口普查区数据成功地构建了一个纽约市县表。

函数列表

• ST_Centroid(geometry)：返回一个点几何，表示输入几何的质心。
• ST_PointOnSurface(geometry)：返回保证在输入几何内部的点几何。
• ST_Buffer(geometry, distance)：对于几何：返回一个几何，它表示与此几何的距离小于或等于距离的所有点。计算在此几何的空间参考系中。对于地理：使用平面变换包装器。
• ST_Intersection(geometry A, geometry B)：返回一个几何图形，表示 geomA 和 geomB 的共享部分。地理实现对几何进行转换以进行交叉，然后转换回 WGS84。
• ST_Union()：返回表示几何的点集联合的几何。
• ST_AsText(text)：返回没有 SRID 元数据的几何/地理的众所周知的文本 (WKT) 表示。
• substring(string [from int] [for int]) : PostgreSQL 字符串函数，用于提取匹配 SQL 正则表达式的子字符串。
• sum(expression)：PostgreSQL 聚合函数，返回一组记录中记录的总和。

练习

📌练习：有多少人口普查区块不包含自己的质心？

SELECT Count(*)
FROM nyc_census_blocks
WHERE NOT ST_Contains(geom, ST_Centroid(geom));

📌练习：将所有人口普查块合并为一个输出。它是什么样的几何图形？它有几个部分？

CREATE TABLE nyc_census_blocks_merge AS
SELECT ST_Union(geom) AS geom
FROM nyc_census_blocks;

SELECT ST_GeometryType(geom) FROM nyc_census_blocks_merge;

SELECT ST_NumGeometries(geom)
FROM nyc_census_blocks_merge;

📌练习：原点周围一单位缓冲区的面积是多少？它与您的预期有何不同？为什么？

SELECT ST_Area(ST_Buffer('POINT(0 0)', 1));

单位圆（半径为 1 的圆）的面积应为 pi，3.1415926… 差异是由于缓冲区边缘的线性描边。缓冲区具有有限数量的边。增加缓冲区中的边数将使值更接近 pi，但由于线性化，它总是会更小。

📌练习：‘Park Slope’ 和 ‘Carroll Gardens’ 的布鲁克林街区即将开战！构建一个多边形，在社区之间的边界上划定一个 100 米宽的 DMZ。非军事区的面积是多少？

CREATE TABLE brooklyn_dmz AS
SELECT
	ST_Intersection(ST_Buffer(ps.geom, 50), ST_Buffer(cg.geom, 50)) AS geom
FROM nyc_neighborhoods ps, nyc_neighborhoods cg
WHERE ps.name = 'Park Slope'
AND cg.name = 'Carroll Gardens';

SELECT ST_Area(geom) FROM brooklyn_dmz;

更多空间连接

创建人口普查表

在 Workshop 的\data\目录中，是一个包含属性数据但不包含几何图形的文件nyc_census_sociodata.sql，该表包括有关纽约的有趣社会经济数据（通勤时间、收入和教育程度）。

但是，数据需按人口普查区汇总，没有人口普查区空间数据！

创建人口普查表

可以通过汇总键的子字符串从人口普查块构建更高级别的几何blkid。为了获得人口普查区，需要对 . 的前 11 个字符进行汇总分组blkid。

360610001001001 = 36 061 000100 1 001

36     = State of New York
061    = New York County (Manhattan)
000100 = Census Tract
1      = Census Block Group
001    = Census Block

使用ST_Union聚合创建新表

-- 制作 tracts 表
CREATE TABLE nyc_census_tract_geoms AS
SELECT
  ST_Union(geom) AS geom,
  SubStr(blkid,1,11) AS tractid
FROM nyc_census_blocks
GROUP BY tractid;

-- 索引 tractid
CREATE INDEX nyc_census_tract_geoms_tractid_idx
ON nyc_census_tract_geoms (tractid);

将属性加入空间数据

使用常规属性将区域几何表连接到区域属性表

-- 制作 tracts 表
CREATE TABLE nyc_census_tracts AS
SELECT
  g.geom,
  a.*
FROM nyc_census_tract_geoms g
JOIN nyc_census_sociodata a
ON g.tractid = a.tractid;

-- 几何索引
CREATE INDEX nyc_census_tract_gidx
ON nyc_census_tracts USING GIST (geom);

计算纽约研究生学位比例排列前10的街区

SELECT
	n.name,
	n.boroname,
	100.0 * Sum(t.edu_graduate_dipl) / Sum(t.edu_total) AS graduate_pct
FROM nyc_neighborhoods n
JOIN nyc_census_tracts t
ON ST_Intersects(n.geom, t.geom)
WHERE t.edu_total > 0
GROUP BY n.name, n.boroname
ORDER BY graduate_pct DESC
LIMIT 10;

多边形/多边形连接

使用ST_Intersects(geometry_a,geometry_b)函数来确定在每个邻域中包含哪些人口普查区域多边形。

如果一块土地落在两个社区之间的边界上怎么办？它将与两者相交，因此将包含在两者的统计统计信息中。

通过以下两种方法避免重复计算

• 简单的方法：使用ST_Centroid(geometry)，确保每个区域只属于一个统计区域；
• 复杂的方法：使用ST_Intersection(geometry,geometry),在边界划分交叉区域。

📌示例：在研究生教育查询中使用简单方法避免重复计算的示例

SELECT
    n.name,
    n.boroname,
    100.0 * Sum(t.edu_graduate_dipl) / Sum(t.edu_total) AS graduate_pct
FROM nyc_neighborhoods n
JOIN nyc_census_tracts t
ON ST_Contains(n.geom, ST_Centroid(t.geom))
WHERE t.edu_total > 0
GROUP BY n.name, n.boroname
ORDER BY graduate_pct DESC
LIMIT 10;

Flatbush 社区的情况

Flatbush 社区已经从名单上掉了下来，可以通过更仔细地查看表中的 Flatbush 社区地图来了解。

正如数据源所定义的，Flatbush 社区并不是传统意义上的真正社区，因为它只是覆盖了 Prospect Park 的区域。该地区的人口普查区自然记录了零居民，然而，社区边界确实刮掉了与公园北侧接壤的昂贵人口普查区之一。当使用多边形/多边形测试时，这个单一区域被添加到原本空的 Flatbush 中，导致该查询的得分非常高。

大半径距离连接

一个有趣的问题是地铁站附近（500 米以内）人的通勤时间与远离地铁站的人的通勤时间有何不同？这个问题遇到了一些重复计算的问题：许多人将在多个地铁站的500米范围内。

纽约人口数量

SELECT Sum(popn_total)
FROM nyc_census_blocks;

📌示例：在地铁站 500 米范围内的纽约人口

SELECT Sum(popn_total)
FROM nyc_census_blocks census
JOIN nyc_subway_stations subway
ON ST_DWithin(census.geom, subway.geom, 500);

靠近地铁的人比总人口数量还多！

显然，我们的简单 SQL 犯了一个很大的重复计算错误，您可以查看缓冲地铁的图片看到问题。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NUOI7JKk-1646384052735)(https://postgis.net/workshops/postgis-intro/_images/subways_buffered.png)]

解决方案是确保在将它们传递到查询的摘要部分之前，只有不同的人口普查块。可以通过将查询分解成一个子查询来找到不同的块，并包装在一个返回我们答案的汇总查询中：

WITH distinct_blocks AS (
    SELECT DISTINCT ON (blkid) popn_total
    FROM nyc_census_blocks census
    JOIN nyc_subway_stations subway
    ON ST_DWithin(census.geom, subway.geom, 500)
)
SELECT Sum(popn_total)
FROM distinct_blocks;

纽约一半以上的人口位于地铁 500m 以内（步行约 5-7 分钟）。

有效性

定义

有效性定义了有界区域并且需要大量的结构。

符合 OGC SFSQL标准的有效性规则

• 任意类别
  • 环不能自相交（它们既不能接触也不能交叉）；
  • 戒指不得接触其他戒指，除非在某个点；
  • 多面体的元素不能相互接触；
• 多边形
  • 多边形环必须闭合；
  • 定义孔的环应位于定义外部边界的环内；

因为几何计算算法依赖于输入中一致结构，所以有效性规则很重要。

📌示例：此多边形无效，为什么结构很重要

POLYGON((0 0, 0 1, 2 1, 2 2, 1 2, 1 0, 0 0));

您可以在此图中更清楚地看到无效性：

图形程序成功地渲染了多边形填充，在视觉上它似乎是一个区域，外圈呈现一个8字形，中间有一个自相交点。

计算一下多边形的面积

SELECT ST_Area(ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))'));

计算面积的算法假定环不会自相交。一个表现良好的环将始终在边界线的一侧具有边界区域（内部）。然而在图形中，有界区域位于一条线的右侧，另一侧位于线的左侧。这会导致为每个瓣计算的面积抵消，因此得出零面积的结果。

检测

📌示例：使用ST_IsValid(geometry)函数检测表中数据的有效性

SELECT ST_IsValid(
	ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))')
);

📌示例：使用ST_IsValidReason(geometry)函数找出无效的原因

SELECT ST_IsValidReason(
	ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))')
);

📌示例：使用ST_IsValid(geometry)函数来测试表格中的数据

-- 查找所有无效多边形及其问题所在
SELECT
	name,
	boroname,
	ST_IsValidReason(geom) AS reason
FROM nyc_neighborhoods
WHERE NOT ST_IsValid(geom);

修复

修复无效性涉及将多边形分解为最简单的结构（环），确保环遵循有效性规则，然后构建遵循环封闭规则的新多边形。

结果通常是直观的，但在输入行为极其恶劣的情况下，有效输出可能不符合对它们外观的直觉。

📌示例：香蕉多边形是一个经典的无效性图形

POLYGON((0 0, 2 0, 1 1, 2 2, 3 1, 2 0, 4 0, 4 4, 0 4, 0 0))

在多边形上运行ST_MakeValid会返回一个有效的 OGC 多边形，由一个外环和内环在一个点接触。

SELECT ST_AsText(
    ST_MakeValid(
    	ST_GeometryFromText('POLYGON((0 0, 2 0, 1 1, 2 2, 3 1, 2 0, 4 0, 4 4, 0 4, 0 0))')
    )
);

香蕉多边形的有效性，在ESRI 模型认为接触的环是无效的；而 OGC 模型正好相反。
两者都不是 正确的，它们只是模拟相同情况的不同方法。

批量修复

📌示例：在向表中添加修复版本时标记无效几何以供查看

-- 旧无效表格的列
ALTER TABLE nyc_neighborhoods
ADD COLUMN geom_invalid geometry
DEFAULT NULL;

-- 修复无效并保存原件
UPDATE nyc_neighborhoods
SET geom = ST_MakeValid(geom), invalid_geom = geom
WHERE NOT ST_IsValid(geom);

-- 审核无效案例
SELECT geom, ST_IsValidReason(geom_invalid)
FROM nyc_neighborhoods
WHERE geom_invalid IS NOT NULL;

OpenJump是可视化修复无效几何的一个好工具，在Tools->QA->Validate Selected Layers下包含一个验证示例。

函数列表

• ST_IsValid(geometry A)：返回一个布尔值，指示几何图形是否有效。
• ST_IsValidReason(geometry A)：返回一个文本字符串，其中包含无效原因和无效坐标。
• ST_MakeValid(geometry A)：返回重新构造的几何图形以遵守有效性规则。

相等

相等

确定几何图形是否相等，是比较困难的。

为了说明这些功能，将使用以下多边形。

加载测试数据：

CREATE TABLE polygons (id integer, name varchar, poly geometry);

INSERT INTO polygons VALUES
(1, 'Polygon 1', 'POLYGON((-1 1.732,1 1.732,2 0,1 -1.732, -1 -1.732,-2 0,-1 1.732))'),
(2, 'Polygon 2', 'POLYGON((-1 1.732,-2 0,-1 -1.732,1 -1.732, 2 0,1 1.732,-1 1.732))'),
(3, 'Polygon 3', 'POLYGON((1 -1.732,2 0,1 1.732,-1 1.732, -2 0,-1 -1.732,1 -1.732))'),
(4, 'Polygon 4', 'POLYGON((-1 1.732,0 1.732, 1 1.732,1.5 0.866, 2 0,1.5 -0.866,1 -1.732,0 -1.732,-1 -1.732,-1.5 -0.866, -2 0,-1.5 0.866,-1 1.732))'),
(5, 'Polygon 5', 'POLYGON((-2 -1.732,2 -1.732,2 1.732, -2 1.732,-2 -1.732))');

完全相等

完全相等是通过逐个比较两个几何图形顶点的方式，确保它们的位置相同。

📌示例：查询了图形是否完全相等

SELECT 
	a.name,
	b.name,
	CASE WHEN ST_OrderingEquals(a.poly, b.poly) THEN 'Exactly Equal' ELSE 'Not Exactly Equal' END
FROM polygons AS a, polygons AS b;

在此示例中，多边形仅与自身相等，而不与其他看似等价的多边形（如多边形 1 到 3 的情况一样）。对于多边形 1、2 和 3，顶点的位置相同，但定义的顺序不同。多边形 4 在六边形边上具有共线（因此是冗余的）顶点，导致与多边形 1 不相等。

空间相等

📌示例：完全相等不考虑几何的空间性质，通过ST_Equals函数可用于检测几何空间相等

SELECT 
	a.name,
	b.name,
	CASE WHEN ST_Equals(a.poly, b.poly) THEN 'Spatially Equal' ELSE 'Not Equal' END
FROM polygons AS a, polygons AS b;

这些结果更符合我们对相等的直观理解。

注意：多边形的方向、定义多边形的起点、使用的点数在这里都不重要；重要的是多边形包含相同的空间。

等界

在最坏的情况下，精确相等需要比较几何中的每个顶点以确定相等。这可能很慢，并且可能不适合比较大量几何图形。为了加快比较速度，提供了~=等界运算符。这仅在边界框（矩形）上运行，确保几何图形占据相同的二维范围，但不一定是相同的空间。

📌示例：所有空间相等的几何图形也具有相同的界限。

SELECT 
	a.name,
	b.name,
	CASE WHEN a.poly ~= b.poly THEN 'Equal Bounds' ELSE 'Non-equal Bounds' END
FROM polygons AS a, polygons AS b;

但是 Polygon 5 与其他几何图形共享相同的边界框，测试中也返回为相等。

适用场景为：~=使用空间索引，在连接或过滤数据时，可快速比较大量数据。

线性参考

线性参考（动态分割）是一种表示特征的方法，可以通过引用一组基本的线性特征来描述。

以下示例使用线性参考建模的特征

• 高速公路资产，使用沿高速公路网络的英里数进行引用；
• 道路维护操作，被称为在一对英里测量之间沿着道路网络发生；
• 水生清单，其中鱼类的存在被记录为在一对上游里程测量之间存在；
• 溪流的水文特征，用起点和终点记录。

优点：依赖空间观测，不需要与基础观测分开记录，可以在知道依赖观测将自动跟踪新几何的情况下对基础观测层进行更新。

Esri 线性参考：
具有线性空间要素的基表和事件的非空间表，其中包括对空间要素的外键参考和沿参考要素的度量;
使用术语事件表来指代构建的非空间表。

创建线性参考

使用ST_LineLocatePoint函数引用线性网络的表，该函数接受点、线，并返回沿线可以找到点的比例。

📌示例：沿线定位一个点的简单示例

SELECT ST_LineLocatePoint('LINESTRING(0 0, 2 2)', 'POINT(1 1)');

-- 如果点不在线上，把它投影到最近点
SELECT ST_LineLocatePoint('LINESTRING(0 0, 2 2)', 'POINT(0 2)');

📌示例：使用ST_LineLocatePoint将 nyc_subway_stations 转换为相对于街道的事件表 。

-- 下面的所有 SQL 都有助于创建新的事件表
CREATE TABLE nyc_subway_station_events AS
-- 首先需要得到一组候选的可能最近的街道，按 id 和距离排序......
WITH ordered_nearest AS (
SELECT
  ST_GeometryN(streets.geom,1) AS streets_geom,
  streets.gid AS streets_gid,
  subways.geom AS subways_geom,
  subways.gid AS subways_gid,
  ST_Distance(streets.geom, subways.geom) AS distance
FROM nyc_streets streets
  JOIN nyc_subway_stations subways
  ON ST_DWithin(streets.geom, subways.geom, 200)
ORDER BY subways_gid, distance ASC
)
-- 使用 'distinct on' 功能来获取每个唯一街道 gid 的第一条街道（最近的）。 
-- 然后将那条街道连同其候选地铁站一起传递到 ST_LineLocatePoint 以计算度量。
SELECT
  DISTINCT ON (subways_gid)
  subways_gid,
  streets_gid,
  ST_LineLocatePoint(streets_geom, subways_geom) AS measure,
  distance
FROM ordered_nearest;

-- 主键对可视化软件很有用
ALTER TABLE nyc_subway_station_events ADD PRIMARY KEY (subways_gid);

一旦有了一个事件表，将它转回空间视图是很有趣的，就可以可视化相对于它们派生的原始点的事件。

要从一个度量到一个点，使用ST_LineInterpolatePoint函数。

📌示例：这是之前颠倒的简单示例

SELECT ST_AsText(ST_LineInterpolatePoint('LINESTRING(0 0, 2 2)', 0.5));

可以将nyc_subway_station_events表连接回nyc_streets表并使用measure属性生成空间事件点，而无需参考原始nyc_subway_stations表。

-- 将事件变回空间对象的新视图
CREATE OR REPLACE VIEW nyc_subway_stations_lrs AS
SELECT
  events.subways_gid,
  ST_LineInterpolatePoint(ST_GeometryN(streets.geom, 1), events.measure)AS geom,
  events.streets_gid
FROM nyc_subway_station_events events
JOIN nyc_streets streets
ON (streets.gid = events.streets_gid);

查看街道的原始（红色星形）和事件（蓝色圆圈）点，可以看到事件如何直接捕捉到最近的街道线。

线性参考函数的一个令人惊讶的用途与线性参考模型无关。
如上所示，可以使用这些函数将点捕捉到线性要素。
对于 GPS 轨迹或其他预期参考线性网络的输入等用例，捕捉是一个方便的功能。

函数列表

• ST_LineInterpolatePoint(geometry A, double measure)：返回沿直线插值的点。
• ST_LineLocatePoint(geometry A, geometry B)：返回一个介于 0 和 1 之间的浮点数，表示 LineString 上离给定点最近的点的位置。
• ST_LineSubstring(geometry A, double from, double to)：返回一个线串，它是输入的子串，以二维总长度的给定分数开始和结束。
• ST_LocateAlong(geometry A, double measure)：返回派生几何集合值，其中包含与指定度量匹配的元素。
• ST_LocateBetween(geometry A, double from, double to)：返回一个派生的几何集合值，其中包含与指定的度量范围匹配的元素。
• ST_AddMeasure(geometry A, double from, double to)：返回一个派生几何，测量元素在起点和终点之间线性插值。如果几何没有测量维度，则添加一个。

9交叉模型

DE9IM(Dimensionally Extended 9-Intersection Model) 是一个用于对空间对象交互进行建模的框架。

多边形都由内部（由环包围）、边界（环本身）、外部（环外的部分）组成

其空间特征由下图所示，内部（由端点界定的线的一部分）；边界（线性特征的末端），外部（平面中的其他所有内容）。

对于点: 内部就是点，边界是空集，外部是平面中的所有其他内容。

对于任意两个空间特征之间的关系，都可以使用内部、边界、外部之间的九个可能交点的维度来表示。

对于上面示例中的多边形，内部的交集是一个二维区域，因此矩阵的该部分用2填充；边界仅在零维的点处相交，因此矩阵的该部分用0填充。

当组件之间没有交集时，矩阵的正方形用F填充。

这是另一个线串部分进入多边形的示例

交互的 DE9IM 矩阵是这样的

注意：两个对象的边界实际上根本不相交（线的端点与多边形内部交互，而不是边界，反之亦然），因此 B/B 单元格填充了F。

前面的示例可以使用简单的框和线进行简化，与多边形和线串具有相同的空间关系：

ST_Relate函数可以生成 DE9IM 矩阵信息：

SELECT ST_Relate('LINESTRING(0 0, 2 0)','POLYGON((1 -1, 1 1, 3 1, 3 -1, 1 -1))');

答案为（1010F0212）和我们肉眼计算的一样，但是返回的是一个 9 个字符的字符串，表格的第一行、第二行和第三行附加在一起。

然而，DE9IM 矩阵的强大之处不在于生成它们，而在于将它们用作匹配键来查找彼此之间具有特定关系的几何图形。

CREATE TABLE lakes ( id serial primary key, geom geometry );
CREATE TABLE docks ( id serial primary key, good boolean, geom geometry );

INSERT INTO lakes ( geom ) VALUES ( 'POLYGON ((100 200, 140 230, 180 310, 280 310, 390 270, 400 210, 320 140, 215 141, 150 170, 100 200))');
INSERT INTO docks ( geom, good ) VALUES
('LINESTRING (170 290, 205 272)',true),
('LINESTRING (120 215, 176 197)',true),
('LINESTRING (290 260, 340 250)',false),
('LINESTRING (350 300, 400 320)',false),
('LINESTRING (370 230, 420 240)',false),
('LINESTRING (370 180, 390 160)',false);

📌示例：假设有一个包含湖、码头的数据模型，并进一步假设码头必须在湖内，且必须在一端触及包含湖的边界。

能在数据库中找到所有遵守该规则的码头吗？

符合条件的码头具有以下特点

• 码头的内部与湖内部有一个线性 (1D) 交叉点；
• 码头的边界与湖内部有一个点 (0D) 交叉点；
• 码头的边界也与湖泊边界有一个点 (0D) 相交；
• 码头的内部与湖的外部没有交集（F）。

符合条件的码头的 DE9IM 矩阵是这样的

要找到所有的合法码头，我们需要找到与湖相交的所有码头（我们用于连接键的潜在候选者的超集），然后找到该集合中具有法律相关模式的所有码头。

SELECT
	docks.*
FROM docks 
JOIN lakes 
ON ST_Intersects(docks.geom, lakes.geom)
WHERE ST_Relate(docks.geom, lakes.geom, '1FF00F212');

📌示例：如果没有在示例图形中包含的一个可能的码头是与湖泊边界具有二维交叉点的码头

INSERT INTO docks ( geom, good ) VALUES ('LINESTRING (140 230, 150 250, 210 230)',true);

如果我们要将这个案例包含在我们的合法停靠集中，我们需要更改查询中的关联模式；码头内湖边界的交点现在可以是 1（我们的新案例）或 F（我们的原始案例）。所以我们在模式中使用了*。

SELECT
	docks.*
FROM docks 
JOIN lakes 
ON ST_Intersects(docks.geom, lakes.geom)
WHERE ST_Relate(docks.geom, lakes.geom, '1*F00F212');

确认上一个示例中更严格的 SQL不会返回新的停靠栏。

数据质量验证

例如：人口普查区块不应与任何其他人口普查区块重叠

SELECT 
	a.gid,
	b.gid
FROM nyc_census_blocks a, nyc_census_blocks b
WHERE ST_Intersects(a.geom, b.geom)
AND ST_Relate(a.geom, b.geom, '2********')
AND a.gid != b.gid
LIMIT 10;

同样，我们希望道路数据都是端节点的，期望相交只发生在线的末端，而不是中点。

可以通过寻找相交的街道（因此我们有一个连接）但边界之间的交点不是零维（即端点不接触）来验证这一点：

SELECT 
	a.gid,
	b.gid
FROM nyc_streets a, nyc_streets b
WHERE ST_Intersects(a.geom, b.geom)
AND NOT ST_Relate(a.geom, b.geom, '****0****')
AND a.gid != b.gid
LIMIT 10;

函数列表

ST_Relate(geometry A, geometry B)：返回一个文本字符串，表示几何之间的 DE9IM 关系。

指数聚类

数据库需要尽快从磁盘中获取信息，小型数据库将完全浮动到 RAM 缓存中，并摆脱物理磁盘的限制；但对于大型数据库，访问物理磁盘将限制磁盘访问速度。

数据是随机写入磁盘的，因此存储在磁盘上的订单数据与应用程序访问或组织数据的方式之间不一定有任何关联。

聚类：确保在同一结果集中一起检索的记录位于硬盘盘片上相似的物理位置。

通用规则：索引为数据定义了一种自然排序方案，类似于检索数据时使用的访问模式。

在某些情况下，以与索引相同的顺序对磁盘上的数据进行排序可以提供速度优势。

R-Tree 上的聚类

基于空间索引的聚类对于将要通过空间查询访问的空间数据是有意义的：相似的事物往往具有相似的位置。

📌示例：根据nyc_census_blocks的空间索引对进行聚类

-- 根据空间索引对块进行聚类
CLUSTER nyc_census_blocks USING nyc_census_blocks_geom_gist;

nyc_census_blocks该命令按照空间索引定义的顺序重新写入nyc_census_blocks_geom_gist。

在空间数据上增量构建的 R-Tree 可能不具有叶子的高空间一致性。例如，请参阅不列颠哥伦比亚省道路索引的空间索引叶的可视化。

平衡 R-Tree 是空间更紧凑的树进行聚类

在 PostGIS 中没有可用的平衡 R-Tree 算法，但可以通过ST_GeoHash()函数使空间数据置于空间自相关顺序排列。

GeoHash 上的聚类

要在ST_GeoHash()函数上进行聚类，首先需要对数据进行 geohash 索引；geohash 算法仅适用于地理（经度/纬度）坐标中的数据，需要在散列的同时转换几何（到 EPSG:4326，即经度/纬度）。

CREATE INDEX nyc_census_blocks_geohash 
ON nyc_census_blocks(ST_GeoHash(ST_Transform(geom, 4326)));

一旦有了一个 geohash 索引，在它上面的聚类使用与 R-Tree 聚类相同的语法。

CLUSTER nyc_census_blocks USING nyc_census_blocks_geohash;

现在的数据以空间相关的顺序排列！

函数列表

ST_GeoHash(geometry A)：返回一个文本字符串，表示对象边界的 GeoHash。

三维

三维几何

PostGIS 支持所有几何类型的附加维度，Z维度用于添加高度信息，M维度用于每个坐标的附加维度信息（通常是时间、或道路英里或上游距离信息）。

对于三维/四维几何，额外的维度作为几何中每个顶点的额外坐标添加，并且几何类型被增强。

添加额外的维度会为每个几何图元产生三种额外的可能几何类型：

• Point（二维类型）由 PointZ、PointM 和 PointZM 类型连接；
• Linestring（一种二维类型）由 LinestringZ、LinestringM 和 LinestringZM 类型连接；
• Polygon（二维类型）由 PolygonZ、PolygonM 和 PolygonZM 类型连接；
• 等等。

WKT中更高维几何的格式由 ISO SQL/MM 规范给出，额外的维度信息只是简单地添加到类型名称之后的文本字符串中，而额外的坐标则添加在 X/Y 信息之后。

例如：

• POINT ZM (1 2 3 4)
• LINESTRING M (1 1 0, 1 2 0, 1 3 1, 2 2 0)
• POLYGON Z ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0))

ST_AsText() 函数将在处理 3-D 和 4-D 几何图形时返回上述表示。

WKB中更高维几何的格式由 ISO SQL/MM 规范给出，该格式的BNF。

除了标准类型的高维形式之外，PostGIS 还包括一些在三维空间中有意义的新类型：

• TIN 类型：允许您将三角形网格建模为数据库中的行；
• POLYHEDRALSURFACE：允许您对数据库中的体积对象进行建模。

由于这两种类型都用于建模 3-D 对象，因此使用 Z 变体才真正有意义。

📌示例：1 个单位的 POLYHEDRALSURFACE Z 立方体：

POLYHEDRALSURFACE Z (
  ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
  ((0 0 0, 0 1 0, 0 1 1, 0 0 1, 0 0 0)),
  ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
  ((1 1 1, 1 0 1, 0 0 1, 0 1 1, 1 1 1)),
  ((1 1 1, 1 0 1, 1 0 0, 1 1 0, 1 1 1)),
  ((1 1 1, 1 1 0, 0 1 0, 0 1 1, 1 1 1))
)

三维函数

有许多函数用于计算 3-D 对象之间的关系：

• ST_3DClosestPoint — 返回 g1 上最接近 g2 的 3 维点。这是 3D 最短线的第一个点。
• ST_3DDistance — 对于几何类型返回投影单位中两个几何之间的 3 维笛卡尔最小距离（基于空间参考）。
• ST_3DDWithin — 对于 3d (z) 几何类型如果两个几何 3d 距离在单位数内，则返回 true。
• ST_3DDFullyWithin — 如果所有 3D 几何都在彼此之间的指定距离内，则返回 true。
• ST_3DIntersects — 如果几何在 3d 中“空间相交”，则返回 TRUE - 仅适用于点和线串
• ST_3DLogestLine — 返回两个几何图形之间的 3 维最长线
• ST_3DMaxDistance — 对于几何类型返回投影单位中两个几何之间的 3 维笛卡尔最大距离（基于空间参考）。
• ST_3DShortestLine — 返回两个几何之间的 3 维最短线

📌示例：使用 ST_3DDistance 函数计算单位立方体与点之间的距离：

-- 这实际上是顶角和点之间的距离
SELECT ST_3DDistance(
  'POLYHEDRALSURFACE Z (
    ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
    ((0 0 0, 0 1 0, 0 1 1, 0 0 1, 0 0 0)),
    ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
    ((1 1 1, 1 0 1, 0 0 1, 0 1 1, 1 1 1)),
    ((1 1 1, 1 0 1, 1 0 0, 1 1 0, 1 1 1)),
    ((1 1 1, 1 1 0, 0 1 0, 0 1 1, 1 1 1))
  )'::geometry,
  'POINT Z (2 2 2)'::geometry
);

-- 所以这里有一个较短的形式
SELECT ST_3DDistance(
  'POINT Z (1 1 1)'::geometry,
  'POINT Z (2 2 2)'::geometry
);

-- 两者都按预期返回 1.73205080756888 == sqrt(3)

ND 索引

一旦有更高维度的数据，索引它可能是有意义的。在应用多维索引之前，应该仔细考虑数据在所有维度中的分布。

索引仅在通过 WHERE 条件大幅减少返回行数时才有用，为了使更高维度的索引有用，数据必须涵盖该维度的广泛范围。

• DEM 点可能不适合三维索引，因为查询通常会提取二维点框，并且很少尝试选择 Z 切片点；
• 如果 GPS 轨迹在所有维度上经常相互重叠（例如，在不同时间反复行驶相同的路线），则X/Y/T 空间中的 GPS 轨迹可能是三维索引的良好候选者，因为数据集的所有维度都会有很大的可变性。

可以对任何维度（甚至是混合维度）的数据创建多维索引。

📌示例：在nyc_streets表上创建多维索引

CREATE INDEX nyc_streets_gix_nd ON nyc_streets
USING GIST(geom gist_geometry_ops_nd);

该gist_geometry_ops_nd参数告诉 PostGIS 使用 ND 索引而不是标准的二维索引。

建立索引后，查询中使用&&&索引运算符，&&&具有与边界框交互相同的语义&&，但使用输入几何的所有维度应用这些语义，具有不匹配维度的几何图形不会相互作用。

-- 返回 true（在零平面上都是 3D）
SELECT 'POINT Z (1 1 0)'::geometry &&& 'POLYGON ((0 0 0, 0 2 0, 2 2 0, 2 0 0, 0 0 0))'::geometry;

-- 返回 false（一个 2D 一个 3D）
SELECT 'POINT Z (3 3 3)'::geometry &&& 'POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))'::geometry;

-- 返回 true（线串周围的体积与点交互）
SELECT 'LINESTRING Z(0 0 0, 1 1 1)'::geometry &&& 'POINT(0 1 1)'::geometry;

要nyc_streets使用 ND 索引搜索表，只需将常用的&&二维索引运算符替换为&&&运算符即可。

-- ND索引
SELECT gid, name
FROM nyc_streets
WHERE geom &&& ST_SetSRID('LINESTRING(586785 4492901,587561 4493037)'::geometry,26918);

-- 二维索引
SELECT gid, name
FROM nyc_streets
WHERE geom && ST_SetSRID('LINESTRING(586785 4492901,587561 4493037)'::geometry,26918);

结果是一样的。

通常 ND索引 比 二维索引 慢得多，因此仅在确定 ND 查询会提高查询的选择性的情况下使用 ND 索引。

邻域搜索

邻域搜索

一个经常提出的空间查询是：离 [查询特征] 最近的 [候选特征] 是什么？

与距离搜索不同，领域搜索不包括任何限制候选几何可能距离的测量，任何距离的特征都将被接受，只要它们是最近的。

通过引入按距离排序 (<->) 运算符来解决邻域问题，使得数据库使用索引来加速排序的返回集。使用按距离排序运算符，邻域查询只需添加排序并将结果集限制为 N 个条目即可返回 N 个最近的特征。

按距离排序运算符适用于几何和地理坐标类型。

• 几何：返回结果取决于所使用的空间参考系统的单位；
• 地理：返回结果是球体距离。

📌示例：查询距离 Broad St 地铁站最近的 3 条街道

-- 获取 Broad St 的几何图形
SELECT 
	ST_AsEWKT(geom, 1)
FROM nyc_subway_stations
WHERE name = 'Broad St';

-- 将几何图形插入邻域查询
SELECT
	streets.gid, streets.name,
	ST_Transform(streets.geom, 4326),
	streets.geom <-> 'SRID=26918;POINT(583571.9 4506714.3)'::geometry AS dist
FROM nyc_streets streets
ORDER BY dist
LIMIT 3;

检查邻域查询的输出是一个好主意EXPLAIN，可以从非索引 SQL 中获得正确的答案，并且在表的大小扩大之前，缺少索引可能并不明显。

📌示例：使用 EXPLAIN 输出索引扫描

EXPLAIN SELECT
	streets.gid, streets.name,
	ST_Transform(streets.geom, 4326),
	streets.geom <-> 'SRID=26918;POINT(583571.9 4506714.3)'::geometry AS dist
FROM nyc_streets streets
ORDER BY dist
LIMIT 3;

邻域连接

运算符的索引辅助排序有一个主要缺点：仅适用于运算符一侧的单个几何字符；对于查找离一个查询对象最近的对象来说很好，但对于空间连接没有帮助，空间连接的目标是为一组完整的候选对象中的每一个找到最近的邻居。

LATERAL join允许在循环中重复运行一个查询。

📌示例：查询离每个地铁站最近的街道

SELECT subways.gid AS subway_gid,
       subways.name AS subway,
       streets.name AS street,
       streets.gid AS street_gid,
       streets.geom::geometry(MultiLinestring, 26918) AS street_geom,
       streets.dist
FROM nyc_subway_stations subways
CROSS JOIN LATERAL (
  SELECT 
	streets.name, 
	streets.geom, 
	streets.gid, 
	streets.geom <-> subways.geom AS dist
  FROM nyc_streets AS streets
  ORDER BY dist
  LIMIT 1
) streets;

注意：CROSS JOIN LATERAL充当由地铁表驱动的循环的内部部分的方式，地铁表中的每条记录都被输入横向子查询，一次一条，因此您可以获得每条地铁记录的最接近的结果。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PJDtmgE8-1646623661713)(https://postgis.net/workshops/postgis-intro/_images/knn4.png)]

使用 EXPLAIN 输出索引扫描

日志审计（触发器实现）

日志审计是生产数据库一个常见要求，数据在两个日期之间发生了怎样的变化，谁进行了更改，以及它们发生在哪里？一些 GIS 系统通过在客户端界面中包含变更管理来跟踪变更，但这增加了编辑工具的复杂性。

使用数据库和触发器，可将历史跟踪添加到任何表，同时保持对主表的简单直接编辑访问。

历史跟踪的工作原理是为每次编辑保留一个记录历史记录表

• 创建记录，添加时间和创建者；
• 删除记录，何时删除以及由谁删除；
• 更新记录，则添加删除记录（用于旧状态）和创建记录（用于新状态）。

TSTZRANGE

历史表使用时间戳范围类型来实时记录、存储历史记录，历史表中时间戳范围可预期、相邻、不重叠的。

新记录的范围将开始于now()并有一个开放的结束点，该范围涵盖了从当前时间到未来的所有时间。

SELECT tstzrange(current_timestamp, NULL);

               tstzrange
------------------------------------
 ["2021-06-01 14:49:40.910074-07",)

搜索时间范围比搜索一对时间戳要简单得多，时间范围从起点到无穷大的所有时间，包含范围使用@>运算符。

-- 十分钟前的范围包括现在吗？
-- 这应该！ :)
SELECT tstzrange(current_timestamp - '10m'::interval, NULL) @> current_timestamp;

构建历史表

使用此信息，可以在任何时间点重建编辑表的状态。在此示例中，我们将向nyc_streets表添加历史跟踪。

1️⃣ 添加一个新的nyc_streets_history表，将用来存储所有历史编辑信息的表；

DROP TABLE IF EXISTS nyc_streets_history;
CREATE TABLE nyc_streets_history (
    -- 主键
    hid SERIAL PRIMARY KEY,
    gid INTEGER,
    id FLOAT8,
    name VARCHAR(200),
    oneway VARCHAR(10),
    type VARCHAR(50),
    geom GEOMETRY(MultiLinestring,26918),
    -- 记录活动的时间范围
    valid_range TSTZRANGE,
    -- 创建用户
    created_by VARCHAR(32),
    -- 删除用户
    deleted_by VARCHAR(32)
);

CREATE INDEX nyc_streets_history_geom_x
ON nyc_streets_history USING GIST (geom);

CREATE INDEX nyc_streets_history_tstz_x
ON nyc_streets_history USING GIST (valid_range);

2️⃣ 将活动表 nyc_streets 的当前状态导入到历史表中，这样就有了跟踪历史的起点。

注意：填写创建时间和创建用户，但将时间范围的结束和删除的信息保留为NULL。

INSERT INTO nyc_streets_history(gid, id, name, oneway, type, geom, valid_range, created_by)
SELECT gid, id, name, oneway, type, geom, tstzrange(now(), NULL),current_user
FROM nyc_streets;

3️⃣ 现在需要在活动表上使用三个触发器，用于 INSERT、DELETE 和 UPDATE 操作。

首先我们创建触发器函数，然后将它们作为触发器绑定到表。

对于插入，只需在历史表中添加一条带有创建时间/用户的新记录。

CREATE OR REPLACE FUNCTION nyc_streets_insert() RETURNS trigger AS
$$
    BEGIN
    INSERT INTO nyc_streets_history(gid, id, name, oneway, type, geom, valid_range, created_by) VALUES(NEW.gid, NEW.id, NEW.name, NEW.oneway, NEW.type, NEW.geom,tstzrange(current_timestamp, NULL), current_user);
    RETURN NEW;
    END;
$$
LANGUAGE plpgsql;

CREATE TRIGGER nyc_streets_insert_trigger
AFTER INSERT ON nyc_streets
FOR EACH ROW EXECUTE PROCEDURE nyc_streets_insert();

对于删除，只需将当前活动的历史记录（删除时间为 NULL 的记录）标记为已删除。

CREATE OR REPLACE FUNCTION nyc_streets_delete() RETURNS trigger AS
$$
    BEGIN
    UPDATE nyc_streets_history SET valid_range = tstzrange(lower(valid_range), current_timestamp), deleted_by = current_user
    WHERE valid_range @> current_timestamp AND gid = OLD.gid;
    RETURN NULL;
    END;
$$
LANGUAGE plpgsql;

CREATE TRIGGER nyc_streets_delete_trigger
AFTER DELETE ON nyc_streets
FOR EACH ROW EXECUTE PROCEDURE nyc_streets_delete();

对于更新，首先将活动历史记录标记为已删除，然后为更新状态插入一条新记录。

CREATE OR REPLACE FUNCTION nyc_streets_update() RETURNS trigger AS
$$
    BEGIN

    UPDATE nyc_streets_history SET valid_range = tstzrange(lower(valid_range), current_timestamp), deleted_by = current_user
    WHERE valid_range @> current_timestamp AND gid = OLD.gid;

    INSERT INTO nyc_streets_history(gid, id, name, oneway, type, geom, valid_range, created_by) VALUES(NEW.gid, NEW.id, NEW.name, NEW.oneway, NEW.type, NEW.geom, tstzrange(current_timestamp, NULL), current_user);

    RETURN NEW;

    END;
$$
LANGUAGE plpgsql;

CREATE TRIGGER nyc_streets_update_trigger
AFTER UPDATE ON nyc_streets
FOR EACH ROW EXECUTE PROCEDURE nyc_streets_update();

编辑表格

启用历史表，可以在主表上进行编辑并观察日志条目出现在历史表中。

注意：无论使用什么工具进行编辑，无论是 SQL 命令行、基于 Web 的 JDBC 工具，还是 QGIS 等桌面工具，历史记录都会得到一致的跟踪。

SQL 编辑

📌示例：把两条名为Cumberland Walk的街道改成更时尚的Cumberland Wynde

更新这两条街道会导致在历史表中将原来的街道标记为已删除，删除时间为现在，并添加两条新街道名称，添加时间为现在。

查询历史表

查询特定时间T的历史表，结果包括：

• 在 T 之前创建但尚未删除的所有记录；
• 在 T 之前创建但在T之后删除的所有记录。

📌示例：创建10分钟前历史表视图

-- 有效范围包括 10 分钟前的记录是当时有效的记录
CREATE OR REPLACE VIEW nyc_streets_ten_min_ago AS
SELECT * FROM nyc_streets_history 
WHERE valid_range @> (now() - '10min'::interval)

📌示例：创建特定用户添加的历史表视图

CREATE OR REPLACE VIEW nyc_streets_postgres AS
SELECT * FROM nyc_streets_history
WHERE created_by = 'postgres';

另见

• QGIS开源GIS
• PostgreSQL 触发器
• PostgreSQL 范围类型

基本 PostgreSQL 调优

为了确保 PostgreSQL 能够在许多不同的环境中正常运行，默认配置非常保守，不太适合高性能生产数据库；再加上地理空间数据库具有不同的使用模式，可以看到默认配置并不完全适合。

所有这些配置参数都可以在 postgresql.conf 数据库配置文件中进行编辑，修改后重启数据库生效。

共享缓冲区 shared_buffers

用于共享内存缓冲区的内存量，默认值（32MB）对于生产环境不合适，推荐配置为：约 75% 的数据库内存，最大约 2GB。

有效缓存大小 effective_cache_size

PostgreSQL 考虑操作系统可能缓存了多少磁盘块作为其虚拟文件系统的一部分。

对于具有大量内存的系统，这可能非常大。大约是机器上的effective_cache_size内存量，减去shared_buffers，减去work_mem预期连接数的倍数，减去机器上运行的任何其他进程所需的任何内存，其他随机操作系统需要大约 1GB。数据库不会直接使用额外的缓存，但它会计算预期操作系统已将文件系统数据缓存在大约那么多内存中的计划。

默认值：通常为 4GB; 推荐值：在普通操作条件下预期存在的任意数量的“空闲”内存

工作记忆 work_mem

在数据保存到磁盘文件之前，内部排序操作、索引操作、哈希表会消耗的内存量；该值定义了每个操作的可用内存；复杂的查询可能有多个并行运行的排序或散列操作，并且每个连接的会话可能正在执行一个查询。

必须在增加此值之前考虑有多少连接和预期查询的复杂性，增大后可以在不产生磁盘写入的情况下完成更多这些操作的处理，包括 ORDER BY 和 DISTINCT 子句、合并和哈希连接、基于哈希的聚合和基于哈希的子查询处理。

默认值：1MB; 推荐值：32MB

工作记忆维护 maintenance_work_mem

用于维护操作的内存量，包括清理、索引和外键创建，由于这些操作不是很常见，因此较高的值只会产生偶然的成本，并且可能会大大加快维护活动的速度，在执行多个CREATE INDEX 或 VACUUM 调用之前，可以在单个会话中交替增加此参数。

SET maintenance_work_mem TO '128MB';
VACUUM ANALYZE;
SET maintenance_work_mem TO '16MB';

默认值：16MB；推荐值：128MB

预写日志缓冲区 wal_buffers

设置用于预写日志 (WAL) 数据的内存量，预写日志提供了一种确保数据完整性的高性能机制。在每个更改命令期间，更改的效果首先写入 WAL 文件并刷新到磁盘。只有在刷新 WAL 文件后，才会将更改写入数据文件本身。这允许以最佳、异步的方式将数据文件写入磁盘，同时确保在发生崩溃时，可以从 WAL 中恢复所有数据更改。

这个缓冲区的大小只需要大到足以容纳单个典型事务的 WAL 数据。虽然默认值通常足以满足大多数数据，但地理空间数据往往要大得多。

默认值：64kB；推荐值：1MB

检查点部分 checkpoint_segments

该值设置在自动 WAL 检查点之间可以填充的最大日志文件段数（通常为 16MB）。WAL 检查点是 WAL 事务序列中的一个点，在该点可以保证数据文件已使用检查点之前的所有信息进行了更新。此时所有脏数据页都被刷新到磁盘，并且检查点记录被写入日志文件。这允许崩溃恢复过程找到最新的检查点记录并应用所有后续日志段来完成数据恢复。

因为检查点进程需要将所有脏数据页刷新到磁盘，所以会产生很大的 I/O 负载。上面的相同论点也适用；地理空间数据大到足以使非地理空间优化失衡。增加此值将防止过多的检查点，但它可能会导致服务器在崩溃时重新启动更慢。

默认值：3；推荐值：6

随机页成本 random_page_cost

这是一个无单位值，表示从磁盘访问随机页面的成本。该值与许多其他成本参数相关，包括顺序页面访问和 CPU 操作成本。虽然这个值没有灵丹妙药，但默认值通常是保守的，并且适用于在旋转介质上运行的数据库。SSD 的随机访问成本应该设置得更低。

可以使用SET random_page_cost TO 2.0命令在每个会话的基础上设置此值，这对于测试它如何影响查询计划很有用。

默认值：4.0；推荐值：2.0 用于旋转介质，1.0 用于 SSD

序列页成本 seq_page_cost

这是控制顺序页面访问成本的参数。这个值一般不需要调整，但是这个值之间的差异random_page_cost会极大地影响查询计划器的选择。该值也可以在每个会话的基础上设置。

默认值：1.0；推荐值：1.0

重新加载配置

进行这些更改后，保存更改并重新加载配置；最简单的方法是重新启动 PostgreSQL 服务。

• 在 Windows 服务 ( services.msc) 中，右键单击PostgreSQL并选择重新启动。

模式机制

生产数据库不可避免地有大量的表和视图，并且在一个模式中管理它们会很快变得笨拙。

PostgreSQL使用 Schema 模式机制，就像文件夹，可以保存表、视图、函数、序列和其他关系。每个数据库都从一个模式开始，即public模式。

在该模式中，PostGIS 的默认安装创建geometry_columns和元数据关系，geography_columns以及spatial_ref_sys 使用的所有类型和函数。

所以 PostGIS 的用户总是需要访问公共模式。在公共模式中，还可以看到迄今为止创建的所有表。

模式的优势

使用模式有优势：

1. 在架构中管理的数据更容易应用批量操作
   a. 备份单独模式中的数据更容易：因此易失性数据可以有与非易失性数据不同的备份计划。
   b. 恢复单独模式中的数据更容易：因此可以单独恢复和备份面向应用程序的模式，以进行时间旅行和恢复。
   c. 当应用程序数据位于模式中时，管理应用程序差异会更容易：因此新版本的软件可以处理新模式中的表结构，并且切换涉及对模式名称的简单更改。
2. 可以将用户的工作限制为单一模式，以允许将分析表和测试表与生产表隔离。

因此，出于生产目的，将应用程序数据在模式中分开可以改善管理；出于用户目的，将用户保持在单独的模式中可以防止他们相互干扰。

创建数据模式

让我们创建一个新模式并将一个表移动到其中。首先，在数据库中创建一个新模式：

CREATE SCHEMA census;

接下来，把nyc_census_blocks表移动到census模式中：

ALTER TABLE nyc_census_blocks SET SCHEMA census;

可以通过两种方式访问模式中的表：

• schema.table通过使用符号引用；
• 通过将架构添加到您的search_path

显式引用很容易，但一段时间后输入会很累：

SELECT * FROM census.nyc_census_blocks LIMIT 1;

操作search_path是提供对多个模式中的表的访问而无需大量额外输入的好方法。

可以使用以下命令设置search_path运行时：SET

SET search_path = census, public;

这可确保在census和public模式中搜索对关系和函数的所有引用。

每次连接时设置搜索路径也会很累，可以设置永久搜索路径：

ALTER USER postgres SET search_path = census, public;

现在 postgres 用户将始终census在其搜索路径中包含架构。

创建用户模式

使用 SQL 的分析操作需要临时表来实现可视化或临时结果，因此空间 SQL 往往要求用户拥有比普通数据库工作负载更多的 CREATE 权限。

默认情况下，Oracle 中的每个角色都有一个个人模式。这对于 PostgreSQL 用户来说也是一个很好的实践，并且很容易使用 PostgreSQL 角色、模式和搜索路径进行复制。

创建具有表创建权限的新用户，然后使用该用户作为授权创建模式：

CREATE USER myuser WITH ROLE postgis_writer;
CREATE SCHEMA myuser AUTHORIZATION myuser;

如果以该用户身份登录，会发现search_path默认设置实际上是这样的：

show search_path;

  search_path
----------------
 "$user",public

搜索路径上的第一个模式是用户的命名模式！所以现在存在以下条件：

• 用户存在，具有创建空间表的能力。
• 用户的命名模式存在，并且用户拥有它。
• 用户的搜索路径首先具有用户模式，因此会自动在那里创建新表，并且查询会首先自动搜索那里。

备份与恢复

备份数据库的方法有很多种，而选择的方法很大程度上取决于使用数据库的方式。

• 对于相对静态的数据库，可以使用基本的 pg_dump/pg_restore 工具对数据进行定期快照；
• 对于频繁更改的数据，使用热备方案可以将更新连续存档到安全位置。

热备是高可用性复制和备用系统的基础，特别是对于 PostgreSQL >= 9.0 的版本。

布局数据

确保生产数据始终存储在单独的模式中是管理数据的一个非常重要的最佳实践，有两个原因：

• 在模式中备份和恢复数据比管理要单独备份的表列表要简单得多；
• 将数据表排除在公共模式之外可以更轻松地进行升级。

基本备份和恢复

使用pg_dump命令行工具可以轻松备份完整的数据库，配合脚本轻松实现自动备份。

要备份nyc数据库，可以使用 GUI，只需右键单击要备份的数据库：

输入要创建的备份文件的名称。

请注意，有三种备份格式选项：compress、tar 和 plain。

• Plain：只是一个文本 SQL 文件。这是最简单的格式，在许多方面也是最灵活的，因为它可以轻松编辑或更改，然后加载回数据库，允许对所有权或其他全局信息等内容进行离线更改。
• Tar：使用 UNIX 归档格式将转储的组件保存在单独的文件中。使用 tar 格式允许pg_restore实用程序有选择地恢复部分转储。
• Compress：类似于 Tar 格式，但单独压缩内部组件，允许有选择地恢复它们而无需解压缩整个存档。

将检查压缩选项并保存备份文件。

可以使用命令行完成相同的操作：

pg_dump --file=nyc.backup --format=c --port=54321 --username=postgres nyc

由于备份文件是压缩格式，可以使用pg_restore命令列出清单来查看内容。在 PgAdmin GUI 中，查看是面板中的一个选项。

当查看清单时，您可能会注意到其中一件事是其中有很多FUNCTION签名。

这是因为pg_dump实用程序会转储数据库中的每个非系统对象，其中包括 PostGIS 函数定义。

PostgreSQL 9.1↑ 包含一个EXTENSION功能，允许像 PostGIS 这样的附加包作为注册的系统组件安装，因此从pg_dump输出中排除。PostGIS 2.0 及更高版本支持使用此扩展系统进行安装。

可以直接使用pg_restore从命令行查看相同的清单：

pg_restore --list nyc.backup

PostGIS 函数签名的转储文件的问题在于，真的想要转储数据，而不是系统函数。

由于每个对象都在转储文件中，可以恢复到空白数据库并获得全部功能。这样做时，要恢复到的系统具有与转储的系统完全相同的 PostGIS 版本（因为函数签名定义引用了 PostGIS 共享库的特定版本）。

从命令行恢复看起来像这样：

createdb --port 54321 nyc2
pg_restore --dbname=nyc2 --port 54321 --username=postgres nyc.backup

在没有函数签名的情况下，仅转储数据是在模式中包含数据的方便之处，因为有一个命令行标志可以仅转储特定模式：

pg_dump --port=54321 -format=c --schema=census --file=census.backup

当列出转储的内容时，只看到想要的数据表：

pg_restore --list census.backup

;
; Archive created at Thu Aug  9 11:02:49 2012
;     dbname: nyc
;     TOC Entries: 11
;     Compression: -1
;     Dump Version: 1.11-0
;     Format: CUSTOM
;     Integer: 4 bytes
;     Offset: 8 bytes
;     Dumped from database version: 8.4.9
;     Dumped by pg_dump version: 8.4.9
;
;
; Selected TOC Entries:
;
6; 2615 20091 SCHEMA - census postgres
146; 1259 19845 TABLE census nyc_census_blocks postgres
145; 1259 19843 SEQUENCE census nyc_census_blocks_gid_seq postgres
2691; 0 0 SEQUENCE OWNED BY census nyc_census_blocks_gid_seq postgres
2692; 0 0 SEQUENCE SET census nyc_census_blocks_gid_seq postgres
2681; 2604 19848 DEFAULT census gid postgres
2688; 0 19845 TABLE DATA census nyc_census_blocks postgres
2686; 2606 19853 CONSTRAINT census nyc_census_blocks_pkey postgres
2687; 1259 20078 INDEX census nyc_census_blocks_geom_gist postgres

只有数据表很方便，可以存储到安装了任何版本的 PostGIS 的数据库中。

备份用户

pg_dump实用程序一次操作一个数据库（或者一次操作一个模式或表）。但是，

有关用户的信息存储在整个集群中，而不是存储在任何一个数据库中，要备份用户信息，使用带有–globals-only标志的pg_dumpall实用程序。

pg_dumpall --globals-only --port 54321

还可以在默认模式下使用pg_dumpall来备份整个集群，但请注意，与pg_dump一样，最终将备份 PostGIS 函数签名，因此必须针对相同的软件安装恢复转储，它不能用作升级过程的一部分。

热备和恢复

热备和恢复允许管理员保持一组非常最新的备份文件，而无需重复转储整个数据库的开销。如果数据库承受频繁的插入和更新负载，那么热备可能比基本备份更可取。

PostgreSQL 热备详细说明请参考手册中有关连续归档和时间点恢复的相关部分。

工作机制

PostgreSQL 最初将更改存储在预写日志（WAL）中，而不是不断地写入主数据表。这些日志是对数据库所做的所有更改的完整记录。热备获取数据库主数据表的副本，然后获取从那时起生成的每个 WAL 的副本。

当需要恢复到新数据库时，系统会在主数据副本上启动，然后将所有 WAL 文件重做到数据库中。最终结果是恢复的数据库与上次接收 WAL 时的原始状态相同。

因为无论如何都在编写 WAL，并且将副本传输到存档服务器的计算成本很低，热备是一种保持系统非常最新的备份的有效方法，而无需求助于密集的定期完整转储。

归档 WAL 文件

设置热备的第一件事是创建归档方法。PostgreSQL 归档方法具有最大的灵活性：PostgreSQL 后端只需调用archive_command配置参数中指定的脚本。

这意味着存档可以像将文件复制到网络安装的驱动器一样简单，也可以像加密和通过电子邮件将文件发送到远程存档。

编辑postgresql.conf打开归档，首先打开 WAL 归档：

wal_level = archive
archive_mode = on

然后设置archive_command将存档文件复制到安全位置（根据需要更改目标路径）：

# Unix
archive_command = 'test ! -f /archivedir/%f && cp %p /archivedir/%f'

# Windows
archive_command = 'copy "%p" "C:\\archivedir\\%f"'

存档命令不要覆盖现有文件很重要，unix 命令包含一个初始测试以确保文件不存在。如果复制过程失败，该命令返回非零状态也很重要。

进行更改后，可以重新启动 PostgreSQL 以使其生效。

进行基础备份

归档过程到位后，需要进行基本备份。

将数据库置于备份模式（这不会改变查询或数据更新的操作，它只是强制检查点并写入一个标签文件，指示何时进行备份）。

SELECT pg_start_backup('/archivedir/basebackup.tgz');

对于标签，使用备份文件的路径是一种很好的做法，因为它可以帮助跟踪备份的存储位置。

将数据库复制到存档位置：

# Unix
tar cvfz /archivedir/basebackup.tgz ${PGDATA}

然后告诉数据库备份过程完成。

SELECT pg_stop_backup();

当然，所有这些步骤都可以为常规基础备份编写脚本。

从存档恢复

这些步骤取自 PostgreSQL连续归档和时间点恢复手册。

• 停止服务器，如果它正在运行。
• 如果您有空间这样做，请将整个集群数据目录和任何表空间复制到一个临时位置，以备日后需要它们时使用。请注意，此预防措施将要求您的系统上有足够的可用空间来保存现有数据库的两个副本。如果您没有足够的空间，您至少应该保存集群的 pg_xlog 子目录的内容，因为它可能包含在系统关闭之前未归档的日志。
• 删除集群数据目录下以及您正在使用的任何表空间的根目录下的所有现有文件和子目录。
• 从文件系统备份中恢复数据库文件。确保以正确的所有权（数据库系统用户，而不是 root！）和正确的权限恢复它们。如果您正在使用表空间，您应该验证 pg_tblspc/ 中的符号链接是否已正确恢复。
• 删除 pg_xlog/ 中存在的所有文件；这些来自文件系统备份，因此可能已过时而不是当前。如果您根本没有归档 pg_xlog/，则使用适当的权限重新创建它，如果您之前以这种方式设置它，请注意确保您将其重新建立为符号链接。
• 如果您在步骤 2 中保存了未归档的 WAL 段文件，请将它们复制到 pg_xlog/。（最好是复制它们，而不是移动它们，这样如果出现问题，你仍然有未修改的文件，你必须重新开始。）
• 在集群数据目录中创建恢复命令文件 recovery.conf（参见第 26 章）。您可能还想临时修改 pg_hba.conf 以防止普通用户连接，直到您确定恢复成功。
• 启动服务器。服务器将进入恢复模式并继续读取它需要的存档 WAL 文件。如果由于外部错误而终止恢复，则可以简单地重新启动服务器并继续恢复。恢复过程完成后，服务器会将 recovery.conf 重命名为 recovery.done（以防止以后意外重新进入恢复模式），然后开始正常的数据库操作。
• 检查数据库的内容以确保您已恢复到所需的状态。如果没有，请返回步骤 1。如果一切正常，请通过将 pg_hba.conf 恢复为正常来允许您的用户连接。

链接

• pg_dump
• pg_dumpall
• pg_restore
• PostgreSQL 高可用性
• PostgreSQL 高可用性持续归档和 PITR

升级

因为 PostGIS 驻留在 PostgreSQL 中，所以每个 PostGIS 安装实际上都包含两个版本的软件：PostgreSQL 版本和 PostGIS 版本。

作为一般原则，理论上每个版本的 PostGIS 都可以在多个版本的 PostgreSQL 中运行，反之亦然。

实际上，可用的确切版本对将由构建您的 PostgreSQL 发行版的打包者决定。大多数 Linux 软件包都为每个 PostgreSQL 版本版本包含几个 PostGIS 版本，允许根据您的偏好独立或同时升级这些部分。

升级可以考虑升级每个组件。

升级 PostgreSQL

PostgreSQL升级场景有两种：

• 当软件版本在补丁级别增加时的小升级。例如，从 8.4.3 到 8.4.4，或从 9.0.1 到 9.0.3。增加一个以上的补丁版本就好了。次要升级修复了错误，但不添加任何新功能或改变行为。
• 主要或次要版本增加时的主要升级。例如，从 8.4.5 到 9.0.0，或从 9.0.5 到 9.1.1。重大升级添加了新功能并改变了行为。

PostgreSQL 小升级

对于小升级，不需要特殊的过程，只需安装新软件，然后重新启动服务器。

主要 PostgreSQL 升级

对于重大升级，有两种方式进行升级。

转储/恢复

转储和恢复涉及在转储时将所有数据转换为平台中性格式（文本表示），并在恢复时转换回本机表示，因此它可能非常耗时且占用大量 CPU。但是，如果您要迁移到新的架构或操作系统，这是一个必需的过程。这也是一个经过时间考验且易于理解的升级路径，因此如果您的数据库不是太大，没有理由不坚持下去。

• pg_dumpall从旧数据库中转储数据。
• 安装新版本的 PostgreSQL 和您在旧数据库中使用的相同版本的 PostGIS。您需要匹配 PostGIS 版本，以便转储文件函数定义引用 PostGIS 库的预期版本。
• initdb使用新软件中的程序初始化新数据区。
• 在新数据区启动新服务器。
• 使用 恢复转储文件pg_restore。

升级

pg_upgrade实用程序允许升级 PostgreSQL 数据目录而不需要转储/恢复步骤。该实用程序不能处理对数据文件本身的更改，但会处理在 PostgreSQL 主要升级中发生的对系统表的更常见和频繁的更改。

运行升级过程的完整说明在 PostgreSQL 站点的pg_upgrade网页中。

pg_upgrade程序希望能够访问它正在使用的 PostgreSQL 的两个版本，旧版本和新版本，因此必须同时安装它们。

• 安装您将使用的新版 PostgreSQL。
• 将您在旧 PostgreSQL 中使用的相同版本的 PostGIS 安装到新 PostgreSQL 中。
• 用新的副本初始化新的 PostgreSQL 数据区initdb。
• 确保新旧 PostgreSQL 服务器均已关闭。
• 运行pg_upgrade，确保使用新软件安装中的二进制文件。

pg_upgrade
  --old-datadir "/var/lib/postgres/12/data"
  --new-datadir "/var/lib/postgres/13/data"
  --old-bindir "/usr/pgsql/12/bin"
  --new-bindir "/usr/pgsql/13/bin"

• 如果pg_upgrade生成了任何.sql文件，请立即运行它们。
• 启动新服务器。

升级 PostGIS

PostGIS 通过EXTENSION机制处理次要和升级。如果您使用空间启用数据库，则可以使用相同的功能更新数据库。

• 首先，安装新软件，使其可用于数据库；
• 然后，运行 SQL 以升级 PostGIS 扩展。

-- 如果您从 PostGIS 2.5 或更高版本升级并想要最新安装的版本
SELECT postgis_extensions_upgrade();

-- 如果您要从早期版本升级，则必须专门打开您想要的版本
ALTER EXTENSION postgis UPDATE TO '2.5.5';

高级几何构造

nyc_subway_stations表提供了许多有趣的示例，但其中有一些令人震惊的地方：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-a8xhZ9an-1646640988502)(https://postgis.net/workshops/postgis-intro/_images/adv_geom0.jpg)]

虽然它是一个包含所有车站的数据库，但它不允许简单的路线可视化！在本章中，我们将使用 PostgreSQL 和 PostGIS 的高级功能，从地铁站的点层构建一个新的线性路线层。

任务因两个问题而变得特别困难：

• nyc_subway_stations的routes列在每一行中有多个路线标识符，因此可能出现在多条路线中的车站在表中只出现一次。
• 车站表中没有路线排序信息，因此虽然可以找到特定路线中的所有车站，但无法使用属性来确定火车通过车站的顺序.

给定路线中的一组无序点，如何对它们进行排序以匹配实际路线。

以下是Q列车的站点

SELECT
	s.gid, s.geom
FROM nyc_subway_stations s
WHERE (strpos(s.routes, 'Q') <> 0);

在这张图片中，停靠点标有其唯一的gid主键。

如果从终点站一侧开始，则线路上的下一站似乎总是最近的。只要从搜索中排除所有先前找到的电台，就可以每次重复该过程。

有两种方法可以在数据库中运行这样的迭代例程：

• 使用过程语言，例如PL/PgSQL。
• 使用递归公用表表达式。

公用表表达式 (CTE) 的优点是不需要函数定义即可运行。这是计算Q列车路线的 CTE，从最北端（gid304 号）开始。

WITH RECURSIVE next_stop(geom, idlist) AS (
    (
        SELECT
          geom,
          ARRAY[gid] AS idlist
        FROM nyc_subway_stations
        WHERE gid = 304
    )
    UNION ALL
    (
        SELECT
          s.geom,
          array_append(n.idlist, s.gid) AS idlist
        FROM nyc_subway_stations s, next_stop n
        WHERE strpos(s.routes, 'Q') != 0
        AND NOT n.idlist @> ARRAY[s.gid]
        ORDER BY ST_Distance(n.geom, s.geom) ASC
        LIMIT 1
    )
)
SELECT geom, idlist FROM next_stop;

CTE 由结合在一起的两半组成：

• 前半部分表达式建立了一个起点。我们使用“gid”304的记录获取初始几何图形并初始化访问标识符数组。
• 后半部分迭代，直到找不到更多记录。在每次迭代中，它通过对next_stop的自引用获取上一次迭代的值。搜索 Q 线 (strpos(s.routes,'Q')) 上尚未添加到访问列表 ( NOT n.idlist @> ARRAY[s.gid]) 上的每个站点，并按距离对它们进行排序从上一点开始，只取第一个（最近的）。

除了递归 CTE 本身之外，这里还使用了许多高级 PostgreSQL 数组特性：

• 正在使用数组！PostgreSQL 支持任何类型的数组。在这种情况下，有一个整数数组，但我们也可以构建一个几何数组，或者任何其他 PostgreSQL 类型。
• 正在使用array_append来构建我们的访问标识符数组。
• 正在使用@>数组运算符（数组包含）来查找已经访问过哪些 Q 火车站。@>运算符在两边都需要 ARRAY 值，因此必须使用 ARRAY[] 语法将各个gid数字转换为单项数组。

运行查询时，会按照找到的顺序（即路线顺序）获取每个几何图形，以及已访问的标识符列表。将几何图形包装到 PostGIS ST_MakeLine聚合函数中，将几何图形集转换为单个线性输出，按提供的顺序构造。

WITH RECURSIVE next_stop(geom, idlist) AS (
    (
        SELECT
          geom,
          ARRAY[gid] AS idlist
        FROM nyc_subway_stations
        WHERE gid = 304
    )
    UNION ALL
    (
        SELECT
          s.geom,
          array_append(n.idlist, s.gid) AS idlist
        FROM nyc_subway_stations s, next_stop n
        WHERE strpos(s.routes, 'Q') != 0
        AND NOT n.idlist @> ARRAY[s.gid]
        ORDER BY ST_Distance(n.geom, s.geom) ASC
        LIMIT 1
    )
)
SELECT ST_MakeLine(geom) AS geom FROM next_stop;

看起来像这样：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eofs8Lpt-1646640988503)(https://postgis.net/workshops/postgis-intro/_images/adv_geom3.jpg)]

还有两个问题

• 这里只计算一条地铁路线，要计算所有路线。
• 查询包括一条_先验_知识，即作为构建路线的搜索算法的种子的初始站点标识符。

首先解决难题，找出路线上的第一个站点，而无需手动查看构成该路线的站点集。

Q列车站可以作为起点。路线的终点站有什么特点？

一个答案是它们是最北和最南的车站。然而，想象一下如果Q列车从东向西行驶。条件还能成立吗？

终点站的方向性较小的特征是它们是距离路线中间最远的站。有了这个特征，路线是南北还是东西都无关紧要，只是它或多或少是一个方向，特别是在末端。

由于没有 100% 的方法来确定终点，让我们尝试第二个排除方法。

“离中间最远”规则的一个明显失效模式是环形线，如英国伦敦的环形线。

要计算出每条路线的终点站，首先要弄清楚有哪些路线！找到了不同的路线。

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
)
SELECT * FROM routes;

注意两个高级 PostgreSQL ARRAY 函数的使用：

• string_to_array接收一个字符串并使用分隔符将其拆分为一个数组。PostgreSQL 支持任何类型的数组，因此可以构建字符串数组，就像本例一样，也可以构建几何和地理，我们将在本示例后面看到。
• unnest接收一个数组并为数组中的每个条目构建一个新行。其效果是将“水平”数组嵌入到单行中，并将其转换为“垂直”数组，每个值对应一行。

结果是所有唯一地铁路线标识符的列表。

 route
-------
 1
 2
 3
 4
 5
 6
 7
 A
 B
 C
 D
 E
 F
 G
 J
 L
 M
 N
 Q
 R
 S
 V
 W
 Z
(24 rows)

可以在此结果的基础上将其连接回nyc_subway_stations表以创建一个新表，其中对于每条路线，该路线上的每个站点都有一行。

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
)
SELECT * FROM stops;

 gid |                      geom                      | route
-----+----------------------------------------------------+-------
   2 | 010100002026690000CBE327F938CD21415EDBE1572D315141 | 1
   3 | 010100002026690000C676635D10CD2141A0ECDB6975305141 | 1
  20 | 010100002026690000AE59A3F82C132241D835BA14D1435141 | 1
  22 | 0101000020266900003495A303D615224116DA56527D445141 | 1
                            ...etc...

现在可以通过将每条路线的所有站点收集到一个多点中来找到中心点，并计算该多点的质心。

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
)
SELECT * FROM centers;

‘Q’ 列车站集合的中心点如下所示：

因此，最北端的站点（终点）似乎也是离中心最远的站点。让我们计算每条路线的最远点。

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
),
stops_distance AS (
  SELECT s.*, ST_Distance(s.geom, c.geom) AS distance
  FROM stops s JOIN centers c
  ON (s.route = c.route)
  ORDER BY route, distance DESC
),
first_stops AS (
  SELECT DISTINCT ON (route) stops_distance.*
  FROM stops_distance
)
SELECT * FROM first_stops;

这次添加了两个子查询：

• stop_distance将中心点连接回站表并计算每条路线的站和中心之间的距离。对结果进行排序，使得每条路线的记录成批出现，最远的站点作为批次的第一个记录。
• first_stops通过仅获取每个不同组的第一条记录来过滤stop_distance输出。由于我们订购stops_distance的方式，第一条记录是最远的记录，这意味着它是我们想要用作我们构建每条地铁路线的起始种子的车站。

现在知道每条路线，并且我们（大约）知道每条路线从哪个站点开始：已准备好生成路线线！

但首先，需要将递归 CTE 表达式转换为可以使用参数调用的函数。

CREATE OR REPLACE function walk_subway(integer, text) returns geometry AS
$$
WITH RECURSIVE next_stop(geom, idlist) AS (
    (SELECT
      geom AS geom,
      ARRAY[gid] AS idlist
    FROM nyc_subway_stations
    WHERE gid = $1)
    UNION ALL
    (SELECT
      s.geom AS geom,
      array_append(n.idlist, s.gid) AS idlist
    FROM nyc_subway_stations s, next_stop n
    WHERE strpos(s.routes, $2) != 0
    AND NOT n.idlist @> ARRAY[s.gid]
    ORDER BY ST_Distance(n.geom, s.geom) ASC
    LIMIT 1)
)
SELECT ST_MakeLine(geom) AS geom
FROM next_stop;
$$
language 'sql';

现在我们准备好了！

CREATE TABLE nyc_subway_lines AS
-- 不同的路线标识符
WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
-- Joined back to stops! Every route has all its stops!
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
-- 按路线收集站点并计算质心！
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
),
-- 计算每条路线中每个停靠点的停靠点/中心距离。
stops_distance AS (
  SELECT s.*, ST_Distance(s.geom, c.geom) AS distance
  FROM stops s JOIN centers c
  ON (s.route = c.route)
  ORDER BY route, distance DESC
),
-- 过滤掉最远的停止/中心对。
first_stops AS (
  SELECT DISTINCT ON (route) stops_distance.*
  FROM stops_distance
)
-- 将路线/停靠信息传递给线性路线生成函数！
SELECT
  ascii(route) AS gid, -- QGIS likes numeric primary keys
  route,
  walk_subway(gid, route) AS geom
FROM first_stops;

-- Do some housekeeping too
ALTER TABLE nyc_subway_lines ADD PRIMARY KEY (gid);

这是最终表格在 QGIS 中的可视化效果：

像往常一样，我们对数据的简单理解存在一些问题：

• 实际上有两列“S”（短途“穿梭”）列车，一列在曼哈顿，一列在洛克威，我们将它们连接在一起，因为它们都被称为“S”；
• ‘4’ 列车（和其他一些）在一条线的末端分成两个终点，因此“跟随一条线”假设中断，结果在末端有一个有趣的钩子。

希望这个示例提供了一些复杂数据操作的体验，这些操作可以结合 PostgreSQL 和 PostGIS 的高级功能。

另见

• PostgreSQL 数组
• PostgreSQL 数组函数
• PostgreSQL 递归通用表表达式
• PostGIS ST_MakeLine

附录

函数

构造

• ST_MakePoint(Longitude, Latitude)：返回一个新点。注意坐标的顺序（经度然后纬度）。
• ST_GeomFromText(WellKnownText, srid)：从标准 WKT 字符串和 srid 返回一个新几何。
• ST_SetSRID(geometry, srid)：更新几何上的 srid。返回相同的几何图形。这不会改变几何的坐标，它只会更新 srid。此功能对于调节没有 srid 创建的几何图形很有用。
• ST_Expand(geometry, Radius)：返回一个新几何，它是输入几何的扩展边界框。此功能对于创建用于索引搜索的信封很有用。

输出

• ST_AsText(geometry)：以人类可读的文本格式返回几何。
• ST_AsGML(geometry)：返回标准 OGC GML格式的几何图形。
• ST_AsGeoJSON(geometry)：将几何返回为标准GeoJSON格式。

测量

• ST_Area(geometry)：以空间参考系为单位返回几何图形的面积。
• ST_Length(geometry)：以空间参考系统的单位返回几何的长度。
• ST_Perimeter(geometry)：以空间参考系统为单位返回几何的周长。
• ST_NumPoints(linestring)：返回线串中的顶点数。
• ST_NumRings(polygon)：返回多边形中的环数。
• ST_NumGeometries(geometry)：返回几何集合中的几何数。

关系

• ST_Distance(geometry, geometry)：以空间参考系为单位返回两个几何图形之间的距离。
• ST_DWithin(geometry, geometry, radius)：如果几何在彼此的半径距离内，则返回 true，否则返回 false。
• ST_Intersects(geometry, geometry)：如果几何不是不相交的，则返回 true，否则返回 false。
• ST_Contains(geometry, geometry)：如果第一个几何完全包含第二个几何，则返回 true，否则返回 false。
• ST_Crosses(geometry, geometry)：如果一条线或多边形边界穿过另一条线或多边形边界，则返回 true，否则返回 false。

术语

CRS （coordinate reference system）

地理坐标系和投影坐标系的组合。

GDAL （Geospatial Data Abstraction Library）

一个开源光栅访问库，支持大量格式，广泛用于开源和专有软件。

GeoJSON （Javascript Object Notation）

一种在 Javascript 虚拟机中解析速度非常快的文本格式。在空间上，常用的是GeoJSON的扩展规范。

GIS （Geographic information system）

地理信息系统或地理信息系统捕获、存储、分析、管理和呈现与位置相关的数据。

GML （Geography Markup Language）

地理标记语言，GML 是用于表示空间特征信息的OGC标准 XML 格式。

JSON （Javascript Object Notation）

一种在 Javascript 虚拟机中解析速度非常快的文本格式。在空间上，常用的是GeoJSON的扩展规范。

JSTL （JavaServer Page Template Library）

一个用于 JSP 的标记库，它将 JSP中处理的许多标准功能（数据库查询、迭代、条件）封装成简洁的语法。

JSP （JavaServer Pages）

一种用于 Java 服务器应用程序的脚本系统，它允许标记和 Java 过程代码的交错。

KML （Keyhole Markup Language）

Keyhole 标记语言，谷歌地球使用的空间 XML 格式。Google 地球最初是由一家名为Keyhole的公司编写的，因此名称中包含（现在晦涩的）参考。

OGC （Open Geospatial Consortium）

开放地理空间联盟，是一个制定地理空间服务规范的标准组织。

OSGeo （Open Source Geospatial Foundation）

开源地理空间基金会，是一个致力于推广和支持开源地理空间软件的非营利基金会。

SFSQL （Simple Features for SQL）

OGC规范定义了构成标准空间数据库的类型和函数。

SLD （Styled Layer Descriptor）

OGC 规范定义了一种用于描述矢量特征的制图渲染的格式。

SRID （Spatial reference ID）

空间参考ID，分配给特定坐标参考系统的唯一编号。PostGIS 表spatial_ref_sys包含大量众所周知的 srid 值和坐标参考系统的文本表示形式。

SQL （Structured query language）

结构化查询语言，是查询关系数据库的标准手段。

SQL/MM （SQL Multimedia）

SQL多媒体，包括几个关于扩展类型的部分，包括关于空间类型的重要部分。

SVG （Scalable vector graphics）

一系列基于 XML 的文件格式规范，用于描述二维矢量图形，包括静态和动态（即交互式或动画）。

WFS （Web Feature Service）

OGC 的 Web 要素服务规范定义了一个用于在 Web 上读取和写入地理要素的接口。

WMS （Web Map Service）

OGC 的 Web 地图服务规范定义了一个接口，用于在 Web 上请求渲染的地图图像。

WKB （Well-known binary）

指的是 SQL 规范中描述的几何图形的二进制表示。

WKT （Well-known text）

可以引用几何图形的文本表示，以POINT、LINESTRING、POLYGON等开头的字符串。也可以引用 CRS 的文本表示，以PROJCS、GEOGCS等开头的字符串众所周知的文本表示是OGC标准，但没有自己的规范文档。WKT 的第一个描述（用于几何和 CRS）出现在 SFSQL 1.0 规范中。

其他链接

• PostGIS - http://postgis.net/
  • 文档 - http://postgis.net/docs/
• PostgreSQL - http://www.postgresql.org/
  • 文档 - https://www.postgresql.org/docs/
  • 下载 - https://www.postgresql.org/download/
  • JDBC 驱动程序 - https://jdbc.postgresql.org/
  • .Net 驱动程序 - https://www.npgsql.org/
  • Python 驱动程序 - http://www.pygresql.org/
  • C/C++ 驱动程序 - https://www.postgresql.org/docs/current/static/libpq.html
• PgAdmin - https://www.pgadmin.org/
• 开源桌面客户端
  • QGIS - https://qgis.org/
  • OpenJUMP - http://openjump.org/
  • uDig - http://udig.github.io/