30. 光栅¶
PostGIS 支持另一种称为光栅的空间数据类型。光栅数据与几何数据类似，使用笛卡尔坐标和空间参考系统。然而，光栅数据不是矢量数据，而是表示为由像素和波段组成的 n 维矩阵。波段定义了矩阵的数量。每个像素存储与每个波段相对应的值。因此，一个 3 波段光栅（例如 RGB 图像）将为每个像素存储 3 个值，分别对应于红-绿-蓝波段。
虽然像你在电视屏幕上看到的那些漂亮图片是光栅，但光栅可能并不那么令人兴奋。简而言之，光栅是一个固定在坐标系上的矩阵，它具有可以代表任何你想要它们代表的值。
由于光栅存在于笛卡尔空间中，因此光栅可以与几何体交互。PostGIS 提供了许多函数，这些函数以光栅和几何体作为输入。许多应用于光栅的操作将导致几何体。常见的有 ST_Polygon、ST_Envelope、ST_ConvexHull 和 ST_MinConvexHull，如下所示。当你将光栅转换为几何体时，输出的是光栅的 ST_ConvexHull。
光栅格式通常用于存储高程数据、温度数据、卫星数据和主题数据，这些数据代表诸如环境污染、人口密度和环境危害发生等事物。你可以使用光栅来存储任何具有有意义的坐标位置的数值数据。唯一的限制是，对于特定波段中的所有数据，数值数据类型必须相同。
虽然光栅数据可以在 PostGIS 中从头开始创建，但更常见的方法是使用与 PostGIS 一起打包的 raster2pgsql 命令行工具从各种格式加载光栅数据。在所有这些之前，你必须通过运行以下命令在你的数据库中启用光栅支持
CREATE EXTENSION postgis_raster;
30.1. 从几何体创建光栅¶
我们将首先从矢量数据创建栅格数据，然后继续从栅格源加载数据的更令人兴奋的方法。您会发现栅格数据非常丰富，而且通常可以从各种政府网站免费获得。
我们将首先使用 ST_AsRaster 函数将一些几何图形转换为栅格，如下所示。
CREATE TABLE rasters (name varchar, rast raster);
INSERT INTO rasters(name, rast)
SELECT f.word, ST_AsRaster(geom, width=>150, height=>150)
FROM (VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_Letters(word) AS geom;
CREATE INDEX ix_rasters_rast
ON rasters USING gist(ST_ConvexHull(rast));
上面的示例使用 PostGIS 3.2+ 的 ST_Letters 函数从字母形成的几何图形创建了一个表 (rasters)。与几何图形类似，栅格可以利用空间索引。用于栅格的空间索引是函数索引，它索引栅格的几何凸包。
您可以使用以下查询查看栅格的一些有用元数据，该查询利用 postgis 栅格的 ST_Count 函数来计算具有数据的像素数量，以及 ST_MetaData 函数来提供我们栅格的所有有用背景信息。
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
--------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
(2 rows)
注意
栅格函数有两个级别。有一些函数，例如 ST_MetaData，在栅格级别工作，还有一些函数，例如 ST_Count 函数和 ST_BandMetaData 函数，在波段级别工作。PostGIS 栅格中大多数在波段级别工作的函数一次只处理一个波段，并假设您想要的波段是 1。
如果您有一个多波段栅格，并且需要计算除 1 以外的波段中非空数据值的像素，则需要显式指定波段号，如下所示 ST_Count(rast,2)。
注意所有栅格都具有 150x150 的尺寸。这不是理想的。这意味着为了强制执行这一点，我们的栅格以各种方式被压缩。如果我们能看到我们面前的栅格的丑陋之处就好了。
30.2. 使用 raster2pgsql 加载栅格¶
raster2pgsql 是一个命令行工具，通常与 PostGIS 一起打包。如果您使用的是 Windows 并且使用了应用程序堆栈构建器 PostGIS Bundle，您将在 C:\Program Files\PostgreSQL\15\bin 文件夹中找到 raster2pgsql.exe，其中 15 应替换为正在运行的 PostgreSQL 版本。
如果您使用的是 Postgres.App，您将在其他 Postgres.app CLI 工具 中找到 raster2pgsql。
在 Ubuntu 和 Debian 上，您需要
apt install postgis
安装 PostGIS 命令行工具。这可能会安装额外的 PostgreSQL 版本。您可以使用 pg_lsclusters 命令查看 Debian/Ubuntu 中的集群列表，并使用 pg_dropcluster 命令删除它们。
对于本练习和后续练习，我们将使用 nyc_dem.tif，该文件位于 PG 栅格研讨会数据集 https://postgis.net.cn/stuff/workshop-data/postgis_raster_workshop.zip 中。对于一些几何/栅格示例，我们还将使用从前几章加载的 NYC 数据。您可以使用 pg_restore 命令行工具或 pgAdmin 的 还原 菜单，在您的数据库中还原 zip 文件中包含的 nyc_dem.backup，而不是加载 tif 文件。
注意
此栅格数据源自 NYC DEM 1-foot Integer，这是一个 3GB 的 DEM tif 文件，表示相对于海平面的高程，已去除建筑物和水面。然后我们创建了一个较低分辨率的版本。
命令 rasterpgsql 与 shp2gpsql 类似，不同之处在于它不是将 ESRI shapefile 加载到 PostGIS 几何/地理表中，而是将任何 GDAL 支持的栅格格式加载到栅格表中。就像 shp2pgsql 一样，您可以传递源数据的空间参考 ID (SRID)。与 shp2pgsql 不同的是，如果您的源数据具有合适的元数据，它可以推断源数据的空间参考系统。
有关所有可能开关的完整说明，请参阅 raster2pgsql 选项。
raster2pgsql 提供的其他一些值得注意的选项，我们不会介绍：
能够指定源数据的 SRID。相反，我们将依赖 raster2pgsql 的猜测能力。
能够设置 nodata 值，如果未指定，raster2pgsql 会尝试从文件中推断。
能够加载数据库外的栅格。
要加载我们文件夹中的所有 tif 文件并创建概览，我们将运行以下命令。
raster2pgsql -d -e -l 2,3 -I -C -M -F -Y -t 256x256 *.tif nyc_dem | psql -d nyc
-d 用于在表已存在时删除它们
上面的命令使用 -e 立即执行加载，而不是在事务中提交
-C 设置栅格约束，这对于 raster_columns 显示信息很有用。您可能希望将其与 -x 结合使用，以排除范围约束，这是一种检查速度很慢的约束，也会阻碍表中的未来加载。
-M 用于在加载后进行真空和分析，以改进查询计划统计信息
-Y 用于以 50 批的方式使用复制。如果您运行的是 PostGIS 3.3 或更高版本，可以使用 -Y 1000 使复制以 1000 批的方式进行，甚至更多。这将运行得更快，但会使用更多内存。
-l 2,3 用于创建概览表：o_2_ncy_dem 和 o_3_nyc_dem。这对于查看数据很有用。
-I 用于创建空间索引
-F 用于添加文件名，如果您只有一个 tif 文件，这有点没有意义。如果您有多个，这将有助于您了解每行来自哪个文件。
-t 用于设置块大小。请注意，如果您不确定最佳大小，请使用 -t auto，raster2pgql 将使用 tif 中的相同切片。输出将告诉您它选择的块大小。如果看起来很大或很奇怪，请取消。原始文件的大小为 300x7，这不是理想的。
使用 psql 在数据库上运行生成的 sql。如果您想将其转储到文件而不是数据库，请使用 > nyc_dem.sql
对于此示例，我们只有一个 tif 文件，因此我们可以指定完整的文件名，而不是 *.tif。如果文件不在当前目录中，您也可以使用 *.tif 指定文件夹路径。
注意
如果您使用的是 Windows 并且需要引用文件夹，请确保包含驱动器号，例如 C:/workshop/*.tif
您经常会在 PostGIS 行话中听到栅格切片和栅格这两个词，它们在某种程度上可以互换使用。栅格切片实际上对应于栅格列中的特定栅格，它是更大栅格的子集，例如我们刚刚加载的 NYC dem 数据。这是因为当从大型栅格文件加载栅格到 PostGIS 时，它们会被切成许多行以使其易于管理。然后，每行中的每个栅格都是更大栅格的一部分。每个切片覆盖由您指定的块大小表示的相同大小的区域。不幸的是，栅格受到 1GB PostgreSQL TOAST 限制以及缓慢的去吐司过程的限制，因此我们需要进行切片以实现良好的性能，甚至存储它们。
30.3. 在浏览器中查看栅格¶
虽然 pgAdmin 和 psql 目前还没有查看 PostGIS 光栅的机制，但我们有几个选择。对于较小的光栅，最简单的方法是使用内置的 PostGIS 光栅函数（如 ST_AsPNG 或 ST_AsGDALRaster，列在 PostGIS 光栅输出函数 中）输出到 Web 友好的光栅格式（如 PNG）。随着光栅越来越大，您可能需要使用 QGIS 等工具来查看它们的所有细节，或者使用 GDAL 系列命令行工具（如 gdal_translate）将它们导出到其他光栅格式。但请记住，PostGIS 光栅是为分析而构建的，而不是为了生成供您查看的漂亮图片。
需要注意的是，默认情况下所有不同光栅类型的输出都是禁用的。为了使用这些功能，您需要启用驱动程序，全部或部分，具体细节请参见 启用 GDAL 光栅驱动程序
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
如果您不想为每个连接都这样做，可以在数据库级别进行设置，使用
ALTER DATABASE nyc SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
每个新的数据库连接都将使用该设置。
运行以下查询，并将输出复制粘贴到 Web 浏览器的地址栏中。
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello';
对于迄今为止创建的光栅，我们没有指定波段数量，甚至没有定义它们与地球的关系。因此，我们的光栅具有未知的空间参考系统 (0)。
您可以将光栅的外骨骼视为几何图形。一个被几何外壳包裹的矩阵。为了进行有用的分析，我们需要对光栅进行地理配准，这意味着我们希望每个像素（矩形）代表一些有意义的空间区域。
ST_AsRaster 有许多重载表示形式。前面的示例使用了最简单的实现，并接受了默认参数，即 8BUI 和 1 个波段，无数据为 0。如果您需要使用其他变体，则应使用命名参数调用语法，这样您就不会意外地使用错误的函数变体或出现 **函数不唯一** 错误。
如果您从具有空间参考系统的几何图形开始，您将最终得到一个具有相同空间参考系统的光栅。在接下来的示例中，我们将把我们的文字放在纽约，并使用鲜艳的色彩。我们还将使用像素比例而不是宽度和高度，以便我们的光栅像素大小代表 1 米 x 1 米的空间。
INSERT INTO rasters(name, rast)
SELECT f.word || ' in New York' ,
ST_AsRaster(geom,
scalex => 1.0, scaley => -1.0,
pixeltype => ARRAY['8BUI', '8BUI', '8BUI'],
value => CASE WHEN word = 'Hello' THEN
ARRAY[10,10,100] ELSE ARRAY[10,100,10] END,
nodataval => ARRAY[0,0,0], gridx => NULL, gridy => NULL
) AS rast
FROM (
VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_SetSRID(
ST_Translate(ST_Letters(word),586467,4504725), 26918
) AS geom;
如果我们查看它，我们会看到一个没有压缩的彩色几何图形。
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello in New York';
重复光栅操作
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Raster in New York';
更重要的是，如果我们重新运行
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
观察纽约条目的元数据。它们具有纽约州平面米空间参考系统。它们也具有相同的比例。由于每个单位是 1x1 米，因此单词 **光栅** 的宽度现在比 **Hello** 更宽。
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
-------------------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+-------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
Hello in New York | 8786 | 586467 | 4504802.1 | 184 | 78 | 1 | -1 | 0 | 0 | 26918 | 3
Raster in New York | 10544 | 586467 | 4504800.4 | 258 | 76 | 1 | -1 | 0 | 0 | 26918 | 3
(4 rows)
30.4. 栅格空间目录表¶
与几何和地理类型类似，栅格也有一组目录，显示数据库中的所有栅格列。这些是 raster_columns 和 raster_overviews。
30.4.1. raster_columns¶
raster_columns 视图是 geometry_columns 和 geography_columns 的兄弟，提供几乎相同的数据，甚至更多，但针对的是栅格列。
SELECT *
FROM raster_columns;
探索该表，你会发现
r_table_catalog | r_table_schema | r_table_name | r_raster_column | srid | scale_x | scale_y | blocksize_x | blocksize_y | same_alignment | regular_blocking | num_bands | pixel_types | nodata_values | out_db | extent | spatial_index
----------------+----------------+--------------+-----------------+------+---------+---------+-------------+-------------+----------------+------------------+-----------+-------------+---------------+--------+--------+---------------
nyc | public | rasters | rast | 0 | | | | | f | f | | | | | | t
nyc | public | nyc_dem | rast | 2263 | 10 | -10 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_2_nyc_dem | rast | 2263 | 20 | -20 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_3_nyc_dem | rast | 2263 | 30 | -30 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
(4 rows)
对于 rasters 表，令人失望的是，大部分信息都是空白的。
与几何和地理不同，栅格不支持类型修饰符，因为类型修饰符的空间太有限，而关键属性比类型修饰符所能容纳的要多。
栅格依赖于约束，并作为视图的一部分读取这些约束。
查看我们使用 raster2pgsql 加载的表中的其他行。因为我们使用了 -C 开关，raster2pgsql 为 srid 和其他信息添加了约束，这些信息可以从 tif 文件中读取，或者是我们传入的。使用 -l 开关生成的概述表 o_2_nyc_dem 和 o_3_nyc_dem 也显示出来。
让我们尝试为我们的表添加一些约束。
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
你会收到很多关于栅格数据混乱且无法约束的通知。如果你再次查看 raster_columns，仍然是令人失望的，rasters 表中有很多空白行。
为了应用约束，表中的所有栅格必须至少符合一条规则。
我们可以尝试一下，假设所有数据都在纽约州平面坐标系中。
UPDATE rasters SET rast = ST_SetSRID(rast,26918)
WHERE ST_SRID(rast) <> 26918;
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
SELECT r_table_name AS t, r_raster_column AS c, srid,
blocksize_x AS bx, blocksize_y AS by, scale_x AS sx, scale_y AS sy,
ST_AsText(extent) AS e
FROM raster_columns
WHERE r_table_name = 'rasters';
有了进展
t | c | srid | bx | by | sx | sy | e
----------+------+-------+-----+-----+----+----+------------------------------------------
rasters | rast | 26918 | 150 | 150 | | | POLYGON((0 -0.90000000000..
(1 row)
你可以对所有栅格施加的约束越多，你就会看到填充的列越多，并且你能够对栅格中的所有瓦片执行的操作也越多。请记住，在某些情况下，你可能不想应用所有约束。
例如，如果您计划将更多数据加载到栅格表中，您将希望跳过范围约束，因为这将要求所有栅格都在范围约束的范围内。
30.4.2. raster_overviews¶
栅格概览列同时出现在 raster_columns 元目录和另一个名为 raster_overviews 的元目录中。概览主要用于加速更高缩放级别的查看。它们也可以用于快速估算分析，提供不太准确的统计数据，但速度比应用于原始栅格表快得多。
要检查概览，请运行
SELECT *
FROM raster_overviews;
您将看到以下输出
o_table_catalog | o_table_schema | o_table_name | o_raster_column | r_table_catalog | r_table_schema | r_table_name | r_raster_column | overview_factor
----------------+----------------+--------------+-----------------+-----------------+----------------+--------------+-----------------+-----------------
nyc | public | o_2_nyc_dem | rast | nyc | public | nyc_dem | rast | 2
nyc | public | o_3_nyc_dem | rast | nyc | public | nyc_dem | rast | 3
(2 rows)
raster_overviews 表只为您提供概览因子和父表的名称。所有这些信息都是您可以通过 raster2pgsql 概览命名约定自己推断出来的。
overview_factor 告诉您该行相对于其父级的分辨率。一个 overview_factor 为 2 表示 2x2 = 4 个瓦片可以放入一个 overview_2 瓦片中。类似地，一个 overview_factor 为 1 表示原始的 2x2x2 = 8 个瓦片可以塞入一个 overview_3 瓦片中。
30.5. 常见栅格函数¶
postgis_raster 扩展提供了超过 100 个函数可供选择。PostGIS 栅格功能是根据 PostGIS 几何支持模式化的。您会发现栅格和几何之间有重叠的函数，只要有意义。您将使用的一些常见函数在几何世界中具有等效函数，例如 ST_Intersects、ST_SetSRID、ST_SRID、ST_Union、ST_Intersection 和 ST_Transform。
除了这些重叠的函数之外，它还支持栅格之间以及栅格和几何之间的 && 重叠运算符。它还提供许多与几何一起使用或非常特定于栅格的函数。
您需要一个像 ST_Union 这样的函数来重建一个区域。由于性能会变慢，一个函数需要分析的像素越多，您需要一个快速执行的函数 ST_Clip 来将栅格裁剪到您分析中感兴趣的部分。
最后，您需要 ST_Intersects 或 && 来放大包含您感兴趣区域的栅格瓦片。 && 运算符比 ST_Intersects 更快。两者都可以利用栅格空间索引。我们将首先介绍这些基础函数，然后继续其他部分，在那里我们将它们与其他栅格和几何函数结合使用。
30.5.1. 使用 ST_Union 合并栅格¶
对于栅格的 ST_Union 函数，就像几何等效的 ST_Union 一样，将一组栅格聚合到单个栅格中。但是，与几何一样，并非所有栅格都可以组合在一起，但栅格合并的规则比几何规则更复杂。在几何的情况下，您只需要具有相同的空间参考系统，但对于栅格来说，这还不够。
如果您尝试以下操作
SELECT ST_Union(rast)
FROM rasters;
您将立即受到错误的惩罚
ERROR: rt_raster_from_two_rasters: The two rasters provided do not have the same alignment SQL state: XX000
什么是这种相同的对齐方式，它阻止您合并您宝贵的栅格？
为了将栅格组合在一起，它们需要在同一个网格上，也就是说，它们必须具有相同的像素大小、相同的方位（倾斜）、相同的空间参考系统，并且它们的像素不能相互切割，这意味着它们共享相同的全球像素网格。
如果您尝试相同的查询，但只是用我们仔细放置在纽约的词语。
同样，相同的错误。这些是相同的空间参考系统，相同的像素大小，但它仍然不够好。因为它们的网格是关闭的。
我们可以通过稍微调整左上角的 y 坐标来解决这个问题，然后再次尝试。如果我们的网格从整数级别开始，因为我们的像素大小是整数，那么像素就不会相互切割。
UPDATE rasters SET rast = ST_SetUpperLeft(rast,
ST_UpperLeftX(rast)::integer,
ST_UpperLeftY(rast)::integer)
WHERE name LIKE '%New York';
SELECT ST_Union(rast ORDER BY name)
FROM rasters
WHERE name LIKE '%New York%';
瞧，它起作用了，如果我们要查看，我们会看到类似的东西
注意
如果您不清楚为什么您的栅格没有相同的对齐方式，您可以使用函数 ST_SameAlignment，它将比较两个栅格或一组栅格，并告诉您它们是否具有相同的对齐方式。如果您启用了通知，则通知将告诉您问题栅格的错误之处。 ST_NotSameAlignmentReason，而不是仅仅通知，将输出原因。但是，它一次只能处理两个栅格。
与 ST_Union(raster) 光栅函数不同，ST_Union(geometry) 几何函数允许使用名为 uniontype 的参数。默认情况下，如果未指定该参数，则将其设置为 LAST，这意味着在光栅像素值重叠的情况下，取 **最后一个** 光栅像素值。一般来说，标记为无数据的波段中的像素将被忽略。
与大多数 PostgreSQL 聚合函数一样，您可以在函数调用中添加 ORDER BY 子句，如前面的示例所示。指定顺序可以控制哪个光栅优先。因此，在前面的示例中，Raster 优先于 Hello，因为 Raster 在字母顺序上排在最后。
观察一下，如果您改变顺序
SELECT ST_Union(rast ORDER BY name DESC)
FROM rasters
WHERE name LIKE '%New York%';
那么 Hello 优先于 Raster，因为 Hello 现在是最后一个叠加的。
FIRST 并集类型与 LAST 相反。
但在某些情况下，LAST 可能不是正确的操作。假设我们的光栅代表来自两个不同设备的两组不同的观测值。这些设备测量的是同一件事，我们不确定它们交叉时哪个是正确的，因此我们希望取结果的 MEAN。我们将这样做
SELECT ST_Union(rast, 'MEAN')
FROM rasters
WHERE name LIKE '%New York%';
瞧，它起作用了，如果我们要查看，我们会看到类似的东西
因此，我们不是取优先级，而是将两种力量混合在一起。在 MEAN 并集类型的情况下，指定顺序没有意义，因为结果将是重叠像素值的平均值。
请注意，对于几何体来说，由于几何体是矢量，因此除了存在或不存在之外没有其他值，因此在两个矢量相交时，如何组合它们实际上没有歧义。
光栅 ST_Union 的另一个我们略过的特性是，您应该返回所有波段还是只返回一些波段。当您没有指定要合并的波段时，ST_Union 将合并相同波段编号并使用 LAST 并集策略。如果您有多个波段，这可能不是您想要做的事情。也许您只想合并第二个波段。在这种情况下，是绿色波段，您想要像素值的计数。
SELECT ST_BandPixelType(ST_Union(rast, 2, 'COUNT'))
FROM rasters
WHERE name LIKE '%New York%';
st_bandpixeltype
------------------
32BUI
(1 row)
请注意，对于COUNT联合类型，它计算填充的像素数量并返回该值，结果始终为32BUI，类似于在 SQL 中执行COUNT时，结果始终为bigint，以适应较大的计数。
在其他情况下，波段像素类型不会改变，并设置为最大值或四舍五入，如果数量超过类型的范围。为什么有人会想要计算在某个位置相交的像素呢？假设您的每个栅格都代表警察中队和该地区的逮捕事件。每个值可能代表不同的逮捕原因。您正在统计每个区域的逮捕次数，因此您只关心逮捕次数。
或者，您可能想对所有波段进行操作，但您想要不同的策略。
SELECT ST_Union(rast, ARRAY[(1, 'MAX'),
(2, 'MEAN'),
(3, 'RANGE')]::unionarg[])
FROM rasters
WHERE name LIKE '%New York%';
使用ST_Union函数的unionarg[]变体，还可以重新排列波段的顺序。
30.5.2. 使用 ST_Intersects 剪切栅格¶
ST_Clip函数是 PostGIS 栅格中最常用的函数之一。主要原因是，您需要检查或操作的像素越多，处理速度就越慢。ST_Clip 将栅格剪切到感兴趣的区域，因此您可以将操作隔离到该区域。
此函数也很特殊，因为它利用了几何的力量来帮助栅格分析。为了减少像素数量，ST_Union 必须首先将每个栅格剪切到我们感兴趣的区域。
SELECT ST_Union( ST_Clip(r.rast, g.geom) )
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
此示例展示了几个函数协同工作。所使用的ST_Intersects函数是与postgis_raster一起打包的，可以对两个栅格或一个栅格和一个几何体进行相交。类似于几何ST_Intersects，栅格 ST_Intersects 可以利用栅格或几何表上的空间索引。
30.6. 将栅格转换为几何体¶
栅格可以很容易地转换为几何体。
30.6.1. 使用 ST_Polygon 获取栅格的多边形¶
让我们从之前的示例开始，但使用ST_Polygon函数将其转换为多边形。
SELECT ST_Polygon(ST_Union( ST_Clip(r.rast, g.geom) ))
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
如果您单击 pgAdmin 中的几何体查看器，您可以在没有任何黑客的情况下看到它。
ST_Polygon 会考虑特定波段中所有具有值（非无数据）的像素，并将它们转换为几何图形。与栅格中的许多其他函数一样，ST_Polygon 仅考虑 1 个波段。如果未指定波段，它将仅考虑第一个波段。
30.6.2. 使用 ST_PixelAsPolygons 获取栅格的像素矩形¶
另一个常用的函数是 ST_PixelAsPolygons 函数。您应该很少在大型栅格上使用 ST_PixelAsPolygons，因为您最终会得到数百万行，每行对应一个像素。
ST_PixelAsPolygons 返回一个包含 geom、val、x 和 y 的表。其中 x 是列号，y 是栅格中的行号。
ST_PixelAsPolygons 与其他栅格函数类似，一次处理一个波段，如果未指定波段，则处理波段 1。它还默认仅返回具有值的像素。
SELECT gv.*
FROM rasters AS r
CROSS JOIN LATERAL ST_PixelAsPolygons(rast) AS gv
WHERE r.name LIKE '%New York%'
LIMIT 10;
输出结果为
如果我们使用几何查看器进行检查，我们会看到
如果我们想要所有波段的所有像素，我们需要执行以下操作。请注意此示例与之前示例的不同之处。
1. 设置 exclude_nodata_value 以确保返回所有像素，以便我们的调用集返回相同数量的行。函数输出的行将自然按相同顺序排列。
2. 使用 PostgreSQL ROWS FROM 构造函数，并使用名称为函数输出中的每一组列取别名。例如，波段 1 列（geom、val、x、y）分别重命名为 g1、v1、x1、x2。
SELECT pp.g1, pp.v1, pp.v2, pp.v3
FROM rasters AS r
CROSS JOIN LATERAL
ROWS FROM (
ST_PixelAsPolygons(rast, 1, exclude_nodata_value => false ),
ST_PixelAsPolygons(rast, 2, exclude_nodata_value => false),
ST_PixelAsPolygons(rast, 3, exclude_nodata_value => false )
) AS pp(g1, v1, x1, y1,
g2, v2, x2, y2,
g3, v3, x3, y3 )
WHERE r.name LIKE '%New York%'
AND ( pp.v1 = 0 OR pp.v2 > 0 OR pp.v3 > 0) ;
注意
我们在这些示例中使用 CROSS JOIN LATERAL 是因为我们希望明确说明我们正在做什么。由于这些都是返回集的函数，您可以用 , 代替 CROSS JOIN LATERAL，以简化操作。在下一组示例中，我们将使用 ,
30.6.3. 使用 ST_DumpAsPolygons 导出多边形¶
栅格还引入了另一种名为 geomval 的复合类型。可以将 geomval 视为几何图形和栅格的子类。它包含一个几何图形，也包含一个像素值。
你会发现一些返回 geomvals 的栅格函数。
一个常用的输出 geomvals 的函数是 ST_DumpAsPolygons，它返回一组具有相同值的连续像素，作为多边形。同样，默认情况下，它只会检查波段 1 并排除无数据值，除非你覆盖它。此示例仅选择波段 2 中的多边形。你也可以对值应用过滤器。对于大多数用例，ST_DumpAsPolygons 比 ST_PixelAsPolygons 更好的选择，因为它将返回更少的行。
这将输出 6 行，并返回对应于“Raster”中字母的多边形。
SELECT gv.geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
AND gv.val = 100;
请注意，它不会返回单个几何图形，因为它会找到具有相同值的连续像素集，这些像素集形成一个多边形。即使所有这些值都相同，它们也不连续。
一种常见的生成更复杂几何图形的方法是按值分组并合并。
SELECT ST_Union(gv.geom) AS geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
GROUP BY gv.val;
这将给你返回 2 行，对应于“Raster”和“Hello”这两个词。
30.7. 统计信息¶
关于栅格最重要的理解是，它们是用于在数组中存储数据的统计工具，你可能碰巧能够在屏幕上使其看起来很漂亮。
你可以在 栅格波段统计信息 中找到这些统计函数的菜单。
30.7.1. ST_SummaryStatsAgg 和 ST_SummaryStats¶
想要获取一组栅格的所有统计信息，请使用 ST_SummaryStatsAgg 函数。
此查询大约需要 10 秒，并提供整个表的摘要
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem;
输出
count | sum | mean | stddev | min | max
-----------+------------+------------------+------------------+-----+-----
246794100 | 4555256024 | 18.4577184948911 | 39.4416860598687 | 0 | 411
(1 row)
这告诉我们有很多像素，我们的最大海拔是 411 英尺。
如果你已经构建了概览，并且只需要对最小值、最大值和平均值进行粗略估计，请使用其中一个概览。下一个查询返回与先前查询大致相同的最小值、最大值和平均值，但大约需要 1 秒而不是 10 秒。
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem ;
现在有了这些信息，我们可以提出更多问题。
30.7.2. ST_Histogram¶
通常，您不希望获得整个表格的统计信息，而只希望获得特定区域的统计信息。在这种情况下，您需要使用我们老朋友 ST_Intersects 和 ST_Clip。如果您还需要一个没有聚合版本的栅格统计函数，那么您需要将 ST_Union 作为辅助函数。
在下一个示例中，我们将使用一个不同的统计函数 ST_Histogram，它没有聚合等效项，并且对于此特定变体，是一个返回集合的函数。我们使用与之前示例相同的感兴趣区域，但我们还需要使用几何 ST_Transform 将我们的纽约州平面米几何转换为纽约市州平面英尺栅格。将几何而不是栅格进行转换几乎总是更高效的，尤其是当您的几何只是一个时。
SELECT (ST_Quantile( ST_Union( ST_Clip(r.rast, g.geom) ), ARRAY[0.25,0.50,0.75, 1.0] )).*
FROM nyc_dem AS r
INNER JOIN
ST_Transform(
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ),
2263) AS g(geom)
ON ST_Intersects(r.rast, g.geom);
上面的查询在 60 毫秒内完成并输出
quantile | value
----------+-------
0.25 | 52
0.5 | 57
0.75 | 68
1 | 78
(4 rows)
30.8. 创建派生栅格¶
PostGIS 栅格包含许多用于编辑栅格的函数。这些函数既用于编辑，也用于创建派生栅格数据集。您可以在 栅格编辑器 和 栅格管理 中找到这些函数。
30.8.1. 使用 ST_Transform 转换栅格¶
我们的大多数数据都位于纽约州平面米（SRID：26918）中，但我们的 DEM 栅格数据集位于纽约州平面英尺（SRID：2263）中。为了获得最简单的流程，我们需要我们的核心数据集位于相同的空间参考系统中。
栅格 ST_Transform 是最适合此工作的函数。
为了在纽约州平面米中创建一个新的 nyc dem 数据集，我们将执行以下操作
CREATE TABLE nyc_dem_26918 AS
WITH ref AS (SELECT ST_Transform(rast,26918) AS rast
FROM nyc_dem LIMIT 1)
SELECT r.rid, ST_Transform(r.rast, ref.rast) AS rast, r.filename
FROM nyc_dem AS r, ref;
在我的系统上，上述操作大约需要 1.5 分钟。对于更大的数据集，它将花费更长时间。
上述示例使用了 ST_Transform 栅格函数的两种变体。第一个是获取一个参考栅格，该栅格将用于将其他栅格瓦片进行转换，以确保所有瓦片具有相同的对齐方式。请注意，使用的 ST_Transform 的第二个变体甚至不接受输入 SRID。这是因为 SRID 和所有像素比例和块大小都从参考栅格中读取。如果您使用的是 ST_Transform(rast, srid) 形式，那么您所有的栅格可能会以不同的对齐方式出现，这使得无法对它们应用诸如 ST_Union 之类的操作。
前面提到的 ST_Transform 方法的唯一问题是，当您进行转换时，转换后的结果通常会存在于其他瓦片中。如果您仔细观察上面的输出，通过输出栅格的凸包，在下一个示例中您会看到边界周围令人讨厌的重叠。
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
在 pgAdmin 中查看会类似于
30.8.2. 使用 ST_MakeEmptyCoverage 创建均匀的瓦片栅格¶
一个更好的方法，虽然速度稍慢，是从头开始使用 ST_MakeEmptyCoverage 定义您自己的覆盖瓦片结构，然后找到每个新瓦片的相交瓦片，并对这些瓦片进行 ST_Union，然后使用 ST_Transform(ref, ST_Union…) 创建每个瓦片。
为此，我们将使用之前学过的一些函数。
DROP TABLE IF EXISTS nyc_dem_26918;
CREATE TABLE nyc_dem_26918 AS
SELECT ROW_NUMBER() OVER(ORDER BY t.rast::geometry) AS rid,
ST_Union(ST_Clip( ST_Transform( r.rast, t.rast, 'Bilinear' ), t.rast::geometry ), 'MAX') AS rast
FROM (SELECT ST_Transform(
ST_SetSRID(ST_Extent(rast::geometry),2263)
, 26918) AS geom
FROM nyc_dem
) AS g, ST_MakeEmptyCoverage(tilewidth => 256, tileheight => 256,
width => (ST_XMax(g.geom) - ST_XMin(g.geom))::integer,
height => (ST_YMax(g.geom) - ST_YMin(g.geom))::integer,
upperleftx => ST_XMin(g.geom),
upperlefty => ST_YMax(g.geom),
scalex => 3.048,
scaley => -3.048,
skewx => 0., skewy => 0., srid => 26918) AS t(rast)
INNER JOIN nyc_dem AS r
ON ST_Transform(t.rast::geometry, 2263) && r.rast
GROUP BY t.rast;
重复与之前相同的练习
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
在 pgAdmin 中查看，我们不再有重叠
在我的系统上，这大约花了 10 分钟，返回了 3879 行。在创建表之后，我们希望像往常一样添加空间索引、主键和约束，如下所示
ALTER TABLE nyc_dem_26918
ADD CONSTRAINT pk_nyc_dem_26918 PRIMARY KEY(rid);
CREATE INDEX ix_nyc_dem_26918_st_convexhull_gist
ON nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
SELECT AddRasterConstraints('nyc_dem_26918'::name, 'rast'::name);
ANALYZE nyc_dem_26918;
对于这个数据集，这应该在 2 分钟内完成。
30.8.3. 使用 ST_CreateOverview 创建概览表¶
就像我们的原始数据集一样，拥有概览表将有助于提高某些操作的性能。 ST_CreateOverview 是一个适合此目的的函数。如果您在 raster2pgsql 加载过程中忽略了创建概览表，或者您决定需要更多概览表，也可以使用 ST_CreateOverview 来创建概览表。
我们将创建级别 2 和 3 的概览表，就像我们使用以下代码对原始数据集所做的那样。
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 2);
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 3);
不幸的是，这个过程需要一段时间，而且行数越多，时间越长，所以请耐心等待。对于这个数据集，概览因子 2 大约花了 3-5 分钟，概览因子 3 大约花了 1 分钟。
ST_CreateOverView 函数还会添加必要的约束，以便这些列在 raster_columns 和 raster_overviews 目录中以完整的详细信息显示。但是，它不会向它们添加索引，也不会添加 rid 列。除非您需要主键来进行编辑，否则 rid 列可能不是必需的。您可能需要一个索引，您可以使用以下命令创建它
CREATE INDEX ix_o_2_nyc_dem_26918_st_convexhull_gist
ON o_2_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
CREATE INDEX ix_o_3_nyc_dem_26918_st_convexhull_gist
ON o_3_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
注意
ST_CreateOverview 有一个可选参数用于指定采样方法。如果未指定，它将使用默认的 NearestNeighbor 方法，该方法通常计算速度最快，但可能不是理想的。重采样方法超出了本研讨会的范围。
30.9. 栅格与几何的交集¶
通常使用几个函数来计算栅格与几何的交集。我们已经看到了 ST_Clip 的作用，它将栅格与几何的交集作为栅格返回，但还有其他函数。对于点数据，最常用的函数是 ST_Value，然后是 ST_Intersection，它有几个重载，有些返回栅格，有些返回一组 geomval。
30.9.1. 几何点处的像素值¶
如果您需要根据几个临时几何点的交集从栅格中返回值，您将使用 ST_Value 或其最近的亲属 ST_NearestValue。
SELECT g.geom, ST_Value(r.rast, g.geom) AS elev
FROM nyc_dem_26918 AS r
INNER JOIN
(SELECT id, geom
FROM nyc_homicides
WHERE weapon = 'gun') AS g
ON r.rast && g.geom;
此示例大约需要 1 秒才能返回 2444 行。如果您使用 ST_Intersects 而不是 &&，则该过程将大约需要 3 秒。 ST_Intersects 速度较慢的原因是在某些情况下它会执行额外的重新检查，即逐像素检查。如果您希望所有点都用栅格集中表示，并且您的栅格表示覆盖范围（一组连续的非重叠栅格瓦片），那么 && 通常是一个更快的选择。
如果您的栅格数据没有密集填充，或者您有重叠的栅格（例如，它们代表不同时间段的观测值），或者它们是倾斜的（不是轴对齐的），那么使用 ST_Intersects 来剔除误报是有优势的。
30.9.2. ST_Intersection 栅格样式¶
就像您可以使用 ST_Intersection 计算两个几何的交集一样，您也可以使用 栅格 ST_Intersection 计算两个栅格或栅格与几何的交集。
从这个函数中，您可以得到两种不同类型的结果：
将几何与栅格相交，您将得到一组 geomval 子项。可能只有一个，但通常会有很多。
将两个栅格相交，您将得到一个 raster。
栅格交集和几何交集的黄金法则是在两个参与方都必须具有相同的空间参考系。对于栅格/栅格，它们还必须具有相同的对齐方式。
这里有一个例子，可以回答你可能好奇的问题。如果我们将海拔高度分成 5 个海拔值范围，哪个海拔范围会导致最多的枪支致死事件？根据我们之前的汇总统计数据，我们知道 0 是我们 nyc dem 数据集中海拔的最低值，而 411 是最高值，因此我们将它们用作 width_bucket 函数的最小值和最大值。
SELECT ST_Transform(ST_Union(gv.geom),4326) AS geom ,
MIN(gv.val) AS min_elev, MAX(gv.val) AS max_elev,
count(g.id) AS count_guns
FROM nyc_dem_26918 AS r
INNER JOIN nyc_homicides AS g
ON ST_Intersects(r.rast, g.geom)
CROSS JOIN
ST_Intersection( g.geom,
ST_Clip(r.rast,ST_Expand(g.geom, 4) )
) AS gv
WHERE g.weapon = 'gun'
GROUP BY width_bucket(gv.val, 0, 411, 5)
ORDER BY width_bucket(gv.val, 0, 411, 5);
枪支凶杀案和海拔高度之间是否存在重要关联？可能没有。
让我们看一下栅格/栅格交集。
SELECT ST_Intersection(r1.rast, 1, r2.rast, 1, 'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
我们得到两行包含 NULL 值，如果你将 PostgreSQL 设置为显示通知，你将看到
NOTICE: 提供的两个栅格没有相同的对齐方式。返回 NULL
为了解决这个问题，我们可以使用 ST_Resample 在栅格出厂时将其中一个对齐到另一个。
SELECT ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1,
'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
我们也将其汇总成一个单独的统计记录。
SELECT (
ST_SummaryStatsAgg(
ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1, 'BAND1'),
1, true)
).*
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
输出结果为
count | sum | mean | stddev | min | max
-------+-------+-----------------+------------------+-----+-----
2075 | 99536 | 47.969156626506 | 9.57974836865737 | 33 | 62
(1 row)
30.10. 地图代数函数¶
地图代数是指你可以对像素值进行数学运算的理念。之前介绍的 ST_Union 和 ST_Intersection 函数是地图代数的一种特殊快速情况。然后是 ST_MapAlgebra 函数系列，它允许你定义自己的疯狂数学运算，但会影响性能。
人们习惯于直接使用 ST_MapAlgebra，可能是因为这个名字听起来很酷很复杂。谁不想告诉他们的朋友，“我正在使用一个叫做 ST_MapAlgebra 的函数”。我的建议是，在你使用这个函数之前，先探索其他函数。你的生活会更简单，性能会提高 100 倍，代码也会更短。
在我们展示 ST_MapAlgebra 之前，我们将探索其他属于 地图代数 函数系列的函数，这些函数通常比 ST_MapAlgebra 性能更好。
30.10.1. 使用 ST_Reclass 对栅格进行重分类¶
一个经常被忽视的地图代数函数是 ST_Reclass 函数，它默默无闻地等待着有人发现它强大的功能和速度。
ST_Reclass 做什么？顾名思义，它根据最小范围代数对像素值进行重分类。
让我们重新审视我们的 NYC Dems。也许我们只关心将海拔高度分类为 1) 低、2) 中等、3) 高和 4) 非常高。我们不需要 411 个值，只需要 4 个。因此，让我们进行一些重分类。
分类方案由 重分类表达式 控制。
WITH r AS ( SELECT ST_Union(newrast) As rast
FROM nyc_dem_26918 AS r
INNER JOIN ST_Buffer(ST_Point(586598, 4504816, 26918), 1000 ) AS g(geom)
ON ST_Intersects( r.rast, g.geom )
CROSS JOIN ST_Reclass( ST_Clip(r.rast,g.geom), 1,
'[0-10):1, [10-50):2, [50-100):3,[100-:4','4BUI',0) AS newrast
)
SELECT SUM(ST_Area(gv.geom)::numeric(10,2)) AS area, gv.val
FROM r, ST_DumpAsPolygons(rast) AS gv
GROUP BY gv.val
ORDER BY gv.val;
这将输出
area | val
------------+-----
6754.04 | 1
1753117.51 | 2
1355232.37 | 3
1848.75 | 4
(4 rows)
如果这是我们喜欢的分类方案，我们可以使用 ST_Reclass 创建一个新表来重新计算每个图块。
30.10.2. 使用 ST_ColorMap 为栅格着色¶
ST_ColorMap 函数是另一个地图代数函数，它可以重新分类像素值。不同的是，它是创建波段的。它将单波段栅格（如我们的 DEM）转换为视觉上可呈现的 3 或 4 波段栅格。
如果您不想费心创建颜色映射，可以使用以下内置颜色映射之一。
SELECT ST_ColorMap( ST_Union(newrast), 'bluered') As rast
FROM nyc_dem_26918 AS r
INNER JOIN
ST_Buffer(
ST_Point(586598, 4504816, 26918), 1000
) AS g(geom)
ON ST_Intersects( r.rast, g.geom)
CROSS JOIN ST_Clip(rast, g.geom) AS newrast;
看起来像
颜色越蓝，海拔越低；颜色越红，海拔越高。
