PostGIS Cookbook 读书笔记

关键点

• 将PostGIS与网络融合，并引入OGC标准，例如WMS和WFS。

• 将2D和3D的矢量数据、栅格数据和轨迹数据转化为可用的表。

• 将从PostGIS数据库的数据可视化，使用桌面GIS程序，例如QGIS和OpenJUMP。

• 容易使用的方法，高级分析功能和实际的应用。

你将会学到

• 使用现有的工具从PostGIS数据库，导入和导出地理数据

• 利用PostgreSQL和PostGIS结合提供的功能构造空间数据

• 使用一组PostGIS功能执行基本和高级矢量分析

• 使用Python连接PostGIS

• 学习使用PostGIS周围的编程框架

• 维护、优化和微调空间数据以实现长期生存能力

• 探索PostGIS的3D功能，包括LiDAR点云和基于运动结构（SfM）技术的点云

• 使用X3D标准通过Web分发3D模型

• 使用PostGIS使用开放的地理空间联盟web标准开发强大的GIS web应用程序

• 管理PostGIS栅格

前言

本实用书将帮助你充分掌握基本的和高级的矢量、栅格和路径方法。你将学会使用数据的保持、优化和性能，它们将会帮助你融合PostGIS到大的桌面和网络工具生态系统中。

有这本参考书，你将会掌握所有的工具和指南从而有效地管理你的数据库系统，并且随着你的项目需求的变化而做出更好的决策。

作者简介

  Paolo Corti

是一名环境工程师，在GIS领域拥有20年的经验，目前在哈佛大学地理分析中心担任地理空间工程师研究员。他是开源地理空间技术和Python的倡导者，OSGeo宪章成员，pycsw和GeoNode项目指导委员会成员。他是本书第一版的合著者，也是Packt的《掌握QGIS》一、二版的审稿人。

https://subscription.packtpub.com/book/big_data_and_business_intelligence/9781849518666

https://www.doc88.com/p-2874898517310.html?r=1

https://qcsdn.com/q/a/221718.html

这本书的内容：

第1章，将数据移入和移出PostGIS，介绍了将空间和非空间数据导入和导出PostGIS的过程。这些过程包括使用PostGIS和第三方（如GDAL/OGR）提供的实用程序。

第2章，有用的结构，讨论了如何利用PostgreSQL提供的机制来组织PostGIS数据。这些机制用于规范可能不干净和非结构化的导入数据。

第3章，使用矢量数据-基础，介绍了PostGIS中通常对矢量进行的操作，即PostGIS中的几何和地理。操作包括无效几何图形的处理、确定几何图形之间的关系以及简化复杂几何图形。

第4章，使用矢量数据-高级配方，深入分析几何图形的高级主题。您将学习如何使用KNN过滤器来提高邻近查询的性能，从LiDAR数据创建多边形，以及计算可用于邻域分析的Voronoi单元。

第5章，使用光栅数据，介绍了在PostGIS中操作光栅的实际工作流程。您将学习如何导入光栅、修改光栅、对光栅进行分析以及以标准光栅格式导出光栅。

第6章，介绍了pgRouting扩展，它为PostGIS提供了图形遍历和分析功能。本章中的方法回答了现实世界中的问题，即有条件地从A点导航到B点，并精确地模拟复杂的路线，如水道。

第7章，第n维，重点介绍了PostGIS中处理和分析多维空间数据的工具和技术，包括LiDAR源点云。所涵盖的主题包括将点云加载到PostGIS、从点云创建2.5D和3D几何图形以及若干摄影测量原理的应用。

第8章，PostGIS编程，展示了如何使用Python语言编写在PostGIS上操作并与之交互的应用程序。编写的应用程序包括在PostGIS中读写外部数据集的方法，以及使用OpenStreetMap数据集的基本地理编码引擎。

第9章，PostGIS和Web，介绍了如何使用OGC和rest web服务将PostGIS数据和服务传递到Web。本章讨论如何使用MapServer和GeoServer提供OGC、WFS和WMS服务，以及如何从客户端（如OpenLayers和sliple）使用它们。然后介绍如何使用GeoDjango构建web应用程序，以及如何在Mapbox应用程序中包含PostGIS数据。

第10章，维护、优化和性能调优，从PostGIS中后退了一步，重点介绍PostgreSQL数据库服务器的功能。通过利用PostgreSQL提供的工具，您可以确保空间和非空间数据的长期生存能力，并最大限度地提高各种PostGIS操作的性能。此外，它还探索了PostgreSQL中的地理空间切分和并行性等新特性。

第11章，使用桌面客户端，告诉您如何使用各种开源桌面GIS应用程序使用和操作PostGIS中的空间数据。本文讨论了几种应用程序，以突出显示与空间数据交互的不同方法，并帮助您找到适合该任务的工具。

第12章，位置隐私保护机制简介，介绍了位置隐私的概念，并介绍了两种不同的位置隐私保护机制的实现，这些机制可以包含在商业应用程序中，为用户的位置数据提供基本级别的保护。

如何充分利用这本书
在深入阅读本书之前，您需要安装最新版本的PostgreSQL和PostGIS（分别为9.6或103和2.3或2.41）。如果您喜欢图形化SQL工具，您可能还需要安装pgAdmin（1.18）。对于大多数计算环境（Windows、Linux、macOS X），安装程序和软件包都包含PostGIS的所有必需依赖项。PostGIS所需的最小依赖项是PROJ.4、GEOS、libjson和GDAL。要理解和修改本书中的代码，需要对SQL语言有基本的了解。

下载示例代码文件：https://github.com/PacktPublishing/PostGIS-Cookbook-Second-Edition

使用的约定：这本书中使用了许多文本约定。

CodeInText: Indicates code words in text, database table names, folder names, filenames, file extensions, pathnames, dummy URLs, user input, and Twitter handles. Here is an example: "We will import the firenews.csv file that stores a series of web news collected from various RSS feeds."

文本中的代码：表示文本中的代码，数据库表名，文件夹名，文件名，文件扩展名，路径名，虚拟URL，用户输入，和推特句柄。这里是一个例子：“我们将导入firenews.csv文件，该文件中存储着一系列的网络新闻，从多个RSS订阅中收集的。”

代码块设置如下：

SELECT ROUND(SUM(chp02.proportional_sum(ST_Transform(a.geom,3734), b.geom, b.pop))) AS population
   FROM nc_walkzone AS a, census_viewpolygon as b
  WHERE ST_Intersects(ST_Transform(a.geom, 3734), b.geom)
  GROUP BY a.id;

当我们希望提请您注意代码块的特定部分时，相关行或项以粗体显示：

SELECT ROUND(SUM(chp02.proportional_sum(ST_Transform(a.geom,3734), b.geom, b.pop))) AS population
   FROM nc_walkzone AS a, census_viewpolygon as b
  WHERE ST_Intersects(ST_Transform(a.geom, 3734), b.geom)
  GROUP BY a.id;

任何命令行输入或输出的编写方式如下：

> raster2pgsql -s 4322 -t 100x100 -F -I -C -Y C:postgis_cookbookdatachap5PRISMus_tmin_2012.*.asc chap5.prism | psql -d postgis_cookbook

Moving Data In and Out of PostGIS

本章，我们将介绍：

• 使用PostGIS函数导入非空间表格数据（CSV）
• 使用GDAL导入非空间表格数据（CSV）
• 使用shp2pgsql导入shapefile
• 使用ogr2ogr GDAL命令导入和导出数据
• 处理数据集的批量导入和导出
• 使用pgsql2shp PostGIS命令将数据导出到shapefile
• 使用osm2pgsql命令导入OpenStreetMap数据
• 使用raster2pgsql PostGIS命令导入光栅数据
• 一次导入多个光栅
• 使用gdal_translate和gdalwarp gdal命令导出光栅

介绍
PostGIS是PostgreSQL数据库的开源扩展，支持地理对象；在这本书中，你会发现食谱，将指导你一步一步地探索它提供的不同功能。
这本书的目的是成为一个有用的工具，了解PostGIS的能力，以及如何应用它们在任何时候。每一个菜谱都有一个准备阶段，用以组织你的工作空间，然后是你为实现任务的主要目标而需要执行的一系列步骤，其中包括你将需要的所有外部命令和SQL语句（这些命令和语句已经在Linux、Mac和Windows环境中进行了测试），最后是配方的小结。本书将介绍地理信息系统和基于位置的服务中的大量常见任务，这使它成为您的技术图书馆中的必备书。
在第一章中，我们将向您展示一组配方，涵盖从PostGIS空间数据库导入和导出地理数据的不同工具和方法，因为在GIS中执行的几乎所有常见操作都是从插入或导出地理空间数据开始的。

Importing nonspatial tabular data (CSV) using PostGISfunctions 使用postgis函数导入非空间表格数据（CSV）

有两种方法可以导入逗号分隔值（CSV）文件，该文件在PostGIS中存储属性和几何图形。在这个方法中，我们将使用PostgreSQL COPY命令和一系列PostGIS函数来导入这样一个文件。

准备
我们将导入firenews.csv文件，该文件存储从欧洲森林火灾信息系统（EFFIS）上下文中与欧洲森林火灾相关的各种RSS提要收集的一系列web新闻，可在http://effis.jrc.ec.europa.eu/.
对于每个新闻提要，都有诸如地名、火灾大小（以公顷为单位）、URL等属性。最重要的是，有x和y字段以十进制表示地理本地化新闻的位置（在WGS 84空间引用系统中，SRID=4326）。
对于Windows机器，有必要安装OSGeo4W，这是一组开放源码的地理库，允许对数据集进行操作。链接是：https://trac.osgeo.org/osgeo4w/
此外，在Path环境变量中包含OSGeo4W和Postgres二进制文件夹，以便能够从PC中的任何位置执行命令。

怎么做。。。
完成此配方需要遵循的步骤如下所示：
检查CSV文件firenews.CSV的结构，您可以在图书数据集中找到该文件（如果您在Windows上，请使用诸如Notepad之类的编辑器打开CSV文件）。

$ cd ~/postgis_cookbook/data/chp01/      
$ head -n 5 firenews.csv

上述命令的输出如下所示：

我们使用psql客户端连接到PostgreSQL，但是您可以使用您最喜欢的客户端，例如pgAdmin。

使用psql客户端，我们将不显示主机和端口选项，因为我们将假定您在标准端口上使用本地PostgreSQL安装。如果不是这样，请提供这些选项！

使用Copy命令将CSV文件中的记录复制到PostgreSQL表中（如果您在Windows上，请使用c: emp而不是/tmp这样的输入目录），如下所示：

postgis_cookbook=> COPY chp01.firenews (
         x, y, place, size, update, startdate,
         enddate, title, url
      ) FROM '/tmp/firenews.csv' WITH CSV HEADER;

PostgreSQL 错误: 编码 "UTF8" 的字符 0x0xc9 0x99 在编码 "GBK" 没有相对应值：https://www.zhangbj.com/p/689.html

消息告诉您，不允许Postres读取文件：https://www.codenong.com/54031813/

确保firenews.csv文件位于PostgreSQL进程用户可以访问的位置。例如，在Linux中，将文件复制到/tmp目录中。如果您在Windows上，则很可能需要在复制之前将编码设置为UTF-8:postgis_cookbook=#set client_encoding to'UTF-8'；并记住设置完整路径“c:\tmpfirenews.csv”。

检查所有记录是否已从CSV文件导入PostgreSQL表：

postgis_cookbook=> SELECT COUNT(*) FROM chp01.firenews;

请检查与此新表相关的记录是否在PostGIS几何图形列元数据视图中：

postgis_cookbook=# SELECT f_table_name,f_geometry_column, coord_dimension, srid, type  FROM geometry_columns where f_table_name = 'firenews';

在postgis2.0之前，您必须通过两个不同的步骤创建一个包含空间数据的表；实际上，geometry_columns视图是一个需要手动更新的表。为此，必须使用AddGeometryColumn函数来创建列。例如，这是针对以下配方：

postgis_cookbook=> CREATE TABLE chp01.firenews(
x float8,
y float8,
place varchar(100),
size float8,
update date,
startdate date,
enddate date,
title varchar(255),
url varchar(255))
WITHOUT OIDS;postgis_cookbook=> SELECT AddGeometryColumn('chp01', 'firenews', 'the_geom', 4326, 'POINT', 2);chp01.firenews.the_geom SRID:4326 TYPE:POINT DIMS:2

在postgis2.0中，如果您愿意，仍然可以使用AddGeometryColumn函数；但是，您需要将其use_typmod参数设置为false。

现在，使用ST_MakePoint或ST_PointFromText函数导入几何列中的点（使用以下两个update命令之一）：

postgis_cookbook=> UPDATE chp01.firenews      
SET the_geom = ST_SetSRID(ST_MakePoint(x,y), 4326); 

postgis_cookbook=> UPDATE chp01.firenews       
SET the_geom = ST_PointFromText('POINT(' || x || ' ' || y || ')',     4326); 

检查几何体字段是如何在表中的某些记录中更新的：

postgis_cookbook=# SELECT place, ST_AsText(the_geom) AS wkt_geom FROM chp01.firenews ORDER BY place LIMIT 5; 

