前几章：https://www.cnblogs.com/2008nmj/p/14839911.html
第五章：使用光栅数据
在本章中，我们将介绍以下内容：
- 获取和加载光栅
- 使用基本光栅信息和分析
- 执行简单的映射代数操作
- 将几何图形与光栅结合起来进行分析
- 在光栅和几何图形之间转换
- 用GDAL-VRT处理和加载光栅
- 扭曲和重采样光栅
- 执行高级地图代数操作
- 执行DEM操作
- 通过SQL共享和可视化光栅
简介
在本章中，配方以一步一步的工作流程呈现，您可以在使用光栅时应用该流程。这需要加载光栅，对光栅有基本的了解，对其进行处理和分析，并将其交付给消费者。我们有意在工作流程中添加一些迂回，以反映这样一个现实：光栅的原始形式可能令人困惑，不适合分析。在本章结束时，你应该能够从食谱中吸取教训，并自信地应用它们来解决你的问题。
（栅格数据信息冗余：3601×3601=12967201个像素，如果要进行查询，，，某个坐标所在的位置的属性信息，例如高程、温度等等）
在进一步讨论之前，我们应该描述一下光栅是什么，以及光栅的用途。在最简单的层面上，光栅是一张照片或图像，其中包含描述将光栅放置在地球表面何处的信息。照片通常有三组值：每种原色（红、绿、蓝）一组。光栅也有一组值，通常比照片中的值更多。每组值称为一个频带。因此，一张照片通常有三个波段，而光栅至少有一个波段。像数码照片一样，光栅有各种各样的文件格式。常见的光栅格式包括PNG、JPEG、GeoTIFF、HDF5和NetCDF。由于光栅可以有许多条带，甚至更多的值，因此可以使用它们高效地存储大量数据。由于它们的效率，光栅被用于卫星和航空传感器和模型表面，如天气预报。
本章和PostGIS生态系统中使用了一些需要定义的关键字：
光栅(Raster)：这是用于在PostgreSQL中存储光栅文件的PostGIS数据类型。
切片(Tile)：这是原始光栅文件的一小块，存储在表行的一列中。每个图块都有自己的空间信息集，因此独立于同一表的同一列中的所有其他图块，即使其他图块来自同一原始光栅文件。
覆盖率(Coverage)：它由一个表中单个光栅列的所有平铺组成。
在本章中我们大量使用GDAL。GDAL通常被认为是事实上与光栅工作的瑞士军刀。GDAL不是一个单一的应用程序，而是一个光栅抽象库，包含许多有用的实用程序。通过GDAL，您可以获得光栅的元数据，将该光栅转换为其他格式，并在许多其他功能中扭曲该光栅。为了满足本章的需要，我们将使用三个GDAL实用程序：gdalinfo、gdalbuildvrt和gdal_translate。（栅格数据查询）（ArcGIS Server 中栅格数据的查询）（栅格索引统计）
5.1 Getting and loading rasters
在这个配方中，我们加载了本章中使用的大多数光栅。这些光栅是卫星图像和模型生成曲面的示例，这是两种最常见的光栅源。
5.1.1 Getting ready
如果您还没有这样做，创建一个目录并复制章节的数据集；对于Windows，请使用以下命令：
>>psql -d postgis_sample -c "CREATE SCHEMA chp05" -U postgres
5.1.2 How to do it...
我们将从2016年PRISM月平均最低气温光栅数据集开始，覆盖美国大陆。光栅由俄勒冈州立大学棱镜气候组提供，附加光栅可供使用。在：http://www.prism.oregonstate.edu/mtd/。
在命令行上，导航到PRISM目录，如下所示：
>>cd PRISM
让我们用GDAL实用程序gdalinfo抽查一个棱镜光栅。最好至少检查一个光栅以了解元数据并确保光栅没有任何问题。可以使用以下命令执行此操作：
>>gdalinfo PRISM_tmin_provisional_4kmM2_201703_asc.asc
gdalinfo输出如下：
gdalinfo输出显示光栅没有任何问题，这可以通过角坐标、像素大小、波段和坐标系为非空来证明。
通过对元数据的查阅，发现空间参照系统的元数据表明光栅使用NAD83坐标系。我们可以通过在spatial_ref_sys表中搜索NAD83的详细信息，对此进行双重检查：
SELECT srid, auth_name, auth_srid, srtext, proj4text FROM spatial_ref_sys WHERE proj4text LIKE '%NAD83%'
将srtext的文本与PRISM光栅的元数据空间属性进行比较，我们发现光栅位于EPSG（SRID 4269）中。
可以使用raster2pgsql将棱镜光栅加载到chp05.PRISM表中，该表将以与shp2pgsql命令类似的方式将光栅文件导入数据库：
> raster2pgsql -s 4269 -t 100x100 -F -I -C -Y .PRISM_tmin_provisional_4kmM2_*_asc.asc chp05.prism | psql -d postgis_sample -U postgres
使用以下标志调用raster2pgsql命令：
-s:此标志将SRID4269指定给导入的光栅。
-t:这个标志表示瓷砖的大小。它将导入的光栅分割成更小、更易于管理的部分；添加到表中的每条记录最多为100 x 100像素。
-F:此标志向表中添加一列，并用光栅文件名填充该列。
-I：这个标志在表的光栅列上创建一个GIST空间索引。
-C:此标志应用于表上的标准约束集。标准约束集包括对尺寸、比例、倾斜和左上角的检查坐标和SRID。
-Y：这个标志指示raster2pgsql使用COPY语句而不是INSERT语句。复制通常比插入快。
我们将-F传递给raster2pgsql是有原因的。如果你看棱镜光栅的文件名，你会注意到年份和月份。因此，让我们将filename列中的值转换为表中的日期：
ALTER TABLE chp05.prism ADD COLUMN month_year DATE; UPDATE chp05.prism SET month_year = ( SUBSTRING(split_part(filename, '_', 5), 0, 5) || '-' || SUBSTRING(split_part(filename, '_', 5), 5, 4) || '-01' ) :: DATE;
这是所有需要做的PRISM光栅现在。
现在，让我们导入一个航天飞机雷达地形任务（SRTM）光栅。SRTM光栅来自于美国宇航局喷气推进实验室于2000年2月进行的SRTM。该光栅和其他类似光栅可在以下网址获得：http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/.
将当前目录更改为SRTM目录：
>>cd SRTM
确保使用gdalinfo抽查SRTM光栅，以确保其有效，并具有坐标系的值。选中后，将SRTM光栅导入chp05.SRTM表：
我们还需要导入一个shapefile文件，是City和County of 旧金山提供的旧金山，可以在本书的文件夹下找到，也可以在以下链接下载：https://data.sfgov.org/Geographic-Locations-and-Boundaries/SF-Shoreline-and-Islands/rgcx-5tix。
本书文件中的旧金山边界将用于许多后续配方，必须按如下方式加载到数据库中：
5.1.3 How it works...
在这个配方中，我们导入了其余配方所需的棱镜和SRTM光栅。我们还导入了一个包含旧金山边界的shapefile，用于各种光栅分析。现在，开始玩吧！
5.2 Working with basic raster information and analysis
到目前为止，我们已经检查并将PRISM和SRTM光栅导入postgis_cookbook数据库的chp05 schema。我们现在将继续处理数据库中的光栅。
5.2.1 Getting Ready
在这个食谱中，我们探索了一些功能，这些功能可以提供对postgis_CookBook数据库中的光栅属性和特性的洞察。这样做，我们可以看到数据库中找到的内容是否与访问gdalinfo提供的信息匹配。
5.2.2 How to do it...
PostGIS包含“光栅列”视图，以提供数据库中所有光栅列的高级摘要。此视图在功能和形式上类似于“几何体”视图和“地理”视图。
让我们在raster_columns视图中运行以下SQL查询，以查看prism表中有哪些可用信息：
SELECT 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, ST_AsText(extent) AS extent FROM raster_columns WHERE r_table_name = 'prism';
SQL查询将返回类似以下内容的记录：
如果回顾其中一个棱镜光栅的gdalinfo输出，您将看到比例（像素大小）的值匹配。传递给raster2pgsql的标志，指定了tile大小和SRID，起作用了。
让我们看看单个光栅图块的元数据是什么样子的。我们将使用ST_Metadata（）函数：
SELECT rid, (ST_Metadata(rast)).* FROM chp05.prism WHERE month_year = '2017-03-01'::date LIMIT 1;
输出将类似于以下内容：
使用ST_BandMetadata（）检查记录ID 54处的第一个也是唯一一个光栅分幅带：
SELECT rid, (ST_BandMetadata(rast, 1)).* FROM chp05.prism WHERE rid = 54;
结果表明，波段为32BF像素类型，NODATA值为-9999。NODATA值是分配给空像素的值：