CMD当前代码页修改

最后，为表的几何列创建空间索引：

postgis_cookbook=> CREATE INDEX idx_firenews_geom   ON chp01.firenews USING GIST (the_geom); 

它的工作原理。。。
这个方法向您展示了如何使用COPY PostgreSQL命令在PostGIS中加载非空间表格数据（CSV格式）。
在创建表并将CSV文件行复制到PostgreSQL表之后，您使用PostGIS提供的一个几何构造函数（用于二维点的stu MakePoint和stu PointFromText）更新了几何列。
这些几何体构造函数（在本例中为STïMakePoint和STïPointFromText）必须始终提供空间参考系统标识符（SRID）和点坐标来定义点几何体。
数据库中任何表中添加的每个几何字段都会在geometryu columns PostGIS元数据视图中使用一条记录进行跟踪。在以前的PostGIS版本（<2.0）中，geometryfields视图是一个表，需要手动更新，可能需要使用方便的AddGeometryColumn函数。
出于同样的原因，在以前的PostGIS版本中，为了在删除几何图形列或删除空间表时保持更新的geometryu columns视图，有DropGeometryColumn和DropGeometryTable函数。使用postgis2.0及更新版本，您不需要再使用这些函数，但是可以使用标准的ALTER table、DROP column和DROP table SQL命令安全地删除列或表。
在配方的最后一步中，您已经在表上创建了一个空间索引以提高性能。请注意，与字母数字数据库字段一样，索引只有在使用SELECT命令读取数据时才能提高性能。在本例中，您对表进行了许多更新（INSERT、UPDATE和DELETE）；根据场景的不同，在更新之后删除和重新创建索引可能会花费更少的时间。

Importing nonspatial tabular data (CSV) using GDAL使用GDAL导入非空间表格数据（CSV）

作为上一个方法的替代方法，您将使用ogr2ogr GDAL命令和GDAL OGR虚拟格式将CSV文件导入PostGIS。地理空间数据抽象库（GDAL）是一个用于光栅地理空间数据格式的转换器库。OGR是一个相关的库，它为矢量数据格式提供了类似的功能。

这一次，作为额外的步骤，您将只导入文件中的一部分要素，并将它们重新投影到不同的空间参照系统。

准备
您将从NASA的地球观测系统数据和信息系统（EOSDIS）将Global_24h.csv文件导入PostGIS数据库。
您可以从本章书籍的数据集目录中复制该文件。
此文件表示中分辨率成像光谱仪（MODIS）卫星在过去24小时内探测到的全球活动热点。对于每一行，都有以十进制度数表示的热点坐标（纬度、经度）（在WGS 84空间参考系中，SRID=4326），以及一系列有用的字段，如采集日期、采集时间和卫星类型，仅举几例。
您将只导入标记为T（Terra MODIS）的卫星类型扫描的现役火力数据，并使用球面墨卡托投影坐标系（EPSG:3857；它有时被标记为EPSG:900913，其中900913代表1337年的谷歌，因为它最初被谷歌地图广泛使用。

怎么做。。。
完成此配方所需的步骤如下：
分析Global_24h.csv文件的结构（在Windows中，使用记事本等编辑器打开csv文件）：

$ cd ~/postgis_cookbook/data/chp01/      
$ head -n 5 Global_24h.csv

创建一个GDAL虚拟数据源，它只包含一个从Global_24h.csv文件派生的层。为此，请在CSV文件所在的同一目录中创建一个名为global_24h.vrt的文本文件，并按如下方式进行编辑：

<OGRVRTDataSource>
    <OGRVRTLayer name="Global_24h">
        <SrcDataSource>Global_24h.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:4326</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="longitude" y="latitude"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

使用ogrinfo命令，检查GDAL是否正确识别虚拟层。例如，分析层的模式及其第一个特性（fid=1）：

您还可以尝试使用支持GDAL/OGR虚拟驱动程序（如quantumgis（QGIS））的桌面GIS打开虚拟层。在下面的屏幕截图中，将显示Global_24h层以及您可以在本书的数据集目录中找到的国家/地区的shapefile：

现在，使用ogr2ogr GDAL/OGR命令将虚拟层导出为PostGIS中的新表（为了使此命令可用，需要将GDAL安装文件夹添加到OS的PATH变量中）。您需要使用-f选项指定输出格式，-tu srs选项将点投影到EPSG:3857空间引用，-where选项仅加载MODIS Terra卫星类型的记录，-lco layer creation选项提供要存储表的模式：

$ ogr2ogr -f PostgreSQL -t_srs EPSG:3857      PG:"dbname='postgis_cookbook' user='me' password='mypassword'"      -lco SCHEMA=chp01 global_24h.vrt -where "satellite='T'"       -lco GEOMETRY_NAME=the_geom

检查ogr2ogr命令如何创建表，如以下命令所示：

$ pg_dump -t chp01.global_24h --schema-only -U me postgis_cookbook 

现在，检查应该出现在“geometry_columns”元数据视图中的记录：

检查表中导入了多少条记录：

postgis_cookbook=# SELECT count(*) FROM chp01.global_24h; 

注意坐标是如何从EPSG:4326投影到EPSG:3857的：

postgis_cookbook=# SELECT ST_AsEWKT(the_geom)      FROM chp01.global_24h LIMIT 1;

上述命令的输出如下：

 

它是如何工作的。。。
如GDAL文件所述：
“OGR虚拟格式是一个驱动程序，它根据XML控制文件中指定的条件转换从其他驱动程序读取的功能。”
GDAL支持读取和写入存储为CSV文件的非空间表格数据，但我们需要使用虚拟格式从CSV文件中的属性列（每个点的经纬度坐标）派生层的几何图形。为此，您至少需要在驱动程序中指定CSV文件的路径（SrcDataSource元素）、几何体类型（GeometryType元素）、图层的空间参照定义（LayerRS元素）以及驱动程序获取几何信息的方式（GeometryField元素）。
使用OGR虚拟格式还有许多其他选项和理由；如果您有兴趣更好地了解，请参阅以下网址提供的GDAL文档http://www.gdal.org/drv_vrt.html.
在正确创建虚拟格式后，原平面非空间数据集由GDAL和基于GDAL的软件支持。这就是为什么我们可以使用GDAL命令（如ogrinfo和ogr2ogr）以及QGIS等桌面GIS软件来操作这些文件的原因。
一旦我们验证GDAL能够正确地从虚拟驱动程序读取这些特性，我们就可以使用流行的ogr2ogr命令行工具轻松地将它们导入PostGIS中。ogr2ogr命令有很多选项，请参阅其文档http://www.gdal.org/ogr2ogr.htmlfor 更深入的讨论。
在这个食谱中，您刚刚看到了其中的一些选项，例如：

  -其中：用于仅导出原始要素类的选择
  -srs：用于将数据重新投影到不同的空间参考系统
  -lco层创建：它用于提供存储表的模式（没有它，新的空间表将在公共模式中创建）和输出层中几何字段的名称

Importing shapefiles with shp2pgsql使用shp2pgsql导入shapefile

如果需要在PostGIS中导入shapefile，则至少有两个选项，如前面看到的ogr2ogr GDAL命令或shp2pgsql PostGIS命令。
在这个方法中，您将使用shp2pgsql命令在数据库中加载一个shapefile，使用ogrinfo命令分析它，并在QGIS桌面软件中显示它。

怎么做。。。
完成此配方需要遵循的步骤如下：
使用OGR2OGR命令从上一个配方中创建的虚拟驱动程序创建shapefile（注意，在这种情况下，您不需要指定-f选项，因为shapefile是ogr2ogr命令的默认输出格式）：

$ ogr2ogr global_24h.shp global_24h.vrt

使用shp2pgsql命令为shapefile生成SQL转储文件。您将使用-G选项生成使用地理类型的PostGIS空间表，使用-I选项生成几何列上的空间索引：

$ shp2pgsql -G -I global_24h.shp        chp01.global_24h_geographic > global_24h.sql

分析Global_24h.sql文件（在Windows中，使用文本编辑器，如记事本）：

在PostgreSQL中运行global_24h.sql文件：

$ psql -U me -d postgis_cookbook -f global_24h.sql

如果您使用的是Linux，则可以通过以下方式将最后两个步骤中的命令串联在一行中：

$ shp2pgsql -G -I global_24h.shp chp01.global_24h_geographic | psql -U me -d postgis_cookbook

检查元数据记录是否在geography_columns视图中可见（而不是在geometry_columns视图中，与shp2pgsql命令的-G选项一样，我们选择了geography类型）：

postgis_cookbook=# SELECT f_geography_column, coord_dimension, srid, type FROM geography_columns    WHERE f_table_name = 'global_24h_geographic';

上述命令的输出如下：

使用ogrinfo分析新的PostGIS表（使用-fid选项仅显示表中的一条记录）：

$ ogrinfo PG:"dbname='postgis_cookbook' user='me'  password='mypassword'" chp01.global_24h_geographic -fid 1

上述命令的输出如下：

现在，打开QGIS并尝试将新图层添加到地图中。导航到Layer | Add Layer | Add PostGIS layers并提供连接信息，然后将图层添加到地图中，如以下屏幕截图所示：

它的工作原理。。。
PostGIS命令shp2pgsql允许用户在PostGIS数据库中导入shapefile。基本上，它生成一个PostgreSQL转储文件，可以通过在PostgreSQL中运行它来加载数据。
SQL文件通常由以下部分组成：
CREATE TABLE部分（如果没有选择-a选项，在这种情况下，表应该已经存在于数据库中）
INSERT INTO部分（对于要从shapefile导入的每个特征，都有一条INSERT语句）
创建索引部分（如果选择了-I选项）
与ogr2ogr不同的是，无法为要导入的shapefile中的特征进行空间或属性选择（-spat，-where ogr2ogr options）。另一方面，通过shp2pgsql命令，也可以导入特征的mcoordinate（ogr2ogr在编写时仅支持x、y和zat）。
要获得shp2pgsql命令选项及其含义的完整列表，只需在shell中键入命令名（如果在Windows上，则在命令提示符中键入），然后检查输出。

还有更多。。。
PostGIS中有GUI工具来管理数据进出，通常集成到QGIS等GIS桌面软件中。在这本书的最后一章中，我们将看一看最受欢迎的一章。

Importing and exporting data with the ogr2ogr GDAL command使用ogr2ogr GDAL命令导入和导出数据

怎么做。。。
完成此配方所需的步骤如下：
将wborders.zip存档解压缩到您的工作目录。你可以在书的数据集中找到这个。
使用ogr2ogr命令在PostGIS中导入世界各国的形状文件（wborders.shp）。使用ogr2ogr中的一些选项，您将仅导入subsection=2（非洲）的要素以及ISO2和NAME属性，并将featureclass重命名为Africau countries：
$ogr2ogr-f PostgreSQL-sql“SELECT ISO2，NAME AS countryu NAME FROM wborders WHERE REGION=2”-nlt MULTIPOLYGON PG:“dbname='postgisu cookbook'user='me'password='mypassword'”-nln africau countries-lco SCHEMA=chp01-lco GEOMETRYu NAME=the u geom wborders.shp
检查是否在PostGIS中正确导入shapefile，查询数据库中的空间表或在桌面GIS中显示它。
查询PostGIS，从上一个配方中创建的globalu 24h表中获取具有最高亮度温度的100个活动热点（brightu t31字段）的列表：
postgis_cookbook=#选择st_AsText（the_geom）作为_geom，bright_t31 FROM chp01.global_24h ORDER BY bright_t31 DESC LIMIT 100；
上述命令的输出如下：

Handling batch importing and exporting of datasets处理数据集的批量导入和导出

在许多GIS工作流中，有一个典型的场景，PostGIS表的子集必须以文件系统格式（最典型的是shapefile或spatialite数据库）部署给外部用户。通常，也有相反的过程，从不同用户收到的数据集必须上传到PostGIS数据库。

在这个配方中，我们将模拟这两个数据流。您将首先创建用于处理PostGIS中的shapefile的数据流，然后创建用于上载shapefile的反向数据流。

您将使用bash脚本和ogr2ogr命令来完成它。

Exporting data to a shapefile with the pgsql2shp PostGIS command使用pgsql2shp PostGIS命令将数据导出到shapefile

在这个方法中，您将使用PostGIS发行版附带的pgsql2shp命令将PostGIS表导出到shapefile。

Importing OpenStreetMap data with the osm2pgsql command使用osm2pgsql命令导入OpenStreetMap数据

在这个方法中，您将使用osm2pgsql命令将OpenStreetMap（OSM）数据导入PostGIS。

您将首先从OSM网站下载一个示例数据集，然后使用osm2pgsql命令导入它。

您将在GIS桌面软件中添加导入的图层，并生成一个视图以获取子数据集，使用hstore PostgreSQL附加模块根据其标记提取特征。

Importing raster data with the raster2pgsql PostGIS command使用raster2pgsql PostGIS命令导入光栅数据

PostGIS 2.0现在完全支持光栅数据集，并且可以使用raster2pgsql命令导入光栅数据集。
在此配方中，您将使用raster2pgsql命令将光栅文件导入PostGIS。此命令包含在从版本2.0开始的任何PostGIS发行版中，能够为任何GDAL光栅支持的格式生成要加载到PostGIS中的SQL转储（与shp2pgsql命令对ShapeFile所做的方式相同）。
将光栅加载到PostGIS后，您将使用SQL命令（分析数据库中包含的光栅元数据信息）和gdalinfo命令行实用程序（了解输入raster2pgsql参数在PostGIS导入过程中的反映方式）对其进行检查。
最后，您将在桌面GIS中打开光栅，并尝试一个基本的空间查询，混合矢量表和光栅表。

Importing multiple rasters at a time一次导入多个光栅

此配方将指导您一次导入多个光栅。

首先使用raster2pgsql命令将一些不同的单波段光栅导入到唯一的单波段光栅表中。

然后，您将尝试另一种方法，将虚拟光栅中的原始单波段光栅与每个原始光栅的一个波段合并，然后将多波段光栅加载到光栅表中。为此，您将使用GDAL gdalbuildvrt命令，然后使用raster2pgsql将数据加载到PostGIS。

准备
确保所有原始光栅数据集都已用于上一个配方。

它是如何工作的。。。
可以使用raster2pgsql命令在PostGIS中导入光栅数据集。
GDAL PostGIS光栅迄今不支持写入操作；因此，目前，您不能使用GDAL命令，例如GDAL_utranslate和gdalwarp，这将在不久的将来发生变化，因此，当您阅读本章时，您可能会有这样一个额外的选项。
在一个场景中，您有多个光栅在不同的时间代表相同的变量，如本配方中所述，将所有原始的光栅存储为PostGIS中的单个表是有意义的。在这个配方中，我们有相同的变量（平均最高温度），每个月用一个光栅表示。您已经看到，您可以以两种不同的方式进行：
将每个单个光栅（表示不同月份）附加到同一PostGIS单波段光栅表中，并从文件名列中的值（使用-F raster2pgsql选项添加到表中）中导出与月份相关的信息。
使用gdalbuildvrt生成多波段光栅（一个带12个波段的光栅，每个月一个），并使用raster2pgsql命令将其导入单个多波段PostGIS表中。

Exporting rasters with the gdal_translate and gdalwarp GDAL commands使用gdal_translate和gdalwarp gdal命令导出光栅

在本教程中，你将看到几个将PostGIS栅格转为不同的栅格格式的主要的选项。它们都是以命令行工具的形式提供，gdal_translate和gdalwarp，由GDAL提供。

Structures That Work

本章，我们将介绍：

• Using geospatial views使用地理空间视图
• Using triggers to populate the geometry column使用触发器填充几何图形列
• Structuring spatial data with table inheritance利用表继承构造空间数据
• Extending inheritance – table partitioning扩展继承-表分区
• Normalizing imports输入正常化
• Normalizing internal overlays规范化内部覆盖
• Using polygon overlays for proportional census estimates使用多边形覆盖进行比例普查估计

介绍
本章重点介绍如何使用PostgreSQL和PostGIS组合提供的功能来构造数据。这些方法对于结构化和清理导入的数据、将表格数据转换为空间数据以及使用PostgreSQL和PostGIS强大组合特有的功能维护表和数据集之间的关系都是有用的。我们将使用三类技术来利用这些功能：使用视图和触发器自动填充和修改数据，使用PostgreSQL表继承面向对象，以及使用PostGIS函数（存储过程）来重建和规范有问题的数据。
自动填充数据是本章的开始。通过利用PostgreSQL视图和触发器，我们可以创建临时和灵活的解决方案，在表之间和表内创建连接。通过扩展，对于更正式或结构化的情况，PostgreSQL提供了表继承和表分区，这允许表之间的显式层次关系。这在对象继承模型强制数据关系的情况下非常有用，这些关系要么更好地表示数据，从而提高效率，要么随着时间的推移减少维护和访问数据集的管理开销。通过PostGIS扩展该功能，继承不仅可以应用于常用的表属性，还可以应用于利用表之间的空间关系，从而提高对非常大的数据集的查询效率。最后，我们将探讨提供数据输入表规范化的postgissql模式，以便将来自平面文件系统或未规范化的数据集转换为我们在数据库中期望的形式。

Using geospatial views使用地理空间视图

PostgreSQL中的视图允许在交互表单中临时表示数据和数据关系。在这个方法中，我们将使用视图来允许基于表格输入自动创建点数据。我们可以想象这样一种情况：输入的数据流是非空间的，但包括经纬度或其他坐标。我们希望自动将这些数据显示为空间中的点。

准备
我们可以很容易地创建一个视图作为空间数据的表示。创建视图的语法类似于创建表，例如：

CREATE VIEW viewname AS
   SELECT...

在前面的命令行中，SELECT查询为我们处理数据。让我们从一个小数据集开始。在这种情况下，我们将从一些随机点开始，它们可能是真实的数据。
首先，我们创建将从中构造视图的表，如下所示：

-- Drop the table in case it exists
DROP TABLE IF EXISTS chp02.xwhyzed CASCADE;
CREATE TABLE chp02.xwhyzed
-- This table will contain numeric x, y, and z values 
(
  x numeric,
   y numeric,
   z numeric
)
WITH (OIDS=FALSE); 
ALTER TABLE chp02.xwhyzed OWNER TO me; 
-- We will be disciplined and ensure we have a primary key 
ALTER TABLE chp02.xwhyzed ADD COLUMN gid serial; 
ALTER TABLE chp02.xwhyzed ADD PRIMARY KEY (gid);

现在，让我们使用以下查询将测试数据填充到这里：

INSERT INTO chp02.xwhyzed (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 

怎么做。。。
现在，要创建视图，我们将使用以下查询：

-- Ensure we don't try to duplicate the view
DROP VIEW IF EXISTS chp02.xbecausezed;
-- Retain original attributes, but also create a point   attribute from x and y
CREATE VIEW chp02.xbecausezed AS
SELECT x, y, z, ST_MakePoint(x,y)
FROM chp02.xwhyzed;

它的工作原理。。。
我们的视图是使用PostGIS的st_MakePoint函数对现有数据进行简单的转换。st_MakePoint函数接受两个数字的输入来创建一个PostGIS点，在本例中，我们的视图只是使用x和y值来填充数据。每当对表进行更新以添加具有X和Y值的新记录时，视图将填充一个点，这对于不断更新的数据非常有用。
这种方法有两个缺点。首先，我们没有在视图中声明我们的空间参考系，因此任何使用这些点的软件都不知道我们使用的坐标系，即，它是地理坐标系（纬度/经度）还是平面坐标系。我们将很快解决这个问题。第二个问题是，许多访问这些点的软件系统可能无法自动检测和使用表中的空间信息。在使用触发器填充几何列配方中解决了此问题。
空间参考系统标识符（SRID）允许我们为给定的数据集指定坐标系。编号系统是一个简单的整数值，用于指定给定的坐标系。SRID最初来源于欧洲石油调查集团（EPSG），现在由国际石油和天然气生产商协会（OGP）的调查和定位委员会维护。srid的有用工具是空间引用(http://spatialreference.org)和Prj2EPSG(http://prj2epsg.org/search).

还有更多。。。
为了解决如何工作…一节中提到的第一个问题（缺少空间参考系统），我们可以简单地将现有的ST_MakePoint函数包装到另一个指定SRID为ST_SetSRID的函数中，如下查询所示：

-- Ensure we don't try to duplicate the view 
DROP VIEW IF EXISTS chp02.xbecausezed; 
-- Retain original attributes, but also create a point   attribute from x and y 
CREATE VIEW chp02.xbecausezed AS
   SELECT x, y, z, ST_SetSRID(ST_MakePoint(x,y), 3734) -- Add ST_SetSRID
   FROM chp02.xwhyzed;

Using triggers to populate the geometry column使用触发器填充几何图形列

在这个配方中，我们设想我们的数据库中有越来越多的数据，这需要空间表示；但是，在这种情况下，我们希望每次在数据库上插入时都更新硬编码的几何体列，将X和Y值在插入数据库时转换为几何体。
这种方法的优点是，几何图形随后在“几何图形”列视图中注册，因此这种方法可以可靠地处理更多的PostGIS客户端类型，而不是创建新的地理空间视图。这还提供了允许空间索引的优势，空间索引可以显著加快各种查询的速度。

准备
我们将首先创建另一个具有x、y和z值的随机点表，如以下查询所示：

DROP TABLE IF EXISTS chp02.xwhyzed1 CASCADE; 
CREATE TABLE chp02.xwhyzed1 
(
   x numeric,
   y numeric,
   z numeric 
) 
WITH (OIDS=FALSE); 
ALTER TABLE chp02.xwhyzed1 OWNER TO me; 
ALTER TABLE chp02.xwhyzed1 ADD COLUMN gid serial; 
ALTER TABLE chp02.xwhyzed1 ADD PRIMARY KEY (gid);  
INSERT INTO chp02.xwhyzed1 (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed1 (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed1 (x, y, z)
   VALUES (random()*5, random()*7, random()*106); 
INSERT INTO chp02.xwhyzed1 (x, y, z)
   VALUES (random()*5, random()*7, random()*106);

怎么做。。。
现在我们需要一个几何列来填充。默认情况下，几何体列将填充空值。我们使用以下查询填充几何体列：

SELECT AddGeometryColumn ('chp02','xwhyzed1','geom',3734,'POINT',2);

我们现在有一个名为geom的列，SRID为3734；即，二维的点几何图形类型。由于我们有x、y和z数据，原则上，我们可以使用类似的方法填充三维点表。
由于所有几何体值当前都为空，因此我们将使用UPDATE语句填充它们，如下所示：

UPDATE chp02.xwhyzed1
   SET geom = ST_SetSRID(ST_MakePoint(x,y), 3734); 

这里的查询很简单。我们更新xwhyzed1表，并使用ST_MakePoint设置_geom列，使用x和y列构造我们的点，并将其包装在ST_SetSRID函数中，以便应用适当的空间引用信息。到目前为止，我们刚刚摆好了桌子。现在，我们需要创建一个触发器，以便在表被使用后继续填充此信息。触发器的第一部分是使用以下查询的新填充几何体函数：

CREATE OR REPLACE FUNCTION chp02.before_insertXYZ()
   RETURNS trigger AS
$$
BEGIN 
if NEW.geom is null then
    NEW.geom = ST_SetSRID(ST_MakePoint(NEW.x,NEW.y), 3734); 
end if; 
RETURN NEW; 
END;
$$ 
LANGUAGE 'plpgsql'; 

本质上，我们创建了一个函数，它完全执行我们手动所做的操作：使用ST_SetSRID和ST_MakePoint的组合来更新表的geometry列，但只更新插入的新寄存器，而不是所有表。

还有更多。。。
虽然我们已经创建了一个函数，但还没有将它作为触发器应用到表中。让我们这样做如下：

CREATE TRIGGER popgeom_insert
  BEFORE INSERT ON chp02.xwhyzed1
  FOR EACH ROW EXECUTE PROCEDURE chp02.before_insertXYZ();

假设还没有进行常规几何列更新，那么原来的五个寄存器的几何列仍然为空。现在，一旦触发器被激活，表中的任何插入都应该填充新的几何记录。让我们使用以下查询执行测试插入：

INSERT INTO chp02.xwhyzed1 (x, y, z)
  VALUES (random()*5, random()*7, 106),
  (random()*5, random()*7, 107),
  (random()*5, random()*7, 108),
  (random()*5, random()*7, 109),
  (random()*5, random()*7, 110);

检查行以验证geom列是否已使用以下命令更新：

SELECT * FROM chp02.xwhyzed1;

或使用pgAdmin：

应用常规更新后，所有寄存器的geom列上都会有一个值：

进一步扩展。。。
到目前为止，我们已经实现了一个insert触发器。如果某一行的值发生变化怎么办？在这种情况下，我们将需要一个单独的更新触发器。我们将更改原始函数以测试UPDATE情况，并使用触发器中的WHEN将更新约束到要更改的列。
另外，请注意，编写以下函数时假设用户希望始终基于更改的值更新更改的几何图形：

CREATE OR REPLACE FUNCTION chp02.before_insertXYZ()
   RETURNS trigger AS
$$ 
BEGIN 
if (TG_OP='INSERT') then
   if (NEW.geom is null) then
     NEW.geom = ST_SetSRID(ST_MakePoint(NEW.x,NEW.y), 3734);
   end if;
ELSEIF (TG_OP='UPDATE') then
   NEW.geom = ST_SetSRID(ST_MakePoint(NEW.x,NEW.y), 3734);
end if;
RETURN NEW; 
END;
  
$$ 
LANGUAGE 'plpgsql'; 

CREATE TRIGGER popgeom_insert
   BEFORE INSERT ON chp02.xwhyzed1
   FOR EACH ROW EXECUTE PROCEDURE chp02.before_insertXYZ(); 


CREATE trigger popgeom_update 
  BEFORE UPDATE ON chp02.xwhyzed1 
  FOR EACH ROW 
   WHEN (OLD.X IS DISTINCT FROM NEW.X OR OLD.Y IS DISTINCT FROM      NEW.Y) 
  EXECUTE PROCEDURE chp02.before_insertXYZ(); 

另请参见
  使用地理空间视图

Structuring spatial data with table inheritance利用表继承构造空间数据

PostgreSQL数据库的一个不寻常且有用的特性是，它允许对象继承模型应用于表。这意味着我们可以在表之间建立父/子关系，并利用这些关系以有意义的方式构造数据。在我们的例子中，我们将把它应用于水文数据。这些数据可以是点、线、多边形或更复杂的结构，但它们有一个共同点：它们在物理意义上是显式链接的，并且内在地相关；它们都是关于水的。水/水文学是一个很好的自然系统，以这种方式进行建模，因为我们的建模方法在空间上可以是非常混合的，这取决于尺度、细节、数据收集过程和许多其他因素。

准备
我们将使用的数据是从工程蓝线（见以下屏幕截图）修改的水文数据，即非常详细的水文数据，用于接近1:600的比例尺。原始应用程序中的数据作为特征线辅助详细的数字地形建模。

虽然数据本身很有用，但数据被进一步处理，将线特征与面特征分离，并对面特征进行额外的多边形化，如以下屏幕截图所示：

最后，将数据分为基本水道类别，如下所示：

此外，还进行了一个为多边形特征（如流）生成中心线的过程，这些特征实际上是线性特征，如下所示：

因此，我们有三个独立但相关的数据集：

cuyahoga_hydro_polygon 多边形
cuyahoga_hydro_polyline 多线
cuyahoga_river_centerlines 中线

现在，让我们来看看表格数据的结构。从图书库中解压缩水文文件，然后转到该目录。ogrinfo实用程序可以帮助我们实现这一点，如以下命令所示：

> ogrinfo cuyahoga_hydro_polygon.shp -al -so  

输出如下：

在每个shapefile上执行此查询时，我们会看到以下所有shapefile共用的字段：

name
hyd_type
geom_type

通过理解我们的公共字段，我们可以应用继承来完整地构造数据。

怎么做。。。
既然我们知道了公共字段，创建继承模型就很容易了。首先，我们将使用以下查询创建一个父表，其中包含所有表的公共字段：

CREATE TABLE chp02.hydrology (
   gid SERIAL PRIMARY KEY,
   "name"      text,
   hyd_type    text,
   geom_type   text,
   the_geom    geometry
); 

如果您注意到了，您会注意到我们还添加了一个几何字段，因为我们所有的shapefile都隐含着这个公共性。通过继承，插入到任何子表中的每条记录也将保存在父表中，只有这些记录将被存储，而不需要为子表指定额外的字段。

要为给定表建立继承，我们只需要使用以下查询声明子表包含的其他字段：

CREATE TABLE chp02.hydrology_centerlines (
   "length"    numeric 
) INHERITS (chp02.hydrology);

CREATE TABLE chp02.hydrology_polygon (
   area    numeric,
   perimeter    numeric 
) INHERITS (chp02.hydrology);

CREATE TABLE chp02.hydrology_linestring (
   sinuosity    numeric 
) INHERITS (chp02.hydrology_centerlines); 

现在，我们可以使用以下命令加载数据：

shp2pgsql -s 3734 -a -i -I -W LATIN1 -g the_geom cuyahoga_hydro_polygon chp02.hydrology_polygon | psql -U me -d postgis_cookbook
shp2pgsql -s 3734 -a -i -I -W LATIN1 -g the_geom cuyahoga_hydro_polyline chp02.hydrology_linestring | psql -U me -d postgis_cookbook
shp2pgsql -s 3734 -a -i -I -W LATIN1 -g the_geom cuyahoga_river_centerlines chp02.hydrology_centerlines | psql -U me -d postgis_cookbook

如果我们查看父表，我们将看到所有子表中的所有记录。以下是水文领域的截图：

将其与hydrology_linestring中的可用字段进行比较，这些字段将显示感兴趣的特定字段：

它的工作原理。。。
PostgreSQL表的继承允许我们在表之间强制本质上的层次关系。在这种情况下，我们利用继承来实现相关数据集之间的通用性。现在，如果我们想从这些表中查询数据，我们可以直接从父表中进行查询，具体如下所示，这取决于我们是要混合几何体还是只需要一个目标数据集：

SELECT * FROM chp02.hydrology

从任何子表中，我们都可以使用以下查询：

SELECT * FROM chp02.hydrology_polygon 

另请参见
通过结合使用检查约束和继承，可以扩展这个概念，以利用和优化存储和查询。有关更多信息，请参阅扩展继承-表分区配方。

Extending inheritance – table partitioning扩展继承-表分区

表分区是一种特定于PostgreSQL的方法，它将继承扩展到模型表，这些表通常在可用字段中彼此没有差异，但是子表表示基于各种因素（时间、值范围、分类，或者在我们的例子中是空间关系）的数据逻辑分区。分区的优点包括由于较小的索引和有针对性的数据扫描、大容量加载和删除而提高了查询性能，从而绕过了清空的成本。因此，它可以用来将常用的数据放在更快、更昂贵的存储上，而将剩余的数据放在更慢、更便宜的存储上。结合PostGIS，我们得到了空间划分的新功能，这对于大型数据集来说是一个非常强大的功能。

准备
我们可以使用许多大型数据集的例子，它们可以从分区中获益。在本例中，我们将使用等高线数据集。等高线是表示地形数据的有用方法，因为它们已建立并因此得到普遍解释。等高线还可用于将地形数据压缩为线性表示，从而使其能够轻松地与其他数据一起显示。
问题是，等高线数据的存储可能相当昂贵。对于一个美国县来说，两英尺的等高线可能需要20到40GB的空间，而对于一个更大的区域（如一个地区或国家）存储这样的数据，从以高效的方式访问数据集的适当部分的角度来看，可能会变得非常禁止。

怎么做。。。
在这种情况下，第一步可能是准备数据。如果我们有一个名为cuy_contours_2的整体轮廓表，我们可以选择将数据剪辑成一系列矩形，作为我们的表分区；在本例中，使用以下查询chp02.contour_clip：

CREATE TABLE chp02.contour_2_cm_only AS
   SELECT contour.elevation, contour.gid, contour.div_10, contour.div_20, contour.div_50,
   contour.div_100, cc.id, ST_Intersection(contour.the_geom, cc.the_geom) AS the_geom FROM
     chp02.cuy_contours_2 AS contour, chp02.contour_clip as cc     WHERE ST_Within(contour.the_geom,cc.the_geom
       OR
       ST_Crosses(contour.the_geom,cc.the_geom); 

我们在查询中执行两个测试。我们使用的是ST_Within，它测试给定的轮廓是否完全在我们感兴趣的区域内。如果是这样的话，我们执行一个交叉；结果几何体应该只是轮廓的几何体。

ST_Cross函数检查轮廓是否穿过我们正在测试的几何体的边界。这应该捕获所有的几何体，部分地躺在我们的区域内和部分外面。这些是我们将真正相交得到的形状。

在我们的例子中，这更容易，我们不需要这个步骤。我们的轮廓形状已经是剪裁到矩形边界的单个形状文件，如以下屏幕截图所示：

由于数据已经被裁剪到分区所需的块中，所以我们可以继续创建适当的分区。
与继承非常相似，我们首先使用以下查询创建父表：

CREATE TABLE chp02.contours (
   gid serial NOT NULL,
   elevation integer,
   __gid double precision,
   the_geom geometry(MultiLineStringZM,3734),
   CONSTRAINT contours_pkey PRIMARY KEY (gid) 
) 
WITH (
   OIDS=FALSE 
);

在这里，我们再次维护我们的约束，例如主键，并指定几何体类型（MultiLineStringZM），这不是因为这些约束将传播到子表，而是因为任何访问父表的客户机软件都可以预期这些约束。

现在我们可以开始创建从父表继承的表。在此过程中，我们将创建一个检查约束，使用以下查询指定关联几何体的限制：

CREATE TABLE chp02.contour_N2260630
   (CHECK
    (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2260000, 630000, 2260000 635000, 2265000 635000,
                 2265000 630000, 2260000 630000))',3734)
    )
  )) INHERITS (chp02.contours);

我们可以完成表结构，以便用类似的CREATE表查询来划分等高线，如下所示：

CREATE TABLE chp02.contour_N2260635
   (CHECK
     (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2260000 635000, 2260000 640000,
                 2265000 640000, 2265000 635000, 2260000 635000))', 3734)
     )
  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2260640
   (CHECK
    (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2260000 640000, 2260000 645000, 2265000 645000,
                 2265000 640000, 2260000 640000))', 3734)
     )  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2265630
   (CHECK
    (ST_CoveredBy(the_geom,ST_GeomFromText
 ('POLYGON((2265000 630000, 2265000 635000, 2270000 635000,
                 2270000 630000, 2265000 630000))', 3734)
    )
  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2265635
   (CHECK
     (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2265000 635000, 2265000 640000, 2270000 640000,
                 2270000 635000, 2265000 635000))', 3734)
     )
  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2265640
   (CHECK
     (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2265000 640000, 2265000 645000, 2270000 645000,
                 2270000 640000, 2265000 640000))', 3734)
     )
  )) INHERITS (chp02.contours);
CREATE TABLE chp02.contour_N2270630
   (CHECK
     (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2270000 630000, 2270000 635000, 2275000 635000,
                  2275000 630000, 2270000 630000))', 3734)
     )
  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2270635
   (CHECK
     (ST_CoveredBy(the_geom,ST_GeomFromText
      ('POLYGON((2270000 635000, 2270000 640000, 2275000 640000,
                  2275000 635000, 2270000 635000))', 3734)
     )
  )) INHERITS (chp02.contours); 
CREATE TABLE chp02.contour_N2270640
   (CHECK
    (ST_CoveredBy(the_geom,ST_GeomFromText
           ('POLYGON((2270000 640000, 2270000 645000, 2275000 645000,
                  2275000 640000, 2270000 640000))', 3734)
     )
  )) INHERITS (chp02.contours); 

现在我们可以通过替换文件名，使用下面的命令将contours1 ZIP文件中的等高线shapefile加载到每个子表中。如果我们愿意，我们甚至可以在父表上实现触发器，将每个insert放入其正确的子表中，尽管这可能会带来性能成本：

shp2pgsql -s 3734 -a -i -I -W LATIN1 -g the_geom N2265630 chp02.contour_N2265630 | psql -U me -d postgis_cookbook

它的工作原理。。。
CHECK约束与继承结合使用，是构建表分区所需的全部内容。在本例中，我们使用边界框作为检查约束，并简单地从父表继承列。既然我们已经有了这个功能，那么对父表的查询将在使用查询之前先检查检查约束。
这也允许我们将任何使用较少的等高线表放置在更便宜和更慢的存储上，从而允许对大型数据集进行经济高效的优化。这种结构也有利于快速变化的数据，因为更新可以应用于整个区域；该区域的整个表可以有效地删除和重新填充，而不需要跨数据集遍历。

另请参见
有关表继承的更多信息，特别是与子表中交替列的使用相关联的灵活性，请参阅前面的配方，使用表继承构造空间数据。

Normalizing imports归一化输入

数据库规范化：https://www.jianshu.com/p/588f01f565d0

数据库范式：https://blog.csdn.net/dreamrealised/article/details/10474391

通常，在空间数据库中使用的数据是从其他数据源导入的。因此，它可能并不适合我们当前的应用。在这种情况下，写一个函数来帮助将数据导入一个表可能更有用。尤其是当从平文件格式如shapefiles到关系型数据库如PostgreSQL。

shapefile是存储空间数据的事实上的格式规范，可能是空间数据最常见的传递格式。一个shapefile，不管它的名字是什么，它绝不仅仅是一个文件，而是一个文件的集合。它至少由*.shp（包含几何图形）、*.shx（索引文件）和*.dbf（包含shapefile的表格信息）组成。它是一个强大的而且有用的格式，但是作为一个平面文件，它是非关系的。每个几何图形与表中的每一行以一对一的关系相关联。

有许多结构可以作为shapefile中关系存储的代理。我们将在这里探讨一个：一个单独的字段与分隔文本的多重关系。这是一种将多个关系编码到平面文件中的常见黑客行为。另一种常见的方法是创建多个字段来存储关系安排中的单个字段。

准备
我们将使用的数据集是一个trails数据集，它具有公园系统中一组trails的线性范围。数据是来自GIS世界的典型数据；作为平面形状文件，数据中没有显式的关系结构。
首先，解压缩trails.zip文件并使用命令行进入其中，然后使用以下命令加载数据：

shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom trails chp02.trails | psql -U me -d postgis_cookbook

从线性数据来看，我们对使用类型有一些分类：

我们希望保留此信息以及名称。不幸的是，label_namefield是一个混乱的字段，其中有许多相关的名称与一个符号（&）串联在一起，如下查询所示：

SELECT DISTINCT label_name FROM chp02.trails   WHERE label_name LIKE '%&%' LIMIT 10;

它将返回以下输出：

这就是我们表的归一化的开始。

怎么做。。。
我们需要做的第一件事是找到所有没有符号的字段，并将它们用作可用轨迹的唯一列表。在我们的例子中，我们可以做到这一点，因为每个轨迹都至少有一个片段是唯一命名的，并且与另一个轨迹名没有关联。这种方法不适用于所有的数据集，因此在将这种方法应用于数据集之前，请仔细理解您的数据。
要选择不带符号的字段，请使用以下查询：

SELECT DISTINCT label_name, res
   FROM chp02.trails
   WHERE label_name NOT LIKE '%&%'
   ORDER BY label_name, res;

它将返回以下输出：

 

 接下来，我们要搜索所有与这些唯一路径名匹配的记录。这将给我们一个记录列表，作为关系。执行此搜索的第一步是将百分号（%）附加到我们的唯一列表中，以便构建一个字符串，我们可以在该字符串上使用LIKE查询进行搜索：

 psql -U postgres -d postgis_sample

最后，我们将在WITH块的上下文中使用它来进行规范化，这将为我们的第一列中的每个段以及相关的标签列提供一个惟一ID表。为了更好地度量，我们将以创建表过程的形式执行此操作，如下查询所示：

CREATE TABLE chp02.trails_names AS WITH labellike AS 
(
 SELECT '%' || label_name || '%' AS label_name, label_name as   label, res FROM
   (SELECT DISTINCT label_name, res
     FROM chp02.trails
     WHERE label_name NOT LIKE '%&%'
     ORDER BY label_name, res
   ) AS label 
) 
SELECT t.gid, ll.label, ll.res
   FROM chp02.trails AS t, labellike AS ll
   WHERE t.label_name LIKE ll.label_name
   AND
   t.res = ll.res
   ORDER BY gid;

如果我们查看创建的表的第一行，即每个表的名称，那么pgAdmin的输出如下：

现在我们有了一个关系表，我们需要一个和gid相关的几何图形表。相比之下，这非常简单，如以下查询所示：

CREATE TABLE chp02.trails_geom AS
   SELECT gid, the_geom
   FROM chp02.trails; 

它的工作原理。。。
在本例中，我们生成了一个唯一的可能记录列表，并搜索相关记录，以便建立表关系。在一个表中，我们有每个空间记录的几何图形和唯一ID；在另一个表中，我们有与每个唯一id相关联的名称。现在我们可以明确地利用这些关系。
首先，我们需要建立我们的唯一ID作为主键，如下所示：

ALTER TABLE chp02.trails_geom ADD PRIMARY KEY (gid);

现在，我们可以使用以下查询将该主键用作trailsu names表中的外键：

ALTER TABLE chp02.trails_names ADD FOREIGN KEY (gid) REFERENCES chp02.trails_geom(gid); 

这一步不是严格必要的，但会强制执行以下查询的引用完整性：

SELECT geo.gid, geo.the_geom, names.label FROM
  chp02.trails_geom AS geo, chp02.trails_names AS names
  WHERE geo.gid = names.gid;

输出如下：

还有更多。。。
如果有多个字段需要规范化（唯一化？范式？），我们可以为每个字段编写CREATE TABLE查询。
值得注意的是，本配方中的方法并不局限于使用分隔字段的情况。这种方法可以为平面文件的规范化问题提供一种相对通用的解决方案。例如，如果我们有一个例子，其中有多个字段来表示关系信息，例如label1、label2、label3，或者一条记录的类似多属性名称，那么我们可以编写一个简单的查询，将它们连接在一起，然后再将这些信息输入到我们的查询中。

Normalizing internal overlays规范化内部覆盖

来自外部源的数据可能在表结构和拓扑结构中存在问题，这是地理空间数据本身特有的问题。以多边形重叠的数据为例。如果我们的数据集具有与内部覆盖重叠的多边形，那么对面积、周长和其他度量的查询可能不会产生可预测或一致的结果。

有几种方法可以解决具有内部覆盖的多边形数据集的问题。这里介绍的一般方法最初是由折射研究的Kevin Neufeld提出的。

在编写查询的过程中，我们还将生成一个将多边形转换为线字符串的解决方案。

 准备

首先，解压use_area.zip文件并使用命令行进入其中；然后，使用以下命令加载数据集：

shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom cm_usearea_polygon chp02.use_area | psql -U me -d postgis_cookbook

怎么做。。。
现在数据被加载到数据库中的表中，我们可以利用PostGIS来展开并得到多边形的联合，这样我们就可以有一个标准化的数据集。使用此方法的第一步将是将多边形转换为line字符串。然后我们可以链接这些线并将它们转换回多边形，表示所有多边形的联合。我们将执行以下任务：

  将多边形转换为线字符串
  将线串转换回多边形
  找到结果多边形的中心点
  使用结果点查询表格关系

要将多边形转换为线串，我们需要使用ST_ ExteriorRing提取所需的多边形部分，使用ST_DumpPoints将这些部分转换为点，然后将这些点连接回线，就像使用ST_MakeLine连接dotscoloring手册一样。

进一步分解，st_ExteriorRing（几何体）将只获取多边形的外部边界。但是st_ExteriorRing返回多边形，所以我们需要获取输出并从中创建一条线。最简单的方法是使用st_DumpPoints将其转换为点，然后连接这些点。默认情况下，Dump函数返回一个名为geometry_Dump的对象，该对象不仅是简单的几何体，而且是包含整数数组的几何体。仅返回几何体的最简单方法是利用对象表示法仅提取几何体的几何体部分，如下所示：

(ST_DumpPoints(geom)).geom 

使用以下查询将几何体与ST_Exterioring一起进行拼接：

SELECT (ST_DumpPoints(ST_ExteriorRing(geom))).geom 

这将为我们提供一个点列表，按照所有点的外环的顺序排列，我们要从这些点开始使用ST_MakeLine构造线，如下查询所示：

SELECT ST_MakeLine(geom) FROM (
   SELECT (ST_DumpPoints(ST_ExteriorRing(geom))).geom) AS linpoints 

由于前面的方法是我们可能希望在许多其他地方使用的过程，因此使用以下查询从中创建函数可能是谨慎的：

CREATE OR REPLACE FUNCTION chp02.polygon_to_line(geometry)
   RETURNS geometry AS
 $BODY$
    SELECT ST_MakeLine(geom) FROM (
     SELECT (ST_DumpPoints(ST_ExteriorRing((ST_Dump($1)).geom))).geom
    ) AS linpoints
 $BODY$
   LANGUAGE sql VOLATILE;
 ALTER FUNCTION chp02.polygon_to_line(geometry)
   OWNER TO me;

现在我们有了多边形到线的函数，在我们的特殊情况下，我们仍然需要强制连接重叠的线。ST_Union函数将有助于实现这一点，如下查询所示：

SELECT ST_Union(the_geom) AS geom FROM (
   SELECT chp02.polygon_to_line(geom) AS geom FROM
     chp02.use_area
   ) AS unioned; 

现在，让我们将线字符串转换回多边形，为此，我们可以使用ST_polygonize对结果进行多边形化，如以下查询所示：

SELECT ST_Polygonize(geom) AS geom FROM (
   SELECT ST_Union(the_geom) AS geom FROM (
     SELECT chp02.polygon_to_line(geom) AS geom FROM
     chp02.use_area
   ) AS unioned 
) as polygonized;

st_Polygonize函数将创建一个多多边形，因此，如果我们要使用它做任何有用的事情，我们需要将它分解为多个单多边形几何体。在执行此操作时，我们不妨在CREATE TABLE语句中执行以下操作：

CREATE TABLE chp02.use_area_alt AS (
   SELECT (ST_Dump(the_geom)).geom AS the_geom FROM (
     SELECT ST_Polygonize(the_geom) AS the_geom FROM (
       SELECT ST_Union(the_geom) AS the_geom FROM (
         SELECT chp02.polygon_to_line(the_geom) AS the_geom        FROM chp02.use_area
       ) AS unioned
     ) as polygonized 
) AS exploded 
); 

我们将对该几何体执行空间查询，因此我们应该创建一个索引，以确保查询的性能良好，如以下查询所示：

CREATE INDEX chp02_use_area_alt_the_geom_gist
   ON chp02.use_area_alt
   USING gist(the_geom); 

为了从原始几何图形中找到合适的表信息并将其应用到生成的几何图形中，我们将执行多边形中的点查询。为此，我们首先需要计算合成几何体的质心：

CREATE TABLE chp02.use_area_alt_p AS
   SELECT ST_SetSRID(ST_PointOnSurface(the_geom), 3734) AS
      the_geom FROM
     chp02.use_area_alt;
ALTER TABLE chp02.use_area_alt_p ADD COLUMN gid serial; 
ALTER TABLE chp02.use_area_alt_p ADD PRIMARY KEY (gid);

一如既往，使用以下查询创建空间索引：

CREATE INDEX chp02_use_area_alt_p_the_geom_gist
   ON chp02.use_area_alt_p
   USING gist(the_geom);

质心然后使用以下查询构造原始表格信息和结果多边形之间的多边形中的点（ST_Intersects）关系：

CREATE TABLE chp02.use_area_alt_relation AS
 SELECT points.gid, cu.location FROM
   chp02.use_area_alt_p AS points,
   chp02.use_area AS cu
     WHERE ST_Intersects(points.the_geom, cu.the_geom);

如果我们查看表的第一行，我们可以看到它将点的标识符链接到相应的位置：

它的工作原理。。。
我们的基本方法是查看几何体的底层拓扑，并重建一个不重叠的拓扑，然后使用新几何体的质心构造一个查询，以建立与原始数据的关系。

还有更多。。。
在此阶段，我们可以选择使用外键建立引用完整性框架，如下所示：

ALTER TABLE chp02.use_area_alt_relation ADD FOREIGN KEY (gid) REFERENCES chp02.use_area_alt_p (gid); 

Using polygon overlays for proportional census estimates使用多边形覆盖进行比例普查估计

PostgreSQL函数用于聚集表格数据，包括sum、count、min、max等。作为一个框架，PostGIS没有明确的空间等价物，但这并不妨碍我们使用PostgreSQL的聚合函数与PostGIS的空间功能一起构建函数。

在这个配方中，我们将探讨空间摘要与美国人口普查数据。美国人口普查数据本质上是聚合数据。这样做是为了保护公民的隐私。但当涉及到用这些数据进行分析时，数据的聚合性质可能会成为问题。分解数据有一些技巧。其中最简单的是使用比例和，我们将在本练习中进行。

准备
眼下的问题是，为了向公众提供服务，已经绘制了一条拟议的线索。这个例子可以应用于道路建设，甚至为提供服务而寻找商业地产的场地。
首先，解压缩trail_census.zip文件，然后使用解压缩文件夹中的以下命令执行快速数据加载：

shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom census chp02.trail_census | psql -U me -d postgis_cookbook
shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom trail_alignment_proposed_buffer chp02.trail_buffer | psql -U me -d postgis_cookbook
shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom trail_alignment_proposed chp02.trail_alignment_prop | psql -U me -d postgis_cookbook

上述命令将产生以下输出：

如果我们在我们最喜欢的桌面GIS中查看建议的轨迹，我们有以下内容：

在我们的案例中，我们想了解在小道1英里内的人口，假设在小道1英里以内的人最有可能使用它，因此最有可能由它服务。

为了找出这条建议的小径附近的人口，我们覆盖了人口普查区组人口密度信息。下一个屏幕截图显示的是一个1英里长的缓冲区，围绕着提议的小径：

关于这次普查数据，我们可能会注意到的一点是普查密度和普查区组规模的广泛范围。计算人口的一种方法是简单地选择与我们地区相交的所有普查区块，如以下屏幕截图所示：

这是一个简单的程序，它估计有130至288人居住在一英里以内，但从选择的形状来看，我们可以看到，我们通过在估计中取完整的街区来高估人口。
同样，如果我们仅仅使用了在我们提议的路线路线1英里内的块组，我们将低估人口。
相反，我们将作出一些有用的假设。区块群设计为在块组种群分布中中等均匀。假设数据是正确的，我们可以假设对于给定的块组，如果块组的50%在我们的目标区域内，我们可以将该块组的一半人口归因于我们的估计。将此应用于所有块组，并将它们相加，我们有一个精良的估计，这可能比纯相交或质心查询更好。因此，我们采用比例和。

怎么做。。。
由于比例和问题是一个一般性问题，它可以应用于许多问题。我们将把基本比例写为一个函数。函数接受输入并返回值。在本例中，我们希望比例函数采用两种几何图形，即缓冲轨迹和块组的几何图形以及我们希望比例的值，并希望它返回比例值：

CREATE OR REPLACE FUNCTION chp02.proportional_sum(geometry,   geometry, numeric)
   RETURNS numeric AS
 $BODY$
 -- SQL here
 $BODY$
   LANGUAGE sql VOLATILE; 

现在，为了计算的目的，对于缓冲区和块组的任何给定交叉点，我们要找到交叉点在整个块组上的比例。然后这个值应该乘以我们要缩放的值。

在SQL中，函数类似于以下查询：

SELECT $3 * areacalc FROM
   (SELECT (ST_Area(ST_Intersection($1, $2)) / ST_Area($2)):: numeric AS areacalc
   ) AS areac; 

前面的查询的完整形式如下所示：

CREATE OR REPLACE FUNCTION chp02.proportional_sum(geometry,   geometry, numeric) 
  RETURNS numeric AS
 $BODY$
     SELECT $3 * areacalc FROM
       (SELECT (ST_Area(ST_Intersection($1,           $2))/ST_Area($2))::numeric AS areacalc 
      ) AS areac ; 
$BODY$
   LANGUAGE sql VOLATILE; 

它的工作原理。。。
由于我们已经将查询作为函数编写，所以查询使用SELECT语句循环遍历所有可用记录，并为我们提供一个比例的填充。精明的读者会注意到，我们还没有做任何总结工作；我们只研究了问题的比例部分。我们可以使用PostgreSQL的内置聚合函数调用函数来完成摘要。这种方法的优点是，我们不仅需要应用一个和，而且还可以计算其他聚合，例如min或max。在下面的示例中，我们将只应用一个总和：

SELECT ROUND(SUM(chp02.proportional_sum(a.the_geom, b.the_geom, b.pop))) FROM
   chp02.trail_buffer AS a, chp02.trail_census as b
   WHERE ST_Intersects(a.the_geom, b.the_geom)
   GROUP BY a.gid; 

返回的值非常不同（人口96081），这更可能是准确的。

Working with Vector Data – The Basics

本章，我们将介绍：

• 使用GPS数据
• 修复无效的几何图形
• 基于空间连接的GIS分析
• 简化几何图形
• 测量距离
• 使用公共属性合并多边形
• 计算交叉点
• 剪裁几何图形以部署数据
• 用PostGIS拓扑简化几何图形

介绍
在本章中，您将使用一组PostGIS函数和矢量数据集。您将首先了解如何将PostGIS与GPS数据结合使用您将使用ogr2ogr导入此类数据集，然后使用ST_MakeLine函数从点几何组成多段线。
然后，您将看到PostGIS如何使用诸如ST_MakeValid、ST_IsValid、ST_IsValidReason和ST_IsValidDetails之类的函数来帮助您查找和修复无效的几何体。
然后，您将了解空间数据库中最强大的元素之一，空间连接。PostGIS为您提供了一组丰富的运算符，例如ST_Intersects、ST_Contains、ST_Covers、ST_crosss和ST_DWithin。
之后，当您不需要太多细节时，您将使用ST_Simplify和ST_SimplifyPreverveTopology函数来简化（概括）几何体。该函数适用于线性几何体，多边形几何体可能会引入拓扑异常。在这种情况下，您应该考虑使用外部GIS工具，例如GRASS。
然后，您将浏览PostGIS函数以进行距离测量-ST_distance、ST_DistanceSphere和ST_DistanceSpheroid正在进行中。
本章介绍的其中一个方法将指导您通过典型的GIS工作流，基于公共属性合并多边形；为此，您将使用ST_Union函数。
然后，您将学习如何使用STèu交集函数剪裁几何图形，然后再深入了解PostGIS拓扑（在版本2.0中介绍的最后一个方法）。

Working with GPS data使用GPS数据

在此配方中，您将使用GPS数据。这种数据通常保存在.gpx文件中。您将从RunKeeper（一个为跑步者提供的大众社交网络）将一堆.gpx文件导入PostGIS。
如果你在RunKeeper上有帐户，你可以导出你的.gpx文件并按照这个配方中的说明进行处理。否则，您可以使用RunKeeper-gpx.zip文件中包含的RunKeeper.gpx文件，该文件位于本书的代码包中的chp03目录中。
首先使用ogr2ogr创建一个bash脚本，用于将.gpx文件导入PostGIS表。导入完成后，您将尝试编写两个SQL查询，并测试一些非常有用的函数，例如ST琰MakeLine用于从点几何图形生成多段线，ST琰Length用于计算距离，ST琰Intersects用于执行空间连接操作。

准备
将data/chp03/runkeeper-gpx.zip文件解压缩到working/chp03/runkeeper_gpx。如果您还没有阅读第1章，将数据移入和移出PostGIS，请确保PostGIS数据库中有国家数据集。

怎么做。。。
首先，确定需要导入PostGIS的.gpx文件的格式。打开其中一个并检查文件结构每个文件必须是由一个<trk>元素组成的XML格式，该元素只包含一个<trkseg>元素，该元素包含许多<trkpt>元素（从跑步者的GPS设备存储的点）。将这些点导入PostGIS点表：
使用以下命令创建一个名为chp03的新模式，以存储本章中所有配方的数据：

postgis_cookbook=# create schema chp03;

通过执行以下命令行，在PostgreSQL中创建chp03.rk_track_points表：

postgis_cookbook=# CREATE TABLE chp03.rk_track_points
 (
         fid serial NOT NULL,
         the_geom geometry(Point,4326),
         ele double precision,
         "time" timestamp with time zone,
         CONSTRAINT activities_pk PRIMARY KEY (fid)
); 

创建以下脚本以使用GDAL ogr2ogr命令导入chp03.rk_track_points表中的所有.gpx文件。

以下是Linux版本（将其命名为working/chp03/import_gpx.sh）：

         #!/bin/bash
         for f in `find runkeeper_gpx -name *.gpx -printf "%f
"
         do
           echo "Importing gpx file $f to chp03.rk_track_points
                 PostGIS table..." #, ${f%.*}"
           ogr2ogr -append -update  -f PostgreSQL
           PG:"dbname='postgis_cookbook' user='me'
           password='mypassword'" runkeeper_gpx/$f
           -nln chp03.rk_track_points
           -sql "SELECT ele, time FROM track_points"
         done 

以下是macOS的命令（将其命名为working/chp03/import_gpx.sh）：

        #!/bin/bash
         for f in `find runkeeper_gpx -name *.gpx `
         do
           echo "Importing gpx file $f to chp03.rk_track_points
                 PostGIS table..." #, ${f%.*}"
           ogr2ogr -append -update  -f PostgreSQL
          PG:"dbname='postgis_cookbook' user='me'
          password='mypassword'" $f
           -nln chp03.rk_track_points
           -sql "SELECT ele, time FROM track_points"
         done 

以下是Windows版本（将其命名为working/chp03/import_gpx.bat）：

        @echo off
         for %%I in (runkeeper_gpx*.gpx*) do ( 
          echo Importing gpx file %%~nxI to chp03.rk_track_points
                PostGIS table...
           ogr2ogr -append -update -f PostgreSQL
          PG:"dbname='postgis_cookbook' user='me'
          password='mypassword'" runkeeper_gpx/%%~nxI
           -nln chp03.rk_track_points
           -sql "SELECT ele, time FROM track_points"
         ) 

在Linux和macOS中，不要忘记在运行脚本之前为脚本分配执行权限。然后，运行以下脚本：

      $ chmod 775 import_gpx.sh
      $ ./import_gpx.sh
      Importing gpx file 2012-02-26-0930.gpx to chp03.rk_track_points
         PostGIS table...
      Importing gpx file 2012-02-29-1235.gpx to chp03.rk_track_points
         PostGIS table...
      ...
      Importing gpx file 2011-04-15-1906.gpx to chp03.rk_track_points
         PostGIS table...

在Windows中，只需双击.bat文件或使用以下命令从命令提示符运行它：

> import_gpx.bat

现在，使用ST_MakeLine函数创建一个包含单个流道轨迹细节的多段线表。假设每一天，跑步者只有一次训练。在此表中，您应包括轨道详细信息的开始和结束时间，如下所示：

    postgis_cookbook=# SELECT
       ST_MakeLine(the_geom) AS the_geom,
       run_date::date,
       MIN(run_time) as start_time,
       MAX(run_time) as end_time
       INTO chp03.tracks
       FROM (
         SELECT the_geom,
         "time"::date as run_date,
         "time" as run_time
         FROM chp03.rk_track_points
         ORDER BY run_time
       ) AS foo GROUP BY run_date; 

在查询已创建的表之前，不要忘记向这两个表添加空间索引以提高其性能，如下所示：

postgis_cookbook=# CREATE INDEX rk_track_points_geom_idx
       ON chp03.rk_track_points USING gist(the_geom);
postgis_cookbook=# CREATE INDEX tracks_geom_idx
       ON chp03.tracks USING gist(the_geom);

如果您尝试在任何给定的一天打开桌面GIS上的两个空间表，您应该会看到rkU trackU points表中的点在tracks表中构成一个多段线几何体记录，如以下屏幕截图所示：

如果从桌面GIS（如QGIS）打开所有轨迹，将看到以下内容：

现在，查询tracks表以获得跑步者每月的总跑步距离（以公里为单位）报告。为此，请使用ST_Length函数，如下查询所示：

postgis_cookbook=# SELECT
         EXTRACT(year FROM run_date) AS run_year, 
        EXTRACT(MONTH FROM run_date) as run_month,
         SUM(ST_Length(geography(the_geom)))/1000 AS distance
         FROM chp03.tracks
         GROUP BY run_year, run_month ORDER BY run_year, run_month; 

使用tracks和countries表之间的空间连接，并再次使用STïLength函数（如下所示），您将得到每个国家的跑步者的距离报告（以km为单位）：

postgis_cookbook=# SELECT
         c.country_name,
         SUM(ST_Length(geography(t.the_geom)))/1000 AS run_distance
       FROM chp03.tracks AS t
       JOIN chp01.countries AS c
       ON ST_Intersects(t.the_geom, c.the_geom)
      GROUP BY c.country_name
       ORDER BY run_distance DESC; 

它的工作原理。。。
gpx文件将所有点的详细信息存储在WGS 84空间参考系统中；因此，我们用SRID（4326）创建了rku tracku points表。
在创建rkU tracku points表之后，我们使用bash脚本导入了runkeeperu gpx目录中的所有.gpx文件。bash脚本迭代runkeeperu gpx目录中扩展名为*.gpx的所有文件。对于每个文件，脚本都运行ogr2ogr命令，使用gpx GDAL驱动程序将.gpx文件导入PostGIS（有关更多详细信息，请转到http://www.gdal.org/drv_gpx.html).
在GDAL的抽象中，.gpx文件是由以下几层组成的OGR数据源：

 在.gpx文件（OGR数据源）中，只有轨迹和轨迹点图层。作为一种快捷方式，您可以使用ogr2ogr只导入tracks层，但是您需要从tracku points层开始，以便使用一些PostGIS函数生成tracks层本身。这就是为什么在bash脚本的ogr2ogr部分中，我们将tracku points图层中的点几何图形以及一些有用的属性（如高程和时间戳）导入到rku tracku points PostGIS表中。

导入记录后，我们使用子查询输入一个名为tracks的新多段线表，并从rk track points表中选择所有点几何图形及其日期和时间，按日期分组，并使用ST U MakeLine函数聚合几何图形。此函数能够从点几何图形创建线串（有关更多详细信息，请转到http://www.postgis.org/docs/ST_MakeLine.html).

不要忘记按日期时间对子查询中的点进行排序；否则，您将获得一个不规则的线串，从一点跳到另一点，并且没有遵循正确的顺序。

加载tracks表之后，我们测试了两个空间查询。
一开始，你得到了一份逐月的跑步者总距离的报告。为此，您选择了按日期（年和月）分组的所有轨迹记录，总距离是通过将单个轨迹的长度相加得到的（使用ST U长度函数获得）。为了从runu date函数中获取年份和月份，您使用了PostgreSQL EXTRACT函数。请注意，如果您使用WGS 84系统中的几何图形测量距离，您将以度为单位获得距离。因此，必须将几何图形投影到为特定区域设计的平面公制系统，从该区域投影数据。
对于大范围区域，例如在我们的例子中，我们有遍布四周的点

如上一个查询结果所示，一个好的选择是使用postgis1.5引入的geography数据类型。计算速度可能较慢，但比其他系统精确得多。这就是在进行测量之前将几何体强制转换为地理数据类型的原因。
最后一个空间查询使用了一个带有STu Intersects函数的空间连接来获取运行者运行的每个轨迹的国家名称（假设运行者没有运行跨境轨迹）。获取每个国家的总跑步距离只需聚合countryu name字段上的选择，并使用PostgreSQL SUM运算符聚合赛道距离。

Fixing invalid geometries修复无效的几何图形

Into the Nth Dimension

本章，我们将介绍：

• 导入LiDAR数据
• 在LiDAR点云数据上执行3D查询
• 构建和服务2.D建筑物
• 使用ST_Extrude来挤出建筑物
• 为PostGIS创建任意三维对象
• 为网络平台导出X3D模型
• 利用PostGIS三维重建无人机足迹图像
• PostGIS中的无人机摄影测量-点云
• PostGIS中的无人机摄影测量–DSM创建

介绍
在本章中，我们将探讨PostGIS的3D功能。我们将重点讨论三大类：如何将三维数据插入PostGIS，如何使用三维数据分析和执行查询，以及如何从PostGIS中转储三维数据。本章将使用三维点云作为三维数据，包括激光雷达数据和从运动结构（SfM）技术得到的数据。此外，我们将建立一个功能，挤出建筑足迹三维。
需要注意的是，在本章中，我们将讨论postgreSQL点云扩展；点云通常是坐标系中点坐标的三维表示的大型数据集。点云是用来表示物体表面的高度准确，如使用地理激光雷达数据。点云扩展将帮助我们将激光雷达数据存储到数据库中的点云对象中。此外，此扩展还添加了一些函数，允许您将点云对象转换为几何图形，并使用点云数据进行空间过滤。有关此扩展的更多信息，请访问官方GitHub存储库https://github.com/pgpointcloud/pointcloud. 此外，你还可以在http://workshops.boundlessgeo.com/tutorial-lidar/.
下载我们的示例数据集供您使用，请访问http://www.packtpub.com/support.

Importing LiDAR data导入激光雷达数据

光探测和测距（LiDAR）是产生点云数据最常用的设备之一。该系统捕获给定空间中对象或曲面的三维位置和其他特性。这种方法与雷达非常相似，因为它使用电磁波来测量距离和亮度等。然而，激光雷达和雷达的一个主要区别是前者使用激光束技术，而不是微波或无线电波。另一个区别是，激光雷达通常发出一个聚焦脉冲，测量返回脉冲的时间，计算距离和深度。相比之下，雷达在接收返回脉冲之前会发出多个脉冲，因此需要额外的处理来确定每个脉冲的来源。
激光雷达数据在地面和空中应用中已变得相当普遍，有助于地面测量，提高摄影测量工程的自动化程度。激光雷达数据来源广泛，数据量大。
激光雷达数据通常分布在LASor LASerformat上，美国摄影测量和遥感学会（ASPRS）制定了LAS标准。LAS是一种二进制格式，因此将其读入PostGIS数据库并非易事。幸运的是，我们可以使用开源工具PDAL。

准备
我们的源数据将采用LAS格式，我们将使用PDAL库将其插入我们的数据库，可在https://www.pdal.io/. 此工具可用于Linux/UNIX和Mac用户；对于Windows，它与OSGeo4W包一起提供(https://www.pdal.io/workshop/osgeo4w.html).
LAS数据可以包含许多有趣的数据，而不仅仅是X、Y和z值。它可以包括从感测物体返回的强度和物体的分类（地面、植被和建筑物）。当我们在PostGIS数据库中放置LAS文件时，我们可以选择性地收集这些信息。此外，PDAL在内部构造了一个管道来转换数据，以便进行读取、处理和写入。
为此，我们需要创建一个表示PDAL处理管道的JSON文件。对于每个LAS文件，我们创建一个JSON文件来配置读写器以使用postgres pointcloud选项。我们还需要编写数据库连接参数。对于测试文件testu 1.las，代码如下：

 现在，我们可以下载我们的数据。建议从以下位置下载http://gis5.oit.ohio.gov/geodatadownload/or 下载示例数据集供您使用，请访问http://www.packtpub.com/support.

怎么做。。。
首先，我们需要将LAS文件转换为PDAL可以使用的格式。我们创建了dapython脚本，它从LAS文件目录中读取并生成相应的JSON。使用这个脚本，如果我们有一个大目录的文件，我们可以自动生成。另外，我们选择Python是因为它的简单性，并且因为无论您使用的操作系统如何，您都可以执行脚本。要执行脚本，请在控制台中运行以下操作（对于Windows用户，请确保路径变量中包含Python解释器）：
$python insertu files.py-f<lasfilesu path>
这个脚本将读取每个LAS文件，并将存储在一个名为pipelines的文件夹中，所有与LAS文件相关的元数据将插入数据库。
现在，使用PDAL，我们执行for循环，将LAS文件插入Postgres：

$ for file in `ls pipelines/*.json`;
   do
     pdal pipeline $file;
   done 

这个点云数据被分成三个不同的表。如果要合并它们，则需要执行以下SQL命令：

DROP TABLE IF EXISTS chp07.lidar; 
CREATE TABLE chp07.lidar AS WITH patches AS 
 (
   SELECT 
    pa 
  FROM "chp07"."N2210595"
    UNION ALL
   SELECT
     pa
    FROM "chp07"."N2215595"
    UNION ALL
   SELECT
     pa
    FROM "chp07"."N2220595"
 )
 SELECT
   2 AS id,
 PC_Union(pa) AS pa
  FROM patches;

Performing 3D queries on a LiDAR point cloud在激光雷达点云上执行三维查询
在前面的方法中，导入LiDAR数据，我们将LiDAR 3D点云引入PostGIS，从输入创建一个显式3D数据集。有了三维形式的数据，我们就能够对其执行空间查询。在这个方法中，我们将利用3D索引，以便我们的最近邻搜索在数据所在的所有维度上都能工作。

怎么做。。。
我们将使用前面配方中导入的激光雷达数据作为我们选择的数据集。我们把那张表命名为chp07.lidar。要执行最近邻搜索，我们需要在数据集上创建索引。空间索引与普通的数据库表索引非常相似，与book索引相似，因为它们可以帮助我们更快地找到所需内容。通常，这样的索引创建步骤如下所示（我们这次不运行）：

CREATE INDEX chp07_lidar_the_geom_idx 
ON chp07.lidar USING gist(the_geom);

对于二维查询，三维索引的执行速度不如二维索引，因此创建索引查询默认为创建二维索引。在我们的例子中，我们希望强制gist应用于所有三个维度，因此我们将显式地告诉PostgreSQL使用索引的n维版本：

CREATE INDEX chp07_lidar_the_geom_3dx ON chp07.lidar USING gist(the_geom gist_geometry_ops_nd);

请注意，如果我们有一个时间维度或3D+time，那么前面代码中描述的方法也可以工作。让我们加载第二个3D数据集和将在查询中使用的流中心线：

$ shp2pgsql -s 3734 -d -i -I -W LATIN1 -t 3DZ -g the_geom hydro_line chp07.hydro | PGPASSWORD=me psql -U me -d "postgis-cookbook" -h localhost

如下图所示，这些数据与我们的激光雷达点云很好地重叠：

 

Constructing and serving buildings 2.5D构建和服务2.5D建筑物

在第四章LiDARrecipe的详细建筑足迹中，我们使用矢量数据-高级配方，探索了利用LiDAR数据自动生成建筑足迹的方法。当时我们要做的是从三维数据中创建二维数据。但在此次这个食谱中，我们尝试相反的方法，在某种意义上。我们从建筑足迹的二维多边形开始，将它们输入到一个函数中，该函数将它们挤出为三维多边形。

准备
对于此次配方，我们将挤出自己制作的建筑足迹。让我们快速创建一个具有单个建筑示意图的表，以便进行测试，如下所示：

DROP TABLE IF EXISTS chp07.simple_building; 
CREATE TABLE chp07.simple_building AS
SELECT 1 AS gid, ST_MakePolygon(
   ST_GeomFromText(
     'LINESTRING(0 0,2 0, 2 1, 1 1, 1 2, 0 2, 0 0)'
   ) 
) AS the_geom; 

将3D建筑的创建尽可能简单地封装在一个函数中是有益的：

CREATE OR REPLACE FUNCTION chp07.threedbuilding(footprint geometry, height numeric) 
RETURNS geometry AS 
$BODY$

我们的函数需要两个输入：建筑足迹和要拉伸到的高度。我们还可以想象一个函数包含第三个参数：建筑物底部的高度。

要构建建筑墙，我们需要首先将多边形转换为线串，然后将线串进一步分离为各自的两点线段：

WITH simple_lines AS 
(
   SELECT
      1 AS gid,
      ST_MakeLine(ST_PointN(the_geom,pointn), 
    ST_PointN(the_geom,pointn+1)) AS the_geom
   FROM (
     SELECT 1 AS gid,
     polygon_to_line($1) AS the_geom
) AS a
   LEFT JOIN(
     SELECT
        1 AS gid,
        generate_series(1,
          ST_NumPoints(polygon_to_line($1))-1
       ) AS pointn
    ) AS b
   ON a.gid = b.gid ), 

前面的代码返回原始形状的每个两点段。例如，对于简单的构建，输出如下：

现在我们有了一系列单独的线条，我们可以用它们来建造建筑物的墙壁。首先，我们需要使用ST_Force3DZ将我们的2D线条重铸为3D：

threeDlines AS(
   SELECT ST_Force3DZ(the_geom) AS the_geom FROM simple_lines
 ),

输出如下：

 下一步是将这些行中的每一行从MULTILINESTRING分解为许多LINESTRINGS：

explodedLine AS (
   SELECT (ST_Dump(the_geom)).geom AS the_geom FROM threeDLines
), 

其输出如下：

 下一步是构造表示拉伸墙边界的线：

threeDline AS (
   SELECT ST_MakeLine(
     ARRAY[
       ST_StartPoint(the_geom),
       ST_EndPoint(the_geom),
       ST_Translate(ST_EndPoint(the_geom), 0, 0, $2),
       ST_Translate(ST_StartPoint(the_geom), 0, 0, $2),
       ST_StartPoint(the_geom)
     ]
   ) 
  AS the_geom FROM explodedLine
 ), 

现在，我们需要将每个linestring转换为polygon.threeDwall：

threeDwall AS
(
   SELECT ST_MakePolygon(the_geom) as the_geom FROM threeDline
 ), 

最后，使用地板的原始几何图形（强制为3D）和转换为输入高度的原始几何图形副本，将屋顶和地板放置在我们的建筑中：

buildingTop AS
 (
   SELECT ST_Translate(ST_Force3DZ($1), 0, 0, $2) AS the_geom
 ),
 -- and a floor
 buildingBottom AS
 (
   SELECT ST_Translate(ST_Force3DZ($1), 0, 0, 0) AS the_geom
 ), 

我们将墙、屋顶和地板放在一起，在过程中，将其转换为3D MULTIPOLYGON：

wholeBuilding A
 (
   SELECT the_geom FROM buildingBottom
     UNION ALL
   SELECT the_geom FROM threeDwall
     UNION ALL
   SELECT the_geom FROM buildingTop
 ), 
-- then convert this collecion to a multipolygon 
multiBuilding AS (
   SELECT ST_Multi(ST_Collect(the_geom)) AS the_geom FROM 
     wholeBuilding
 ), 

虽然我们可以把我们的几何学作为一个MULTIPOLYGON多面体，但我们会做得很好，并对POLYHEDRALSURFACE多面体曲面进行非正式的转换。在我们的例子中，我们已经被有效地格式化为多面体曲面，所以我们只需将几何体转换为文本ST_AsText，然后将单词替换为多面体曲面，然后将文本转换回几何体，并使用ST_GeomFromText：

textBuilding AS 
(
   SELECT ST_AsText(the_geom) textbuilding FROM multiBuilding
),
textBuildSurface AS 
(
   SELECT ST_GeomFromText(replace(textbuilding, 'MULTIPOLYGON',          'POLYHEDRALSURFACE')) AS the_geom FROM textBuilding 
) 
SELECT the_geom FROM textBuildSurface 

最后，整个功能是：

CREATE OR REPLACE FUNCTION chp07.threedbuilding(footprint geometry,
    height numeric) 
RETURNS geometry AS
 $BODY$
  -- make our polygons into lines, and then chop up into individual line segments 
WITH simple_lines AS (
   SELECT 1 AS gid, ST_MakeLine(ST_PointN(the_geom,pointn),     ST_PointN(the_geom,pointn+1)) AS the_geom
   FROM (SELECT 1 AS gid, polygon_to_line($1) AS the_geom ) AS a
   LEFT JOIN   (SELECT 1 AS gid, generate_series(1,      ST_NumPoints(polygon_to_line($1))-1) AS pointn
    ) AS b
   ON a.gid = b.gid ),
 -- convert our lines into 3D lines, which will set our third  coordinate to 0 by default
 threeDlines AS (
   SELECT ST_Force3DZ(the_geom) AS the_geom FROM simple_lines 
), 
-- now we need our lines as individual records, so we dump them out using ST_Dump, and then just grab the geometry portion of the dump
 explodedLine AS (
   SELECT (ST_Dump(the_geom)).geom AS the_geom FROM threeDLines 
), 
-- Next step is to construct a line representing the boundary of the  extruded "wall" 
threeDline AS
 (
   SELECT ST_MakeLine(
     ARRAY[
     ST_StartPoint(the_geom),
     ST_EndPoint(the_geom),
     ST_Translate(ST_EndPoint(the_geom), 0, 0, $2),
     ST_Translate(ST_StartPoint(the_geom), 0, 0, $2),
     ST_StartPoint(the_geom)
     ]
   ) 
AS the_geom FROM explodedLine
 ), 
-- we convert this line into a polygon
 threeDwall AS
 ( 
  SELECT ST_MakePolygon(the_geom) as the_geom FROM threeDline
 ),
 -- add a top to the building buildingTop AS ( 
  SELECT ST_Translate(ST_Force3DZ($1), 0, 0, $2) AS the_geom ), 
-- and a floor
 buildingBottom AS
 (
   SELECT ST_Translate(ST_Force3DZ($1), 0, 0, 0) AS the_geom
 ),
 -- now we put the walls, roof, and floor together
 wholeBuilding AS
 (
   SELECT the_geom FROM buildingBottom
     UNION ALL
   SELECT the_geom FROM threeDwall
     UNION ALL
   SELECT the_geom FROM buildingTop
 ),
 -- then convert this collecion to a multipolygon 
multiBuilding AS (
   SELECT ST_Multi(ST_Collect(the_geom)) AS the_geom FROM wholeBuilding
 ), 
-- While we could leave this as a multipolygon, we'll do things properly and munge an informal cast 
-- to polyhedralsurfacem which is more widely recognized as the appropriate format for a geometry like 
-- this. In our case, we are already formatted as a polyhedralsurface, minus the official designation,
-- so we'll just convert to text, replace the word MULTIPOLYGON with POLYHEDRALSURFACE and then convert 
-- back to geometry with ST_GeomFromText 
 textBuilding AS (
   SELECT ST_AsText(the_geom) textbuilding FROM multiBuilding
 ),
 textBuildSurface AS (
   SELECT ST_GeomFromText(replace(textbuilding, 'MULTIPOLYGON',      'POLYHEDRALSURFACE')) AS the_geom FROM textBuilding 
) 
SELECT the_geom FROM textBuildSurface ;
 $BODY$
   ANGUAGE sql VOLATILE
   COST 100; 
ALTER FUNCTION chp07.threedbuilding(geometry, numeric)
OWNER TO me; 

怎么做。。。

现在我们有了一个三维建筑拉伸功能，我们可以很容易地用我们的封装良好的功能拉伸我们的建筑足迹：

DROP TABLE IF EXISTS chp07.threed_building; 
CREATE TABLE chp07.threed_building AS  
SELECT chp07.threeDbuilding(the_geom, 10) AS the_geom  
FROM chp07.simple_building; 

我们可以将此函数应用于真实的建筑迹线数据集（在我们的数据目录中提供），在这种情况下，如果我们有一个高度字段，我们可以根据它进行拉伸：

shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom building_footprintsc
hp07.building_footprints | psql -U me -d postgis-cookbook 
-h <HOST> -p <PORT> 
DROP TABLE IF EXISTS chp07.build_footprints_threed; 
CREATE TABLE chp07.build_footprints_threed AS  
SELECT gid, height, chp07.threeDbuilding(the_geom, height) AS the_geom  
FROM chp07.building_footprints;

生成的输出为我们提供了一组漂亮的拉伸建筑示意图，如下图所示：

 第四章详细介绍了LiDARrecipe中的建筑物足迹，《利用矢量数据-高级配方》，探索了从LiDAR中提取建筑物足迹的方法。可以设想一个完整的工作流程，从激光雷达中提取建筑物足迹，然后使用当前配方重建多边形几何体，从而将点云转换为曲面，将当前配方与前面引用的配方相结合。

Using ST_Extrude to extrude building footprints使用st_Extrude拉伸建筑示意图

postgis2.1为PostGIS带来了许多非常酷的附加功能。对PostGIS光栅类型的操作是PostGIS 2.1中更重要的改进之一。一个更安静和同样有效的游戏改变者是SFCGAL库作为PostGIS的可选扩展的增加。该网站称http://sfcgal.org/SFCGAL是围绕CGAL的C++包装库，目的是支持ISO 19107:2013和OGC简单特性1.2访问3D操作。

从实际的角度来看，这意味着什么？这意味着PostGIS正从几何图形本身的表示和对这些三维几何图形的操作走向一个全功能的三维环境。更多信息请访问http://postgis.net/docs/reference.html#reference_sfcgal.

这和其他一些方法将假设您安装了一个版本的PostGIS，并编译和启用了SFCGAL。这样做将启用以下功能：

ST_Extrude: This extrudes a surface to a related volume将曲面挤出到相关体积

ST_StraightSkeleton: This computes a straight skeleton from a geometry从几何体计算直线骨架

ST_IsPlanar: This checks whether a surface is a planar or not检查曲面是否为平面

ST_Orientation: This determines the surface orientation确定曲面方向

ST_ForceLHR: This forces LHR orientation强制左后定向

ST_MinkowskiSum: This computes the Minkowski sum计算Minkowski和

ST_Tesselate: This performs surface Tessellation执行曲面细分

怎么做。。。
对于这个配方，我们将使用ST_Extrude，就像我们在上一个配方中使用我们自己的定制功能一样，构建和服务建筑2.5D。与前面的方法相比，它的优点是我们不需要在PostGIS中编译SFCGAL库。这个配方的优点是我们对挤压过程有更多的控制；也就是说，我们可以在所有的三维空间中挤出。

st_Extrude返回一个几何体，特别是一个多面体曲面。它需要四个参数：输入几何体和沿X、Y和Z轴的拉伸量：

DROP TABLE IF EXISTS chp07.buildings_extruded; 
CREATE TABLE chp07.buildings_extruded AS  
SELECT
 gid, ST_CollectionExtract(ST_Extrude(the_geom, 20, 20, 40), 3) as the_geom
FROM chp07.building_footprints 

 因此，在2.5DREIPE建筑的建设和服务的帮助下，我们得到了挤压建筑，但具有一定的灵活性。

Creating arbitrary 3D objects for PostGIS为PostGIS创建任意的3D对象

三维信息源不仅是由激光雷达生成的，也不是纯粹由二维几何图形和相关属性合成的，如《构建与服务建筑物2.5D》和《使用ST_Extrude挤出建筑足迹》中所述那样，它们也可以根据计算机视觉原理创建。从图像之间相关关键点的关联计算3D信息的过程称为SfM。

作为一种计算机视觉概念，我们可以利用SfM以类似于人脑在3D中感知世界的方式生成3D信息，并进一步将这些信息存储和处理到PostGIS数据库中。

计算机视觉是计算机科学中的一门学科，主要研究图像和视频的自动分析和推理。它被认为是一个研究领域，发展算法，解释世界的方式类似于人类的视觉。一个优秀的总结可以在http://en.wikipedia.org/wiki/Computer_vision.

许多开源项目已经成熟，可以解决SfM问题。其中流行的是Bundler，可以在http://phototour.cs.washington.edu/bundler/，和VisualSFMathttp://ccwu.me/vsfm/. 二进制文件可用于这些工具的多个平台，包括版本。这类项目的好处在于，一组简单的照片可以用来重建3D场景。

出于我们的目的，我们将使用VisualSFM并跳过此软件的安装和配置。原因是SfM超出了PostGIS书籍的详细介绍范围，我们将重点介绍如何在PostGIS中使用数据。

准备
重要的是要理解，SfM技术虽然非常有效，但在可有效处理为点云的图像类型方面有一定的局限性。这些技术依赖于在随后的图像之间找到匹配，因此在处理平滑的图像、丢失相机嵌入的可交换图像文件格式（EXIF）信息或来自手机相机的图像时会有困难。
EXIF标记是图像的元数据格式。存储在这些标签中的通常是相机设置、相机类型、镜头类型以及与SfM提取相关的其他信息。
我们将开始将一个图像序列处理成一个点云，其中包含一个我们知道基本上可以工作的照片序列，但是当您使用SfM进行实验时，您可以输入您自己的照片序列。关于如何创建一个将产生一个三维模型的照片系列的好技巧可以在https://www.youtube.com/watch?v=IStU-WP2XKs&amp;t=348sand http://www.cubify.com/products/capture/photography_tips.aspx.

怎么做。。。
从下载VisualSFMhttp://ccwu.me/vsfm/. 在控制台终端中，执行以下操作：
Visualsfm<IMAGESu FOLDER>
VisualSFM将开始渲染3D模型，并使用包含图像的文件夹作为输入。这需要几个小时来处理。然后，当它完成时，它将返回一个点云文件。
我们可以在MeshLabat之类的程序中查看这些数据http://meshlab.sourceforge.net/. 有关使用MeshLab查看点云的良好教程，请访问http://www.cse.iitd.ac.in/~mcs112609/Meshlab%20Tutorial.pdf.
下图显示了在MeshLab中查看点云时的外观：

 在VisualSFM输出中，有一个扩展名为.ply的文件，例如长颈鹿.ply（包含在本章的源代码中）。如果在atext编辑器中打开此文件，它将显示如下内容：

这是我们文件的头部分。它指定.ply格式、编码格式ascii 1.0、顶点数，然后为返回的所有数据指定列名：x、y、z、nx、ny、nz、red、green和blue。
对于导入PostGIS，我们将导入所有字段，但将重点放在x、y和Z上，以获得点云，以及查看颜色。为了我们的目的，此文件指定相对x、y和z坐标，以及通道中每个点的颜色（红色、绿色和蓝色）。这些颜色是24位颜色，因此它们可以有介于0到255之间的整数。
对于该配方的其余部分，我们将创建一个PDAL管道，将JSON结构读取器修改为.ply文件。请检查本章中导入LiDAR数据的配方，以了解如何创建PDAL管道：
{  "pipeline": [{    "type": "readers.ply",   

                           "filename": "/data/giraffe/giraffe.ply"  },

                     {    "type": "writers.pgpointcloud",   

                          "connection": "host='localhost' dbname='postgis-cookbook' user='me'      password='me' port='5432'",   

                           "table": "giraffe",   

                           "srid": "3734",   

                           "schema": "chp07"  }]

}

然后我们在终端执行：
$pdal pipeline giraffe.json“
这个输出将作为下一个配方的输入。

Exporting models as X3D for the web将模型导出为web的X3D

如果我们没有能力将数据以某种可用的形式提取出来，那么在PostGIS数据库中输入3D数据就没有那么有趣了。解决这个问题的一种方法是利用PostGIS将3D表写入X3D格式的能力。
X3D是用于显示3D数据的XML标准，在web上运行良好。对于那些熟悉虚拟现实建模语言（VRML）的人来说，X3D是该语言的下一代。
要在浏览器中查看X3D，用户可以选择各种插件，也可以利用JavaScript api来支持查看。我们将执行后者，因为它不需要用户配置就可以工作。我们将使用X3DOM的JavaScript框架来实现这一点。X3DOM演示了HTML5和3D的集成，并使用了Web图形库（WebGL）(https://en.wikipedia.org/wiki/WebGL)允许在浏览器中渲染三维内容并与之交互。这意味着我们的数据将不会在不兼容WebGL的浏览器中显示。

准备
我们将使用上一个示例中的点云以X3D格式提供服务。关于X3D的PostGIS文档包括一个使用STu AsX3D函数输出格式化X3D代码的示例：
复制（用pts作为（从chp07.giraffe中选择PCu Explode（pa）作为pt）选择“<X3D xmlns=”http://www.web3d.org/specifications/x3d-namespaceshowStat=“false”showLog=“false”x=“0px”y=“0px”width=“800px”height=“600px”><Scene><Transform><Shape>'| | ST|u AsX3D（ST|u Union（pt:：geometry））| |'</Shape></Transform></Scene></X3D>'从pts）以CSV格式发送到STDOUT；
我们使用CSV将拷贝包含到STDOUT中，以便在原始代码中进行转储。用户可以将此查询保存为SQL脚本文件，并从控制台执行它，以便将结果转储到文件中。例如：
$psql-U me-d postgis食谱-h localhost-f“x3dU query.sql”>result.html

怎么做。。。
这个示例虽然在纯X3D中已经完成，但是需要额外的代码来允许在浏览器中查看。我们通过包含样式表来实现，相应的X3DOM包含XHTML文档的标题：

<link rel="stylesheet" type="text/css" href="http://x3dom.org/x3dom/example/x3dom.css" />

<script type="text/javascript" src="http://x3dom.org/x3dom/example/x3dom.js">

</script>
The full query to generate the XHTML of X3D data is as follows:
COPY(WITH pts AS ( SELECT PC_Explode(pa) AS pt FROM chp07.giraffe)

SELECT regexp_replace(' <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta http-equiv="X-UA-Compatible" content="chrome=1" />

<meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>

<title>Point Cloud in a Browser</title>

<link rel="stylesheet" type="text/css" href="http://x3dom.org/x3dom/example/x3dom.css" />

<script type="text/javascript"

src="http://x3dom.org/x3dom/example/x3dom.js">     

</script>   

</head>   

<body>     

<h1>Point Cloud in the Browser</h1>     

<p>        Use mouse to rotate, scroll wheel to zoom, and control         (or command) click to pan.      </p>     

<X3D xmlns="http://www.web3d.org/specifications/x3d-namespace       showStat="false" showLog="false" x="0px" y="0px" width="800px"       height="600px">       

<Scene>         

<Transform>           

<Shape>' ||  ST_AsX3D(ST_Union(pt::geometry)) || '</Shape>         

</Transform>       

</Scene>     

</X3D>   

</body> 

</html>', E'[\n\r]+','', 'g')

从pts）到STDOUT；
如果在我们最喜欢的浏览器中打开.html文件，我们将得到以下内容：

还有更多。。。
人们可能希望将此X3D转换用作函数，将几何图形输入函数并获得一个页面作为回报。这样，我们可以轻松地为其他表重用代码。X3D转换在函数中体现如下：