Python-地理空间分析精要-全-

Python 地理空间分析精要（全）

原文：zh.annas-archive.org/md5/138d751064dbfff7254887b8a6a73ad4

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

有几个强大的 Python 库用于读取、处理、分析和查看地理空间数据。还有许多网站提供高质量的地理空间数据，您可以在自己的项目中自由使用。这些数据通常将成为您分析的基础，提供国家形状、城市位置、道路轮廓等。结合可用的地理空间库，您将拥有一个强大的工具集，用于使用 Python 进行自己的地理空间分析。

这本书涵盖的内容

第一章, 地理空间分析和技术，引导读者通过下载示例地理空间数据的过程，然后编写一个简单的 Python 程序来读取和分析这些示例数据。

第二章, 地理空间数据，专注于用于地理空间分析的数据：如何获取它、为什么好的数据很重要、地理空间数据可以采用的不同格式以及如何生成您自己的空间数据集。

第三章, 空间数据库，简要介绍了创建地理空间数据库、如何在空间启用数据库中存储数据以及如何对数据进行高效查询。

第四章, 创建地图，探讨了如何使用 Mapnik 库来制作外观精美的地图。

第五章, 地理空间数据分析，指导读者通过使用 Python 编写空间分析程序的过程。基于在第二章（ch02.html "第二章。地理空间数据"）中下载的数据集，并使用主要的 Python 地理空间分析库，本章采用类似菜谱的格式来解决一系列典型的空间分析问题。

第六章, 构建完整的地理空间分析系统，使用前面章节中介绍的各种库和技术来构建一个完整的地理空间分析系统。

您需要为此书准备的

本书中的代码示例使用 Python 2 来分析地理空间数据。任何运行 Windows、Mac OS X 或 Linux 的合理强大的计算机都将是合适的。您需要在您的计算机上下载并安装以下软件：

Python 版本 2.7 或更高版本，不包括 Python 3.x
GDAL/OGR 版本 1.11 或更高版本
GEOS 版本 3.4.2 或更高版本
Shapely 版本 1.5.7 或更高版本
PostgreSQL 版本 9.3 或更高版本
PostGIS 版本 2.1.4 或更高版本
psycopg2 版本 2.5 或更高版本
Mapnik 版本 2.2 或更高版本
PROJ 版本 4.0 或更高版本
PyProj 版本 1.9.4 或更高版本
NetworkX 版本 1.9.1 或更高版本

本书包含下载、安装和使用这些各种工具和库的完整说明。

这本书面向谁

如果您是一位希望掌握地理空间编程的 Python 开发者，或者有特定的空间编程需求，那么这本书就是为您准备的。虽然熟悉安装第三方 Python 库将是一个优势，但不需要具备地理空间编程概念或技术的先验知识。

惯例

在本书中，您将找到许多文本样式，用于区分不同类型的信息。以下是一些这些样式的示例及其含义的解释。

文本中的代码单词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 昵称如下所示：“安装后，您可以通过启动 Python 解释器并输入import osgeo.gdal然后import osgeo.ogr来检查它是否正常工作。”

代码块如下设置：

import osgeo.ogr
shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    feature_name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    geometry_type = geometry.GetGeometryName()
    print i, feature_name, geometry_type

当我们希望引起您对代码块中特定部分的注意时，相关的行或项目将以粗体显示：

from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
dstFile = driver.CreateDataSource("test-shapefile")

任何命令行输入或输出如下所示：

% python readRaster.py
-500 53081919
-84 1
-83 8
-82 9
-81 17
...
5241 1
5295 1
5300 1
5443 1

新术语和重要词汇以粗体显示。屏幕上看到的单词，例如在菜单或对话框中，在文本中如下所示：“点击下一步按钮将您带到下一屏幕。”

注意

警告或重要注意事项如下所示：

小贴士

小技巧和窍门看起来是这样的。

读者反馈

我们欢迎读者的反馈。告诉我们您对这本书的看法——您喜欢或不喜欢什么。读者反馈对我们来说很重要，因为它帮助我们开发出您真正能从中获得最大收益的标题。

要向我们发送一般反馈，只需发送电子邮件至<feedback@packtpub.com>，并在邮件主题中提及本书的标题。

如果您在某个主题上具有专业知识，并且您有兴趣撰写或为本书做出贡献，请参阅我们的作者指南www.packtpub.com/authors。

客户支持

现在您已经是 Packt 图书的骄傲拥有者，我们有一些事情可以帮助您从您的购买中获得最大收益。

下载示例代码

您可以从www.packtpub.com下载您购买的所有 Packt 出版物的示例代码文件。如果您在其他地方购买了这本书，您可以访问www.packtpub.com/support并注册，以便将文件直接通过电子邮件发送给您。

错误清单

尽管我们已经尽一切努力确保我们内容的准确性，但错误仍然会发生。如果您在我们的某本书中发现错误——可能是文本或代码中的错误——如果您能向我们报告这一点，我们将不胜感激。通过这样做，您可以节省其他读者的挫败感，并帮助我们改进本书的后续版本。如果您发现任何勘误，请通过访问www.packtpub.com/submit-errata，选择您的书籍，点击勘误提交表单链接，并输入您的勘误详情来报告它们。一旦您的勘误得到验证，您的提交将被接受，勘误将被上传到我们的网站或添加到该标题的勘误部分下的现有勘误列表中。

要查看之前提交的勘误表，请访问www.packtpub.com/books/content/support，并在搜索字段中输入书籍名称。所需信息将出现在勘误部分。

盗版

在互联网上对版权材料的盗版是一个跨所有媒体的持续问题。在 Packt，我们非常重视保护我们的版权和许可证。如果您在互联网上发现任何形式的我们作品的非法副本，请立即提供位置地址或网站名称，以便我们可以寻求补救措施。

请通过<copyright@packtpub.com>与我们联系，并提供疑似盗版材料的链接。

我们感谢您在保护我们作者和我们提供有价值内容的能力方面的帮助。

询问

如果您对本书的任何方面有问题，您可以通过<questions@packtpub.com>联系我们，我们将尽力解决问题。

第一章：地理空间分析和技术

在本章的介绍中，我们将通过了解您通常将执行的任务类型来开始我们对地理空间分析的探索，然后查看空间数据以及您可以用来处理它的 Python 库。我们将通过编写一个用于分析一些地理空间数据的 Python 示例程序来结束。

随着您在本章的学习，您将：

熟悉地理空间分析将帮助解决的问题类型
了解各种类型的地理空间数据以及与位置数据相关的一些重要概念
设置您的计算机以使用您开始使用 Python 分析地理空间数据所需的第三方库
获取一些基本的地理空间数据以开始
学习如何使用 GDAL/OGR 库读取 shapefile 并提取每个要素的属性和几何形状
学习如何使用 Shapely 操作和分析地理空间数据
编写一个简单但完整的程序来识别邻近国家

让我们先看看通常使用地理空间分析解决的问题和任务类型。

关于地理空间分析

地理空间分析是读取、操作和总结地理空间数据以产生有用和有趣结果的过程。很多时候，您将回答如下问题：

苏萨利托和棕榈泉之间的最短驾驶距离是多少？
法国和比利时之间的边界总长度是多少？
新西兰沿海的每个国家公园的面积是多少？

这些问题的答案通常是一个或一系列数字。其他类型的地理空间分析将涉及根据现有数据计算新的地理空间数据集。例如：

计算从加利福尼亚州洛杉矶到新墨西哥州阿尔伯克基的 USA 66 的高程剖面。
展示我巴西赤道以北的部分。
突出显示如果海洋上升 2 米，拉罗汤加可能被淹没的区域。

在这些情况下，您将生成一组新的地理空间数据，您通常会将其显示在图表或地图上。

要执行此类分析，您需要两样东西：适当的地理空间分析工具和合适的地理空间数据。

我们很快将进行一些简单的地理空间分析。在我们这样做之前，让我们更仔细地看看地理空间数据的概念。

理解地理空间数据

地理空间数据是在地球表面上定位事物的数据。这是一个故意模糊的定义，它涵盖了位置和形状的概念。例如，一个车祸数据库可能包括标识每次事故发生位置的纬度和经度坐标，而一个县轮廓文件将包括每个县的定位和形状。同样，一次旅行的 GPS 记录将包括旅行者在时间上的位置，描绘出他们在旅行中走过的路径。

重要的是要认识到，地理空间数据不仅包括地理空间信息本身。例如，以下轮廓本身并不特别有用：

然而，一旦你添加了适当的元数据，这些轮廓就更有意义了：

因此，地理空间数据包括每个描述项的时空信息（位置和形状）和非时空信息（元数据）。

空间信息通常表示为一系列坐标，例如：

location = (-38.136734, 176.252300)
outline = ((-61.686,17.024),(-61.738,16.989),(-61.829,16.996) ...)

这些数字对你来说可能没有太多意义，但一旦你将这些坐标系列绘制到地图上，数据突然变得容易理解：

地理空间数据有两种基本类型：

栅格数据：这是一种将世界划分为单元格并将值与每个单元格关联的地理空间数据。这与位图图像将图像划分为像素并将颜色与每个像素关联的方式非常相似；例如：

每个单元格的值可能代表在地图上绘制栅格数据时使用的颜色——这通常是为了提供一个栅格底图，其他数据可以在其上绘制——或者它可能代表其他信息，例如海拔、湿度水平或土壤类型。
矢量数据：这是一种由一系列特征组成的地理空间数据。例如，包含国家的 shapefile 将有一个特征对应于每个国家。对于每个特征，地理空间数据集将有一个几何形状，这是与该特征关联的形状，以及任何数量的属性，包含该特征的元数据。

一个特征的几何形状只是位于地球表面上的几何形状。这个几何形状由点、线（有时称为线字符串）和多边形或这些三种基本类型的组合构成：

你可能会遇到的典型栅格数据格式包括：

GeoTIFF 文件，基本上是带有地理参考信息的 TIFF 格式图像文件，用于在地球表面上准确定位图像。
美国地质调查局（USGS）的.dem文件，它以简单的 ASCII 数据格式存储数字高程模型（DEM）。
.png、.bmp和.jpeg格式的图像文件，以及相关的地理参考文件，用于在地球表面上定位图像。

对于矢量格式数据，你可能会遇到以下格式：

Shapefile：这是一种极其常见的文件格式，用于存储和共享地理空间数据。
WKT（已知文本）：这是一种基于文本的格式，常用于将几何形状从一个库或数据源转换为另一个。这也是从数据库检索特征时常用的格式。
WKB（Well-Known Binary）：这是 WKT 格式的二进制等效格式，将几何形状存储为原始二进制数据而不是文本。
GML（Geometry Markup Language）：这是一个基于 XML 的行业标准格式，通常在与其他 Web 服务通信时使用。
KML（Keyhole Markup Language）：这是由 Google 推广的另一种基于 XML 的格式。
GeoJSON：这是为存储和传输几何数据而设计的 JSON 版本。

因为你的分析只能与你要分析的数据一样好，所以获取和使用高质量的地理空间数据是至关重要的。确实，在执行地理空间分析时，一个主要挑战是获取适合工作的正确数据。幸运的是，有几个网站提供免费的高质量地理空间数据。但如果你在寻找更不为人知的数据集，你可能很难找到它。当然，你始终可以选择从头开始创建自己的数据，尽管这是一个极其耗时的工作过程。

我们将在第二章“地理空间数据”中回到地理空间数据的话题，我们将探讨什么使好的地理空间数据，以及如何获取它。

设置你的 Python 安装

要开始使用 Python 分析地理空间数据，我们将利用两个免费可用的第三方库：

GDAL：地理空间数据抽象库使你能够轻松以矢量和栅格格式读取和写入地理空间数据。
Shapely：正如其名所示，这是一个非常棒的库，它使你能够对几何形状执行各种计算。它还允许你操作形状，例如，通过将形状连接在一起或将它们拆分成各自的组成部分。

让我们继续安装这两个库到你的 Python 设置中，这样我们就可以立即开始使用它们了。

安装 GDAL

GDAL，或者更准确地说，GDAL/OGR 库，是由开源地理空间基金会发起的一个项目，旨在提供用于以各种格式读取和写入地理空间数据的库。从历史上看，GDAL 这个名字指的是用于读取和写入栅格格式数据的库，而 OGR 指的是用于访问矢量格式数据的库。这两个库现在已经合并，尽管名字仍然在类和函数名中使用，因此理解两者之间的区别是很重要的。

默认安装的 GDAL/OGR 允许你以 100 种不同的格式读取栅格地理空间数据，并以 71 种不同的格式写入栅格数据。对于矢量数据，GDAL/OGR 允许你以 42 种不同的格式读取数据，并以 39 种不同的格式写入。这使得 GDAL/OGR 成为访问和操作地理空间数据的一个极其有用的工具。

GDAL/OGR 是一个 C++库，具有各种绑定，允许您从其他语言访问它。在您的计算机上安装它后，您通常使用 Python 绑定通过 Python 解释器访问库。以下图解说明了这些各个部分是如何结合在一起的：

让我们继续安装 GDAL/OGR 库。GDAL（和 OGR）的主网站可以在gdal.org找到。

您如何安装它取决于您的计算机正在使用的操作系统：

对于 Windows 机器，您可以使用 FWTools 安装程序安装 GDAL/OGR，该安装程序可以从fwtools.maptools.org下载。

或者，您可以使用 OSGeo 安装程序安装 GDAL/OGR 和 Shapely，该安装程序可以在trac.osgeo.org/osgeo4w找到。
对于 Mac OS X，您可以从www.kyngchaos.com/software/frameworks下载 GDAL 和 OGR 的完整安装程序。
对于 Linux，您可以从 GDAL 的主网站下载 GDAL/OGR 的源代码，并按照网站上的说明从源代码构建它。您可能还需要安装 GDAL 和 OGR 的 Python 绑定。

安装完成后，您可以通过启动 Python 解释器并输入import osgeo.gdal然后import osgeo.ogr来检查它是否正常工作。如果每次 Python 命令提示符都重新出现而没有错误消息，那么 GDAL 和 OGR 已成功安装，您就可以开始使用了：

>>>import osgeo.gdal
>>>import osgeo.ogr
>>>

安装 Shapely

Shapely 是一个几何操作和分析库。它基于几何引擎，开源（GEOS）库，该库在 C++中实现了广泛的地理空间数据处理。Shapely 提供了一个 Python 接口来访问 GEOS，使得您可以直接在 Python 程序中使用这些操作。以下插图显示了您的 Python 代码、Python 解释器、Shapely 和 GEOS 库之间的关系：

Shapely 的主网站可以在pypi.python.org/pypi/Shapely找到。

网站上有您所需的一切，包括如何使用库的完整文档。请注意，要安装 Shapely，您需要下载 Shapely Python 包和底层 GEOS 库。GEOS 库的网站可以在trac.osgeo.org/geos找到。

您如何安装 Shapely 取决于您的计算机正在使用的操作系统：

对于 Windows，您应该使用 Shapely 网站上可用的预构建安装程序之一。这些安装程序包括它们自己的 GEOS 副本，因此不需要安装其他内容。
对于 Mac OS X，您应该使用在www.kyngchaos.com/software/frameworks可用的预构建 GEOS 框架。

小贴士

注意，如果你从前面的网站安装了 GDAL Complete 软件包，那么你电脑上已经安装了 GEOS。

一旦安装了 GEOS，你可以使用 pip，Python 包管理器来安装 Shapely：
```
pip install shapely
```
如果你电脑上没有安装 pip，你可以按照 pip.pypa.io/en/latest/installing.html 中的说明进行安装。
对于 Linux 机器，你可以从 GEOS 网站下载源代码并自行编译，或者安装一个包含 GEOS 的合适 RPM 或 APT 软件包。一旦完成，你可以使用 pip install shapely 来安装 Shapely 库本身。

安装完成后，你可以通过运行 Python 命令提示符并输入以下命令来检查 Shapely 库是否正常工作：

>>> import shapely.geos
>>>

如果你再次得到 Python 命令提示符而没有错误，就像前面的例子一样，那么 Shapely 已经成功安装，你可以开始使用了。

获取一些地理空间数据

对于本章，我们将使用一个简单但仍然非常有用的地理空间数据文件，称为 世界边界数据集。该数据集由单个 shapefile 组成，其中 shapefile 中的每个特征代表一个国家。对于每个国家，相关的几何对象代表该国的轮廓。附加属性包含诸如国家名称、ISO 3166-1 代码、总面积、人口和联合国区域分类等元数据。

要获取世界边界数据集，请访问 thematicmapping.org/downloads/world_borders.php。

滚动到下载部分，并点击文件以下载。确保你下载的是完整版本而不是简化版本——你想要的文件将被称为 TM_WORLD_BORDERS-0.3.zip。

注意，shapefile 以 ZIP 归档的形式提供。这是因为 shapefile 由多个文件组成，如果它们存储在 ZIP 归档中，则更容易分发。下载文件后，双击 ZIP 归档以解压缩它。你将得到一个名为 TM_WORLD_BORDERS-0.3 的目录。在这个目录中应该包含以下文件：

以下表格解释了这些不同的文件以及它们包含的信息：

文件名	描述
`Readme.txt`	这是你的典型 README 文件，包含有关 shapefile 的有用信息。
`TM_WORLD_BORDERS-0.3.shp`	该文件包含每个特征的几何数据。
`TM_WORLD_BORDERS-0.3.shx`	这是一个 `.shp` 文件的索引，使得可以快速访问给定特征的几何形状。
`TM_WORLD_BORDERS-0.3.dbf`	这是一个数据库文件，包含每个特征的各个属性。
`TM_WORLD_BORDERS-0.3.prj`	此文件以纯文本文件的形式描述了数据使用的坐标系统和投影。

将此目录放置在方便的位置。我们将在整本书中广泛使用这个数据集，因此您可能希望在某个地方保留一个备份副本。

解锁 shapefile

最后，我们准备好开始处理一些地理空间数据。打开一个命令行或终端窗口，并将 cd 命令切换到您之前解压缩的 TM_WORLD_BORDERS-0.3 目录。然后输入 python 来启动您的 Python 解释器。

我们将首先加载我们之前安装的 OGR 库：

>>> import osgeo.ogr

接下来，我们想使用 OGR 打开 shapefile：

>>> shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")

执行此语句后，shapefile 变量将包含一个 osgeo.ogr.Datasource 对象，它代表我们打开的地理空间数据源。OGR 数据源可以支持多个信息层，尽管 shapefile 只有一个层。因此，我们接下来需要从 shapefile 中提取（唯一的一个）层：

>>>layer = shapefile.GetLayer(0)

让我们遍历 shapefile 中的各种特征，依次处理每个特征。我们可以使用以下方法：

>>> for i in range(layer.GetFeatureCount()):
>>>     feature = layer.GetFeature(i)

feature 对象，osgeo.ogr.Feature 的一个实例，允许我们访问与特征相关的几何形状以及特征的属性。根据 README.txt 文件，国家的名称存储在一个名为 NAME 的属性中。现在让我们提取这个名称：

>>>    feature_name = feature.GetField("NAME")

注意

注意，属性是大写的。Shapefile 属性是区分大小写的，因此你必须使用确切的字母大小写才能获取正确的属性。使用 feature.getField("name") 会生成错误。

要获取特征几何对象的引用，我们使用 GetGeometryRef() 方法：

>>>     geometry = feature.GetGeometryRef()

我们可以用几何形状做很多事情，但就目前而言，让我们看看我们有什么类型的几何形状。我们可以使用 GetGeometryName() 方法来做到这一点：

>>>>    geometry_type = geometry.GetGeometryName()

最后，让我们打印出我们为此特征提取的信息：

>>>    print i, feature_name, geometry_type

这里是我们编写的完整小程序，用于解锁 shapefile 的内容：

import osgeo.ogr
shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    feature_name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    geometry_type = geometry.GetGeometryName()
    print i, feature_name, geometry_type

如果您第二次按 Return 键来关闭 for 循环，程序将运行，显示从 shapefile 中提取的每个国家的有用信息：

0 Antigua and Barbuda MULTIPOLYGON
1 Algeria POLYGON
2 Azerbaijan MULTIPOLYGON
3 Albania POLYGON
4 Armenia MULTIPOLYGON
5 Angola MULTIPOLYGON
6 American Samoa MULTIPOLYGON
7 Argentina MULTIPOLYGON
8 Australia MULTIPOLYGON
9 Bahrain MULTIPOLYGON
...

注意，一些国家的几何形状是多边形，而其他国家的几何形状是多边形集合。正如其名称所暗示的，多边形集合只是一系列多边形。因为几何形状代表每个国家的轮廓，所以当国家的轮廓可以用单个形状表示时，使用多边形；而当轮廓有多个部分时，使用多边形集合。这种情况最常见于由多个岛屿组成的国家。例如：

如您所见，阿尔及利亚由一个多边形表示，而澳大利亚及其外围岛屿则是一个多边形。

分析数据

在上一节中，我们得到了一个表示每个国家轮廓的osgeo.ogr.Geometry对象。虽然我们可以直接对这个几何对象做很多事情，但在这个例子中，我们将轮廓复制到 Shapely 中，以便我们可以利用 Shapely 的地理空间分析功能。为此，我们必须将几何对象从 OGR 导出，并将其作为 Shapely 对象导入。为此，我们将使用 WKT 格式。仍然在 Python 解释器中，让我们获取单个特征的几何形状并将其转换为 Shapely 对象：

>>> import shapely.wkt
>>> feature = layer.GetFeature(0)
>>> geometry = feature.GetGeometryRef()
>>> wkt = geometry.ExportToWkt()
>>> outline = shapely.wkt.loads(wkt)

因为加载了特征编号0，我们检索到了安提瓜和巴布达的轮廓，如果我们在地图上显示它，看起来会像以下这样：

outline变量以 Shapely MultiPolygon对象的形式保存了这个国家的轮廓。现在我们可以使用这个对象来分析几何形状。以下是一些我们可以用 Shapely 几何形状做的有用事情：

我们可以计算质心，这是几何形状中中心位置的点。
我们可以计算几何形状的边界框。这是一个定义多边形北部、南部、东部和西部边缘的矩形。
我们可以计算两个几何形状之间的交集。
我们可以计算两个几何形状之间的差值。

注意

我们还可以计算诸如每个多边形的长度和面积等值。然而，因为 World Borders Dataset 使用所谓的未投影坐标，所以得到的长度和面积值将以度为单位来衡量，而不是米或英里。这意味着计算出的长度和面积可能不太有用。我们将在下一章探讨地图投影的性质，并找到一种方法来解决这个问题，这样我们就可以为多边形计算有意义的长度和面积值。但这对于我们来说现在太复杂了。

让我们显示我们特征质心的纬度和经度：

>>> print outline.centroid.x, outline.centroid.y
-61.791127517 17.2801365868

因为 Shapely 不知道多边形位于哪个坐标系中，它使用更通用的x和y属性来表示一个点，而不是讨论经纬度值。记住，纬度对应南北方向的位置，即y值，而经度是东西方向的位置，即x值。

我们也可以显示大纲的边界框：

>>> print outline.bounds
(-61.891113, 16.989719, -61.666389, 17.724998)

在这种情况下，返回的值是最小经纬度和最大经纬度（即min_x、min_y、max_x、max_y）。

当然，我们可以用 Shapely 做更多的事情，但这足以证明 Shapely 库正在工作，并且我们可以从 shapefile 中读取地理空间数据并将其转换为 Shapely 几何对象进行分析。

这是我们使用 Python shell 直接进行的极限——shell 对于这种快速实验来说很棒，但当你打字错误时，不得不重新输入行（或使用命令历史）会很快变得令人厌烦。对于任何更严肃的事情，您将想要编写一个 Python 程序。在本章的最后部分，我们将做 exactly that：创建一个基于我们所学知识来解决有用地理空间分析问题的 Python 程序。

一个识别邻国程序

对于我们的第一个真正的地理空间分析程序，我们将编写一个 Python 脚本，用于识别邻国。基本概念是提取每个国家的多边形或复合多边形，并查看哪些其他国家与每个多边形或复合多边形接触。对于每个国家，我们将显示一个与该国接壤的其他国家的列表。

让我们从创建 Python 脚本开始。创建一个名为borderingCountries.py的新文件，并将其放置在您之前下载的TM_WORLD_BORDERS-0.3.shp形状文件所在的同一目录中。然后输入以下内容到这个文件中：

import osgeo.ogr
import shapely.wkt

def main():
    shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
    layer = shapefile.GetLayer(0)

    countries = {} # Maps country name to Shapely geometry.

    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        country = feature.GetField("NAME")
        outline = shapely.wkt.loads(feature.GetGeometryRef().ExportToWkt())

        countries[country] = outline

    print "Loaded %d countries" % len(countries)

if __name__ == "__main__":
    main()

到目前为止，这相当直接。我们正在使用我们之前学到的技术将形状文件的内容读入内存，并将每个国家的几何形状转换为 Shapely 对象。结果存储在countries字典中。最后，请注意，我们将程序逻辑放入了一个名为main()的函数中——这是一个好的实践，因为它允许我们使用return语句来处理错误。

现在运行您的程序以确保它正常工作：

$ python borderingCountries.py
Loaded 246 countries

我们接下来的任务是识别接壤国家。我们的基本逻辑是遍历每个国家，然后找到与这个国家接壤的其他国家。以下是相关代码，您应该将其添加到main()函数的末尾：

    for country in sorted(countries.keys()):
        outline = countries[country]

        for other_country in sorted(countries.keys()):

            if country == other_country: continue

            other_outline = countries[other_country]

            if outline.touches(other_outline):

                print "%s borders %s" % (country, other_country)

如您所见，我们使用touches()方法来检查两个国家的几何形状是否接触。

运行此程序现在将显示相邻的国家：

$ python borderingCountries.py
Loaded 246 countries
Afghanistan borders Tajikistan
Afghanistan borders Uzbekistan
Albania borders Montenegro
Albania borders Serbia
Albania borders The former Yugoslav Republic of Macedonia
Algeria borders Libyan Arab Jamahiriya
Algeria borders Mali
Algeria borders Morocco
Algeria borders Niger
Algeria borders Western Sahara
Angola borders Democratic Republic of the Congo
Argentina borders Bolivia
...

恭喜！您已经编写了一个简单的 Python 程序来分析国家轮廓。当然，还有很多可以改进和扩展这个程序的地方。例如：

您可以添加命令行参数，让用户指定形状文件的名称以及用于显示国家名称的属性。
您可以添加错误检查来处理无效和非存在的形状文件。
您可以添加错误检查来处理无效的几何形状。
您可以使用空间数据库来加快处理过程。程序目前需要大约一分钟来完成，但使用空间数据库将大大加快这个过程。如果您处理的是大量空间数据，正确索引的数据库绝对是必不可少的，否则您的程序可能需要几周时间才能运行。

我们将在本书的后面部分讨论所有这些内容。

摘要

在本章中，我们通过探讨你通常需要解决的问题类型以及你将要处理的数据类型，开始了对地理空间分析的探索。我们发现并安装了两个主要的 Python 库来处理地理空间数据：GDAL/OGR 用于读取（和写入）数据，以及 Shapely 用于执行地理空间分析和操作。然后我们下载了一个包含国家数据的简单但有用的 shapefile，并学习了如何使用 OGR 库读取该 shapefile 的内容。

接着，我们看到了如何将 OGR 几何对象转换为 Shapely 几何对象，然后使用 Shapely 库来分析和操作该几何对象。最后，我们创建了一个简单的 Python 程序，结合了我们所学的一切，将国家数据加载到内存中，然后使用 Shapely 找到相互接壤的国家。

在下一章中，我们将更深入地探讨地理空间数据的话题，了解更多的地理空间数据类型和概念，以及探索一些主要的免费地理空间数据来源。我们还将学习为什么拥有良好的数据来工作很重要——以及如果你没有良好的数据会发生什么。

第二章：地理空间数据

在本章中，我们将重点关注用于地理空间分析的数据。你将了解更多关于地理空间数据性质的知识，并发现一些你可以免费获取地理空间数据集的主要网站。然后我们将探讨使用 Python 读取和写入地理空间数据的方法。

尤其是本章将涵盖以下主题：

为什么高质量的地理空间数据很重要
你可能会遇到的地理空间数据的各种类型
可免费获取的地理空间数据集的主要来源
如何使用 GDAL/OGR 库读取和写入地理空间数据
如何处理空间参考系统
地理空间数据错误及其修复方法

让我们先看看为什么拥有高质量的地理空间数据很重要。

地理空间数据质量

想象一下，你正在编写一个程序，需要将每个城市和镇的位置显示在栅格底图上。你尽职尽责地获取了一个用于底图的优秀栅格数据源，然后在网上搜索城市和镇数据的来源。你选择了国家地理空间情报服务（NGIS）网站下载地名数据库，然后将其绘制到你的地图上。该数据库包括许多其他信息，例如每个地名的纬度和经度：

位置	纬度	经度
Abache	7.3551	7.6407
Abacheke	5.50372	6.729519
Abacher	13.816667	20.816667
Abacheri	14.183333	41.5
Abachi	7.3551	7.6407
...等等

到目前为止一切顺利，但当你完成程序后，当用户放大你的地图时，位置看起来可疑地规则：

如果你在地图上绘制一个网格，你可以清楚地看到问题所在：

如你所见，位置是均匀分布的——尽管纬度和经度值有很多精度，但实际上它们只精确到大约两位小数。在显示荷兰部分地区的前一个图像中，这可能会导致位置偏差近一公里。

这只是如果你不使用高质量的地理空间数据可能会出错的一种类型的事例。另一个例子通常出现在执行点在多边形内计算时——即，在尝试确定一个给定点是在给定多边形内部还是外部时：

如果多边形代表的是，例如，一个国家的轮廓，那么你可以使用点在多边形内计算来确定给定位置是否在国家的边界内。这通常用于地理定位一个给定的点（即，将一个点与一个或多个已知位置匹配）。

现在，当你尝试定位接近多边形边缘的点时，如果你的多边形不够详细，你很容易出现地理定位错误：

在前面的图像中，点代表一个需要地理定位的点。这个点在旧金山是一个有效的位置，但由于 polygon 不够详细，这个点位于旧金山市 polygon 之外，因此地理定位失败。

当使用覆盖较大区域的 polygon，如州或国家轮廓时，这个问题尤为严重。由于这些情况中 polygon 的大小，精度通常会被牺牲。

你可能会认为这个问题的答案是拥有更详细的 polygon——也就是说，使用更多的点来构成 polygon 的轮廓，以便更准确地表示所需的轮廓（在这种情况下，旧金山海岸线）。然而，更多的细节并不总是更好的。一个 polygon 越详细，处理它所需的时间就越长——如果点太多，你的程序可能会因为处理过多的数据而崩溃。

解决这个问题有方法；例如，可以通过将大型多边形分割成更小的部分，或者通过缓冲区多边形以包括靠近边缘的点。但重要的是要认识到，高质量的数据并不总是意味着高精度数据；它意味着数据对于你想要使用的目的来说是合适的。

我们现在将继续通过查看你可能会遇到的多种地理空间数据类型来继续我们对地理空间概念的探索。

地理空间数据类型

在前一章中，我们简要地查看了一些用于存储和传输栅格和矢量地理空间数据的常见格式。现在，让我们看看你可能会遇到的一些重要的地理空间数据类型。

Shapefiles

如前一章所述，shapefile 是磁盘上的一组文件集合，这些文件一起存储了一组地理空间特征及其属性和几何形状。例如，以下插图显示了典型 shapefile 中存储的数据：

由于 shapefile 格式已经存在多年，并且可以追溯到 dBase 时代，一个 shapefile 由几个单独的文件组成。通常，这些文件被组合成一个 ZIP 存档以进行分发。

Shapefiles 之所以非常受欢迎，是因为它们使得存储和分发地理空间数据变得非常容易。实际上，几乎所有处理地理空间数据的 GIS 系统和库都能够理解 shapefile 格式。

然而，shapefiles 确实有一些缺点：

与几乎所有的其他地理空间数据格式不同，shapefile 中属性的名称是区分大小写的。当你的代码在另一种格式（例如，在数据库中）中处理数据时，这可能会引起问题，但当你尝试访问 shapefile 中的属性时，代码突然停止。
单个 shapefile 中存储的所有几何形状必须是同一类型。这意味着，例如，您不能有些特征由线表示，而其他特征由点表示。更严重的是，当处理复合几何形状（如 MultiLines 或 MultiPolygons）或尝试在 shapefile 中存储几何集合时，这会导致问题。
虽然 shapefiles 允许您按顺序读取特征，但通常没有空间索引支持。这意味着您不能根据地球表面上特征的位置进行搜索。例如，回答“哪些国家在伦敦 1,000 公里范围内？”这样的问题，您需要逐个检查每个特征，这并不特别高效。

已知文本

知名文本（WKT）格式通常不用于存储地理空间数据。相反，它用于将几何形状从一个格式传输到另一个格式。我们在上一章中看到了一个例子，我们使用 OGR 库提取了一个几何形状，然后将其转换为 WKT，以便我们可以将其重新创建为一个 Shapely 几何对象。

WKT 是一种非常紧凑且易于使用的格式。例如，以下 WKT 字符串定义了纽约市中央公园中间的一个点几何形状：

POINT(-73.967344 40.782148)

如您所见，点被表示为一个x（经度）值，一个空格，然后是y（纬度）值。相同的通用格式用于表示多边形的坐标。例如，以下是一个 WKT 格式的多边形，这次定义了中央公园的大致轮廓：

POLYGON((-73.973057 40.764356, -73.981898 40.768094, -73.958209 40.800621, -73.949282 40.796853, -73.973057 40.764356))

除了使用 WKT 在不同系统和库之间传输数据外，当您需要快速将几何形状硬编码到 Python 代码中时，您也会发现 WKT 字符串很有用。例如，以下代码展示了如何快速创建一个 Shapely 多边形进行测试：

p = shapely.wkt.loads("POLYGON((23.4 38.9, 23.5 38.9, 23.5 38.8, 23.4 38.9))")

如果您想将几何数据存储在文本文件中，例如将分析输出临时保存到磁盘上以便将其加载到另一个程序进行进一步处理，WKT 格式也非常有用。

已知二进制

知名的二进制（WKB）格式是 WKT 的二进制等价物。WKB 通常仅用于在数据库内部传输和存储地理空间数据。在 WKB 格式中，坐标以双精度浮点数存储，并使用数字代码来表示几何类型。与 WKT 相比，这种格式对计算机的读写速度更快，尽管当然这种格式对人类来说不容易理解。

空间数据库

正如普通数据库用于存储大量数据并允许用户对数据进行高效查询一样，空间数据库是一种设计用于存储几何形状并根据每个几何形状在空间中的位置执行高效查询的数据库。例如，您可以快速找到在给定点 20 英里范围内的所有道路交通事故，或者找到您当前位置最近的岛屿。

设置和使用空间数据库是一项相当复杂的任务。有几种不同类型的空间数据库可供选择。在免费选项中，轻量级的 SpatiaLite 数据库和功能强大但复杂的 PostGIS 数据库是最受欢迎的选择。

创建数据库后，您必须配置数据库以处理空间数据。您还需要使用数据库自己的语法来存储和查询您的空间数据——您如何做这取决于您使用哪种数据库。

无论您使用哪种类型的数据库，尝试检索空间几何体通常将返回 WKT 格式的字符串或 WKB 格式的原始二进制数据。然后您可以将这些转换为不同的格式（例如，Shapely 几何对象）进行处理。

当然，在空间数据库中可用的强大功能，尤其是 PostGIS 中内置的空间操作和查询功能，您可能不需要在数据库本身之外进行任何空间分析。例如，一旦您要求数据库识别给定半径内的所有道路交通事故，您通常会检索这些事故的非空间信息。所有的空间处理都是在数据库内部完成的，一旦您找到了所需的记录集，您就可以像在非空间数据库中一样检索和处理结果。

我们将在第三章 空间数据库 中回到空间数据库的主题。

地理空间微格式

所说的 地理空间微格式通常由 API 用于发送和接收地理空间数据。各种公司和组织已经定义了自己的标准来传输地理空间数据，因此随着时间的推移，已经开发出多种常见格式。如果您使用使用这些微格式之一的 API，您需要熟悉这些数据格式。

我们在这里将探讨两种地理空间微格式：GeoJSON 和 GML。

GeoJSON

GeoJSON (geojson.org) 是一种开放标准，用于将地理空间数据结构表示为 JavaScript 对象表示法（JSON）对象。例如，以下 GeoJSON 格式的字符串用于表示一个点几何体：

{"type": "Point", "coordinates": [-73.967344, 40.782148]}

由于 GeoJSON 建立在 JSON 标准之上，因此你可以使用json标准 Python 库在 GeoJSON 格式字符串和 Python 字典之间进行转换。这使得在 Python 程序中使用 GeoJSON 变得特别容易。

GeoJSON 标准包括对以下内容的支持：

将任何标准几何对象（点、LineString、Polygon、MultiPoint、MultiLineString、MultiPolygon 和 GeometryCollection）表示为 GeoJSON 字符串。
将一个特征作为 GeoJSON 字符串存储，包括特征的几何形状、任意数量的属性（在 GeoJSON 中称为属性）以及可选的空间参考系统。
使用现有的 XML 模式和处理工具。因为 GML 格式基于 XML，所以可以使用现有的 XML 解析器和验证器处理 GML 数据。你还可以创建一个特定于应用程序的 XML 模式，定义你对 GML 标准的扩展，然后使用现有的 XML 库与该模式一起工作。

GeoJSON 格式被处理地理空间数据的软件广泛支持。实际上，GDAL/OGR 库包括对读取和写入 GeoJSON 格式数据的支持，同样，我们稍后将要使用的 Mapnik 库也支持。

GML

地理标记语言（GML）是一种基于 XML 的格式，用于以文本形式存储几何和特征。GML 是一个复杂且精细的标准。因为它基于 XML，所以 GML 格式的数据往往相当冗长。例如，以下 GML 字符串表示一个最小的点几何：

<gml:Point>
    <gml:pos>40.782148 -73.967344</gml:pos>
</gml:Point>

GML 包括以下支持：

表示点、LineString 和多边形几何。

注意

GML 标准的 3.0 版本增加了对栅格格式数据的支持。
定义特征，并为每个特征存储属性。
将多个几何与每个特征关联；例如，一个特征可能有一个轮廓、一个边界框和一个质心，所有这些都定义为与特征关联的几何。
定义几何所使用的空间参考系统。
允许你为特定数据集使用 GML 标准子集的配置文件；例如，GML 简单特征配置文件将数据限制在表示几何及其相关属性。
由于 GML 格式基于 XML，因此可以使用现有的 XML 模式和处理工具；例如，你可以定义一个特定于应用程序的 XML 模式，定义你对 GML 标准的扩展，然后使用 XML 解析器和验证器处理你的 GML 数据。

GML 标准是由开放地理空间联盟（OGC）开发的，现在已被接受为 ISO 标准。GML 被 OGC 定义的各个网络标准广泛使用，并且当你想要访问遵循这些标准之一的 API 时，你将使用 GML，例如 Web 特征服务。

数字高程模型

数字高程模型（DEM）是一种以栅格格式数据表示地球表面曲线和等高线的迷人方式。正如我们在上一章中提到的，栅格格式的地理空间数据将世界划分为单元格，并将信息与每个单元格关联。在 DEM 中，每个单元格包含关于该点地球表面高度的信息。例如，考虑以下来自典型 DEM 文件的高程数据：

2874  2871  2874  2933  2995  3022  3028  3031  3035  3031
2874  2871  2874  2933  2992  3012  3025  3028  3031  3028
2871  2871  2877  2932  2989  3007  3018  3025  3023  3020
2872  2871  2886  2935  2975  2997  3010  3020  3022  3023
2871  2879  2903  2942  2965  2991  3005  3015  3022  3026
2871  2887  2930  2972  2992  2998  3013  3023  3029  3031
2880  2899  2941  2992  3005  3005  3021  3028  3033  3039
2896  2920  2956  3000  3013  3019  3019  3028  3037  3042
2915  2939  2981  3008  3017  3026  3028  3028  3036  3044
2928  2952  2986  3024  3029  3034  3038  3034  3031  3044
2936  2960  3009  3040  3044  3046  3049  3044  3037  3044
2943  2977  3041  3051  3051  3051  3051  3051  3037  3046
2960  3029  3051  3051  3051  3051  3051  3050  3044  3049

这份数据是从俄勒冈州 Forked Horn Butte 的 DEM 文件中提取的。每个数字表示海平面以上的高度，单位为英尺。如果将这些高程值绘制成三维图，就可以揭示地球表面的形状，如下所示：

这只是整体 DEM 文件的一小部分，但它确实展示了 DEM 文件是如何编码地球表面形状的。

注意

DEM 文件还有一个称为无数据值的概念。这是一个特殊的高度值，表示该点没有高程值。无数据值用于你不想记录或显示 DEM 中某些部分的高程值的地方。例如，特定国家的 DEM 文件会使用无数据值来表示该国边界以外的所有地区。

数字高程模型通常用作构建地球表面有用图像的基石。例如：

可以使用称为颜色映射的技术将不同的颜色与不同的高度带相关联。如果选择了正确的颜色组合，结果几乎可以看起来像一张显示不同植被、裸地、岩石和雪的带的照片。
可以从高程数据生成阴影渲染图像。这模仿了光源（例如太阳）照射到地球表面的效果，揭示了深度并产生阴影和亮部，使得生成的图像几乎像是从太空中拍摄的地球照片。
等高线可以通过平滑 DEM 数据并通过 GDAL 库提供的程序gdal_contour生成。

通常，将几个生成的图像合并以产生更逼真的效果。然后，这些派生图像被用作背景地图，在背景地图上叠加地理空间数据。

栅格底图

而不是从 DEM 文件构建你的图像，你可以使用为底图预先生成的图像。这些底图通常非常复杂。例如，水下区域可能使用不同深度的蓝色色图来绘制，以表示水的深度，而海平面以上的区域则使用阴影渲染图像，结合植被和基于高度的颜色来产生逼真的效果。

下图展示了这种类型的典型底图：

这些基础地图仅仅是带有相关地理参照信息的图像文件。地理参照信息标识了基础地图覆盖的地球区域。这通常是通过指定图像的左上角和右下角的纬度和经度来完成的。使用这些点，可以在地球表面上准确放置图像，允许在基础地图上正确绘制地理空间数据，并且根据你希望显示的地球表面区域，允许正确绘制基础地图的相应部分。

多波段栅格文件

如前一章所述，栅格格式的地理空间数据不仅能存储图像。栅格信息可能包括诸如高程（如我们在数字高程模型部分之前所见）、土壤类型、平均降雨量、人口密度等值。

栅格格式的数据文件不仅限于存储单一信息。单个文件可以包含多个波段的栅格数据，如下面的插图所示：

每个波段都有一个单元格的值，因此对于给定的（x，y）位置，波段 1 将包含一个值，波段 2 将包含一个值，依此类推。存储在每个波段中的值的含义完全取决于你使用的栅格文件；你需要参考该栅格文件的文档以了解每个波段中存储了什么。

有时，多个波段可以组合起来产生一种颜色。例如，你可以下载由陆地卫星捕获的栅格数据（有关详细信息，请参阅landsatlook.usgs.gov），这些数据包括三个单独波段中的红色、绿色和蓝色颜色成分。额外的波段包含红外和全色值，在某些情况下也可能很有用。

可免费获取的地理空间数据来源

现在你已经了解了拥有适当地理空间数据的重要性，并且已经了解了你将想要使用的主要数据类型，让我们看看你可以获取所需数据的几个地方。

有一些情况下你可能需要购买地理空间数据集。一个例子是在寻找美国邮政编码边界时；此类信息属于美国邮政服务（USPS）的专有信息，并且只能通过从获得 USPS 许可以销售此类数据的供应商处购买合适的数据集来获取准确版本。然而，这只是一个例外：在几乎所有其他情况下，你可以免费获取、修改和使用地理空间数据。

现在我们来看看在寻找地理空间数据时你将想要使用的几个主要网站。

自然地球数据

自然地球数据网站 (naturalearthdata.com) 是高质量和免费地理空间数据的综合来源。在矢量格式数据方面，文件以高、中、低分辨率提供。提供了两种不同类型的矢量数据：

文化数据：这包括国家、州或省、城市区域和公园轮廓的多边形，以及人口密集地区、道路和铁路的点数据和 LineString 数据。
物理数据：这包括陆地、海岸线、海洋、小岛屿、珊瑚礁、河流和湖泊的多边形和 LineStrings。

在栅格格式数据方面，自然地球数据在 1:100 万和 1:500 万比例尺下提供了五种不同类型的栅格基础地图。

这些图像文件以地理参考的 TIFF 图像提供，这使得它们很容易在您的 Python 程序中用作栅格基础地图。

OpenStreetMap

OpenStreetMap (openstreetmap.org) 是一个巨大的协作努力，旨在创建和提供地理空间地图数据。该网站将其描述为“由像您这样的人制作的整个世界的免费可编辑地图”。它将自己定位为 Google Maps 的直接竞争对手。以下图像显示了基于 OpenStreetMap 数据的新西兰罗托鲁阿市的部分街道地图：

不幸的是，OpenStreetMap 使用其自有的基于 XML 的格式来存储地理空间数据。如果您愿意，您可以下载整个 OpenStreetMap 数据库，称为 Planet.osm，然后使用空间数据库来访问这些信息。然而，在大多数情况下，您可能希望使用已经转换为更标准格式（如 shapefile）的 OpenStreetMap 数据库的提取。

您可以在 wiki.openstreetmap.org/wiki/Planet.osm#Country_and_area_extracts 找到提供 OpenStreetMap 数据提取的网站列表。

如果您想要操作街道地图或使用街道地图作为背景来显示其他地理空间数据，OpenStreetMap 非常有用。

美国人口普查局

美国人口普查局已将大量地理空间数据以 TIGER（拓扑集成地理编码和参考系统）的名义提供。TIGER 数据集包括街道、铁路、河流、湖泊、地理边界以及州、学区、城市边界等法律和统计区域。

TIGER 数据以 shapefile 格式提供，可以从 www.census.gov/geo/maps-data/data/tiger.html 下载。

由于它是由美国政府生产的，TIGER 数据仅覆盖美国及其保护领地（即波多黎各、美属萨摩亚、北马里亚纳群岛、关岛和美国维尔京群岛）。然而，对于这些地区，TIGER 是准确地理空间数据的绝佳来源。

世界边界数据集

我们在上一章中使用了这个数据集。虽然它非常简单，但世界边界数据集(thematicmapping.org/downloads/world_borders.php)以 shapefile 的形式提供了有用的国家轮廓。shapefile 包括国家的名称、相关的 ISO、FIPS 和联合国识别代码、基于联合国的地区和子地区分类，以及国家的人口和面积。

世界边界数据集的简单性使其成为许多需要整个世界基本地图的地理空间应用的理想选择。

GLOBE

全球陆地一公里基础高程（GLOBE）数据集是一个国际努力，旨在生产高质量、中等分辨率的 DEM 数据，覆盖整个世界。在栅格 DEM 文件中的每个单元格代表地球表面上一个正方形的地面高程，该正方形的经度为 30 弧秒，纬度为 30 弧秒。这相当于每边大约一公里的正方形。

GLOBE 项目的官方网站可在www.ngdc.noaa.gov/mgg/topo/globe.html找到。请注意，如果您下载了覆盖地球表面区域的预制的“瓦片”，您还需要下载相关的头文件（.hdr），该头文件为 DEM 数据提供地理参考。这些头文件可以从www.ngdc.noaa.gov/mgg/topo/elev/esri/hdr下载。

注意

由于在本章后面我们需要一些示例 DEM 数据，现在就下载E瓦片。然后访问提供的链接下载相关的头文件。您应该最终得到两个文件，分别命名为e10g和e10g.hdr。

国家高程数据集

国家高程数据集(ned.usgs.gov)为美国大陆、阿拉斯加、夏威夷和其他美国领土提供了高分辨率 DEM 数据。根据您查看的区域，DEM 数据集中的每个单元格对应于 3 到 60 平方米的面积。这比 GLOBE 项目使用的 1 平方公里分辨率要高得多。

如果您想为美国制作自己的阴影高程基础地图，国家高程数据集是一个极佳的选择。所有文件都可用多种格式，包括 GeoTIFF 和 ArcGRID，这些都可以使用 GDAL 进行处理。

使用 Python 读取和写入地理空间数据

由于我们将使用 GDAL/OGR 库来访问地理空间数据，让我们更详细地看看如何使用这个库读取和写入矢量格式和栅格格式数据。

读取矢量数据

在上一章中，我们编写了一个简单的程序，用于从 shapefile 中读取要素。以下是该程序的副本：

import osgeo.ogr
shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    feature_name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    geometry_type = geometry.GetGeometryName()
    print i, feature_name, geometry_type

让我们更详细地看看这个程序是如何工作的，以及更一般地，如何使用 OGR 库读取矢量格式数据。

当读取地理空间数据时，osgeo.ogr.Open() 函数只接受一个参数：要打开的数据集的名称。OGR 库遍历它支持的所有不同的驱动程序，直到找到一个能够读取此数据集的驱动程序。然后，驱动程序创建一个新的 OGRDataSource 对象，该对象提供对数据集内容的访问，Open() 函数将此对象返回给调用者。

所有这些操作的效果是将 shapefile 变量设置为 OGR 数据源。OGR 数据源由一个或多个图层组成，每个图层代表一组不同的数据。例如，Countries 数据源可能有一个表示国家地形的图层，一个包含道路的图层，另一个表示国家城市边界的图层，另一个表示区域边界的图层，等等。

注意

请记住，shapefile 只能有一个图层。为了使用 shapefile 表示这些不同的信息，您必须为这些不同的数据分别创建一个单独的 shapefile。

如前述代码示例所示，您使用 GetLayer() 方法从数据源中检索一个图层；返回的对象是 OGRLayer 类的实例。还有一个方便的 GetLayerCount() 方法，它返回数据源中的图层数量。

每个图层都有一个 空间参考系统，它告诉您如何解释图层中的各个坐标，以及包含实际数据的要素列表。如果您不知道空间参考系统是什么，请不要担心；您将在本章后面的 处理空间参考系统 部分中了解到所有关于它的信息。

我们可以使用 GetFeatureCount() 和 GetFeature() 方法遍历图层中的各个要素。正如您所期望的，每个要素都由 ogr.Feature 类的实例表示。

每个要素都有一个唯一的 ID，可以使用 GetFID() 方法检索，以及一个 几何形状 和一个属性列表。我们使用 GetGeometryRef() 方法检索几何形状（OGRGeometry 的实例），并且可以使用 GetField() 方法访问要素的属性。

使用这些各种类和方法，您可以遍历矢量数据源中的各种功能，依次检索每个功能的几何形状和属性（以及 ID）。然而，所有这些的奇妙之处在于，这并不重要您的数据格式是什么：您使用完全相同的过程从 shapefile 中读取数据，就像您使用 OGR 库的 PostGIS 数据库驱动程序从空间数据库中读取数据一样。OGR 库隐藏了如何读取不同数据格式的所有细节，并为您提供了一个简单的、高级的接口来从任何数据源读取矢量格式数据。

写入矢量数据

将地理空间数据写入矢量格式文件几乎与读取它一样简单。然而，您必须采取一些额外的步骤。让我们编写一个简单的 Python 程序，该程序创建一个 shapefile 并将一些示例数据保存到其中。这个程序将教会您所有关于使用 OGR 编写矢量格式数据所需知道的事情。

使用 OGR 创建矢量格式数据集涉及以下步骤：

首先，我们通过选择一个 OGR 驱动程序并告诉该驱动程序创建新的数据源来创建目标文件：

from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
dstFile = driver.CreateDataSource("test-shapefile")

然后，我们创建一个 空间参考 对象，该对象定义了数据集中坐标应该如何解释：
```
from osgeo import osr
spatialReference = osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("WGS84")
```
正如您所看到的，我们使用 osr 模块来定义空间参考，然后使用代码 WGS84 将其设置为“已知”空间参考。

注意

WGS84 是用于经纬度值的标准。我们将在本章后面的 处理空间参考系统 部分详细讨论这一点。
然后，我们将图层添加到目标文件中，以存储图层的数据：
```
layer = dstFile.CreateLayer("layer", spatialReference)
```
正如您所看到的，每个图层都有自己的空间参考，因此我们必须在创建图层时传递我们之前定义的空间参考。
下一步是定义目标文件将为每个功能存储的各种属性。让我们定义一个名为 NAME 的字段：
```
field = ogr.FieldDefn("NAME", ogr.OFTString)
field.setWidth(100)
layer.CreateField(field)
```
注意

注意，我们使用大写字母定义字段名称。这是因为我们正在写入 shapefile，它对属性名称的大小写敏感，通常使用大写字母定义属性。在 shapefile 中使用大写属性名称将有助于避免以后的问题。

这完成了文件的创建。现在让我们创建一些示例数据并将其保存到文件中。这涉及以下步骤：

让我们定义一个简单的多边形来表示功能的几何形状。我们将使用 WKT 格式来简化这个过程：
```
wkt = "POLYGON((23.4 38.9, 23.5 38.9, 23.5 38.8, 23.4 38.9))"
polygon = ogr.CreateGeometryFromWkt(wkt)
```

我们接下来创建表示功能的 OGR Feature 对象，并按需设置几何形状和属性：

feature = ogr.Feature(layer.GetLayerDefn())
feature.SetGeometry(polygon)
feature.SetField("NAME", "My Polygon")

然后，我们可以将功能添加到图层中：
```
layer.CreateFeature(feature)
feature.Destroy()
```
注意到对 feature.Destroy() 的调用。这看起来可能有些奇怪，但这样做是释放由该功能使用的内存，并且巧合的是，它还会将功能的详细信息写入 shapefile。
最后，我们通过调用 Destroy() 方法来关闭目标文件：
```
dstFile.Destroy()
```
这将关闭目标文件并确保所有内容都已保存到磁盘。

与读取矢量格式数据一样，此代码不仅限于创建 shapefile。OGR 允许您以许多不同的格式创建地理空间数据，并且无论使用什么格式，您都使用相同的类和方法名称。

读取栅格数据

要读取栅格格式的地理空间数据，您使用 GDAL 库。让我们看看如何使用我们之前下载的 E 栅格数据来执行此操作。

确保将 e10g 和 e10g.hdr 文件都放在一个方便的目录中，然后在同一目录中创建一个 Python 脚本。我们将首先在这个脚本中输入以下内容：

from osgeo import gdal
dem_file = gdal.Open("e10g")

如您所见，我们使用 gdal.Open() 函数打开栅格数据文件。正如我们之前提到的，栅格数据源可以由多个数据带组成。要查看文件中有多少个数据带，您可以使用 RasterCount：

num_bands = dem_file.RasterCount

对于这个 DEM 文件，只有一个数据带；我们将使用 GetRasterBand() 方法获取对这个数据带的引用：

band = dem_file.GetRasterBand(1)

注意，数据带编号从 1 开始，而不是通常的 0。结果是 gdal.Band 对象。虽然您可以使用 Band 类的各种方法以原始字节序列的形式读取栅格数据带的内容，但从栅格数据带中提取数据的最简单方法是将其转换为 NumPy 数组：

data = band.ReadAsArray()

注意

如果您使用 GDAL 的 Mac OS X 安装程序，NumPy 将自动安装。在其他平台上，您可能需要自己安装它。NumPy 可以在 numpy.org 找到。

然后，您可以使用 NumPy 数组处理方法从该数组中提取数据。为了了解如何进行，让我们遍历数组并计算 DEM 数据中的高程值直方图：

num_rows,num_cols = data.shape

histogram = {} # Maps elevation to number of occurrences of that elevation.

for row in range(num_rows):
    for col in range(num_cols):
        elevation = int(data[row, col])
        try:
            histogram[elevation] += 1
        except KeyError:
            histogram[elevation] = 1

for elevation in sorted(histogram.keys()):
    print elevation, histogram[elevation]

如您所见，从 NumPy 数组中读取数据相当简单。

注意

我们程序中有一小部分可能让人困惑。请注意，我们使用：

elevation = int(data[row, col])

除了从 NumPy 数组中提取高程之外，我们还将它强制转换为整数。我们这样做是因为 data 是一个 NumPy 数组，它为数组中的每个条目返回一个 numpy.uint16 值。NumPy 会自动将其转换为整数，但这样做会减慢我们的程序速度。由于这些值已经是整数，我们直接将高程转换为常规整数。这提高了我们程序的速度大约一个数量级——当我们处理大量数据时这一点非常重要。

如果您运行此程序，您将看到一个唯一高程值列表以及该高程值在 DEM 文件中出现的频率：

% python readRaster.py
-500 53081919
-84 1
-83 8
-82 9
-81 17
...
5241 1
5295 1
5300 1
5443 1

注意负高程值。其中大部分是因为美国（例如死亡谷）的一些地区低于海平面。然而，有一个 -500 的高程值，这不是真实的高程值。这是我们之前提到的无数据值。

您可以通过在程序中添加以下突出显示的行来避免将其添加到直方图中：

...
histogram = {} # maps elevation to number of occurrences of that elevation.
no_data = int(band.GetNoDataValue())

for row in range(num_rows):
    for col in range(num_cols):
        elevation = int(data[row, col])
        if elevation == no_data: continue
        try:
            ...

以这种方式使用 NumPy，读取栅格格式数据源的内容相对简单。现在让我们看看写入栅格格式数据涉及哪些内容。

写入栅格数据

要写入栅格格式数据，我们需要生成一些样本数据，告诉 GDAL 如何将数据中的每个单元格地理配准到地球表面上的一个点，然后将数据保存到磁盘。让我们一步一步来完成这个任务。

我们首先创建栅格格式数据文件。我们将使用EHdr格式，它是 ESRI 头标签文件的简称——这是我们之前读取 DEM 数据时使用的相同文件格式。

注意

如同往常，GDAL 使得处理不同的数据格式变得容易；无论您选择什么格式，相同的代码都将工作。

下面是创建 EHdr 格式栅格数据文件的代码：
```
from osgeo import gdal
driver = gdal.GetDriverByName("EHdr")
dstFile = driver.Create("Example Raster", 180, 360, 1, gdal.GDT_Int16)
```
Create()方法的参数是文件名、单元格的宽度和高度、栅格波段数以及每个单元格中要存储的数据类型。
接下来，我们需要告诉 GDAL 为文件使用哪个空间参考系统。在这种情况下，我们将使用之前遇到的相同WGS84参考系统；如果您还记得，这意味着我们的坐标由纬度和经度值组成。以下是相关代码：
```
from osgeo import osr

spatialReference = osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("WGS84")

dstFile.SetProjection(spatialReference.ExportToWkt())
```
接下来，我们需要将栅格数据地理配准到地球表面。这是通过使用地理配准变换来完成的。在定义地理配准变换时，您可以使用许多选项，允许您执行诸如翻转栅格数据或旋转数据等复杂操作。然而，在这种情况下，我们只需要告诉 GDAL 左上角单元格应该放置的位置以及每个单元格的大小：
```
originX    = -180
originY    = 90
cellWidth  = 0.25
cellHeight = 0.25

geoTransform = [originX, cellWidth, 0, originY, 0, -cellHeight]
dstFile.SetGeoTransform(geoTransform)
```
在这个示例代码中，我们已将左上角单元格设置为latitude=90、longitude=-180，并定义每个单元格覆盖 1/4 度的纬度和经度。
我们现在准备好创建我们的栅格格式数据并将其保存到文件中。让我们生成一个 360 行和 180 列的数组，其中每个值是介于 1 和 100 之间的随机数：
```
import random

data = []
for row in range(360):
    row_data = []
    for col in range(180):
        row_data.append(random.randint(1, 100))
    data.append(row_data)
```
我们可以将这个数组转换为 NumPy 数组，其中数组中的每个条目都是一个 16 位有符号整数：
```
import numpy
array = numpy.array(data, dtype=numpy.int16)
```
然后，这些数据可以保存到文件中：
```
band.WriteArray(array)
```
最后，让我们定义一个无数据值，并关闭文件以将所有内容保存到磁盘：
```
band.SetNoDataValue(-500)
del dstFile
```

运行此程序将在磁盘上创建一个新的栅格格式文件，包括头部（.hdr）文件以及如何将我们的（随机）数据地理配准到地球表面的信息。除了添加空间参考系统和地理配准变换外，写入地理空间数据的过程几乎与读取一样简单。

注意

实际上，在从文件读取栅格数据时，您也可以使用空间参考系统和地理变换——我们只是跳过了这一步以保持简单。稍后，当我们想要将单元格精确地放置在地球表面的一个点上时，我们将在读取栅格格式数据时同时使用这两个概念。

处理空间参考系统

当你开始处理地理空间数据时，可能会相当令人困惑的一个概念是空间参考系统。想象一下，你正在执行一次搜救行动，并且被提供了一个飞机坠毁的坐标位置，例如：

(-114.93, 12.478)

这些数字代表什么？这些值是纬度和经度，还是可能是从给定参考点起数公里的距离？如果不了解这些坐标如何转换到地球表面上的一个点，你就无法知道该将救援人员送往何处。

注意

空间参考系统有时被称为坐标参考系统。不用担心：这两个术语指的是同一件事。

要理解空间参考系统的概念，你首先需要了解一些关于制图理论的知识。地图是尝试在二维笛卡尔平面上绘制地球的三维表面：

为了将地球表面转换成二维平面，你需要使用一种称为投影的数学过程。问题是，在数学上不可能有一个完美的投影：形状将会变形，面积将会被错误表示，或者点之间的距离将会不正确。

由于这种不完美，多年来已经开发了大量不同的地图投影。某些地图投影对于世界的某些地区来说非常准确，但在其他地方则不准确。其他地图投影在保持大陆形状的同时，错误地表示距离和面积，等等。

无论何时你处理地理空间数据，你都需要回答以下三个问题：

哪种数学模型被用来定义地球的形状？
坐标是否已经投影到地图上？
如果是这样，使用了哪种投影？

知道这些问题的答案将使你知道给定一组坐标所指的是确切位置。正如你可以想象的那样，知道这些问题的答案是任何地理空间分析成功的关键。

空间参考系统封装了这三个问题的答案。让我们看看几个常见的空间参考系统，看看它们是如何工作的。

WGS84

WGS84代表 1984 年世界大地测量系统，是一个全球标准，用于表示地球表面的点。它使用了一个精确的地球形状数学模型，以及定义纬度和经度的标准。综合来看，WGS84 空间参考系统提供了一个完整的系统来描述地球表面的点。

让我们更详细地看看 WGS84 是如何定义纬度和经度值的。给定地球表面上的一个点，纬度和经度是通过从地球中心到所需点的想象线来计算的：

然后，您可以测量纬度作为这条线与延伸到赤道的线之间的北南方向的角：

同样，经度可以通过这条线在东西方向上的角与延伸到零度的线（基于英国格林尼治的位置）来计算：

如您所见，经纬度值是基于地球表面上所需点的位置。WGS84 是一个典型的未投影坐标系统的例子。它是一种非常常见的地理空间数据格式，在许多情况下，你将只处理这种格式的数据。

普通横轴墨卡托

普通横轴墨卡托（UTM）是一个非常常用的标准，用于在平面直角坐标系上表示坐标。UTM 不是一个单一的地图投影，而是一系列称为区域的六十个不同的投影，每个区域覆盖地球表面的一小部分：

对于任何给定的 UTM 区域，坐标被测量为“北向”和“东向”值，这些值对应于给定参考点北或东方向上的米数。参考点被计算得使得北向和东向值始终为正。

由于 UTM 投影基于二维地图，这些是投影坐标系统的例子。

描述空间参考系统

无论何时您处理地理空间数据，您都需要知道您正在使用哪个空间参考系统。通常，在生成地图或读取和写入地理空间数据时，您需要构建一个 osr.SpatialReference 对象或其等效对象来描述您正在使用的空间参考系统。

描述空间参考系统的一种最简单的方法是使用名称。我们之前在创建空间参考对象时使用已知名称时看到了这一点，如下所示：

spatialReference = osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("WGS84")

描述空间参考系统的另一种常见方法是使用 EPSG 代码。EPSG 是一个维护所有已知空间参考系统数据库的标准机构，并为每个系统分配一个唯一的数字代码。您可以在 www.epsg-registry.org 找到 EPSG 网站。例如，WGS84 的 EPSG 代码是 4326，因此您也可以使用以下 Python 代码创建一个 WGS84 空间参考对象：

spatialReference = osr.SpatialReference()
spatialReference.ImportFromEPSG(4326)

最后，您可以使用 WKT 格式字符串来定义一个空间参考系统。GDAL/OGR 库使得使用 WKT 导入和导出空间参考系统变得容易。例如：

>>> spatialReference = osr.SpatialReference()
>>> spatialReference.ImportFromEPSG(4326)
>>> print spatialReference.ExportToWkt()
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

还有一个 ImportFromWkt() 方法，它允许你使用 WKT 定义字符串来定义一个空间参考对象。

坐标转换

除了知道你正在使用哪个空间参考系统之外，有时能够将空间数据从一个空间参考系统转换到另一个也很重要。例如，如果你想使用 Shapely 来计算两个多边形之间的交集，而这两个多边形使用不同的空间参考系统，你需要在交集工作之前将它们转换到同一个空间参考系统。

注意

记住，Shapely 是一个几何操作库。它不了解空间参考系统，所以你需要自己处理这个问题。

要将几何形状从一个空间参考系统转换到另一个，你可以使用osr.CoordinateTransformation类。让我们看看这是如何完成的：

src_spatialReference = osr.SpatialReference()
src_spatialReference.SetWellKnownGeogCS("WGS84")

dst_spatialReference = osr.SpatialReference()
dst_spatialReference.SetUTM(10)

transform = osr.CoordinateTransformation(src_spatialReference, dst_spatialReference)

geometry.Transform(transform)

你首先定义两个空间参考系统，然后创建坐标变换以将一个系统转换到另一个。然后你可以简单地调用Transform()方法将几何形状从源空间参考系统转换到目标空间参考系统。

计算长度和面积

现在我们已经了解了空间参考系统的基础知识以及如何将数据从一个空间参考系统转换到另一个，我们终于可以解决我们在第一章中遇到的问题，即地理空间分析和技术。如果你还记得，当我们查看使用 Shapely 库可以进行的计算时，我们发现我们无法准确计算使用经纬度值的空间数据的长度和面积。

让我们再次审视这个问题，以及我们如何可以使用坐标变换来解决这些问题。

让我们定义一个简单的多边形，它定义了纽约中央公园的大致轮廓：

import shapely.wkt
wkt = "POLYGON((-73.973057 40.764356, -73.981898 40.768094, -73.958209 40.800621, -73.949282 40.796853, -73.973057 40.764356))"
outline = shapely.wkt.loads(wkt)

我们在本章前面关于已知文本的部分使用了这个多边形作为例子。

如果我们要求 Shapely 计算这个几何形状的面积，它将计算这个多边形覆盖的数学面积：

>>> print outline.area
0.000377902804

不幸的是，得到的结果是“平方度”，这是一个没有意义的数字。这是因为 Shapely 不了解地图投影——它只是将每个坐标值视为一个数字。为了计算这个多边形在真实单位中的面积，我们必须将未投影的经纬度坐标转换为所谓的“等面积”地图投影，该投影以米为单位测量坐标。然后我们可以要求 Shapely 计算面积，结果将以平方米为单位。让我们看看我们如何使用 OGR 和 Shapely 的组合来完成这项工作：

首先，我们使用我们的轮廓的 WKT 定义创建一个 OGR 几何对象：
```
from osgeo import ogr
polygon = ogr.CreateGeometryFromWkt(wkt)
```

接下来，我们需要定义一个从 WGS84 到使用米的投影坐标系统的坐标变换。我们将使用世界莫勒维德投影（EPSG 代码 54009），这是一个在全球范围内相当准确的等面积投影：

from osgeo import osr

src_spatialReference = osr.SpatialReference()
src_spatialReference.ImportFromEPSG(4326)

dst_spatialReference = osr.SpatialReference()
dst_spatialReference.ImportFromEPSG(54009)

transform = osr.CoordinateTransformation(src_spatialReference, dst_spatialReference)

然后，我们可以将 OGR 几何体从 WGS84 转换为世界莫勒韦德投影，将其转换回 Shapely 几何体，并最终让 Shapely 计算多边形的面积：
```
polygon.Transform(transform)

outline = shapely.wkt.loads(polygon.ExportToWkt())
print outline.area
```

结果是中央公园面积的准确数值（与原始多边形轮廓允许的精度一样），以平方米为单位。然后你可以将这个面积转换为平方英里或其他任何你想要的单位。

在这个例子中，我们使用了等面积投影。为了准确计算长度，你必须使用一个等距地图投影来覆盖你感兴趣的地球区域。或者，你可以使用PyProj库来计算未投影坐标的距离；我们将在第五章分析地理空间数据中详细探讨PyProj。

地理空间数据错误及其修复方法

当你开始处理地理空间数据时，你很快会发现事情并不总是像你预期的那样运作。当尝试将几何体保存到形状文件中时，OGR 可能会崩溃，或者当计算两个多边形的交集时，Shapely 可能会引起系统错误。虽然这可能会令人沮丧，但一旦你了解了导致这些问题的原因，就有办法解决这些问题。

地理空间数据，以及像 GDAL/OGR 和 Shapely 这样的库，都是基于一个几何结构应该如何构建的数学模型。当你的地理空间数据不符合这个数学理想时，就会出现问题。让我们看看一个数学上正确的几何体是什么样的。

点

虽然坐标只是一对数字，但可接受值的范围是有限的。例如，想象以下使用 WGS84（即纬度和经度坐标）的点几何体：

POINT(-0.076 51.506)
POINT(2.295 48.858)
POINT(37.784 -122.402)

这些点本应代表伦敦塔、埃菲尔铁塔和旧金山莫斯 cone 中心的位置。然而，第三个点几何体被错误地定义了，通过交换纬度和经度。这个位置被设置为经度=37.784 和纬度=-122.402。但是纬度值只能在-90 到+90 的范围内，因此这个点几何体是无效的。

当然，所有几何体都是由坐标组成的，所以例如一个多边形可能只有一个超出范围的坐标，这可能导致你的程序崩溃。在构建或操作几何体时，你有时需要添加代码来检查坐标是否都有效，并在必要时调整几何体。

线段

线段几何体由一系列坐标组成，从一点到下一点绘制一条直线段：

然而，如果你尝试定义只有一个坐标的线段，或者两个坐标恰好相同的线段，那么你的线段几何体在数学上就是无效的，并可能导致你的程序崩溃。

线性环

线性环是起点和终点相同的 LineString：

线性环用于包围空间区域，是多边形几何体的构建块。为了使线性环有效，它必须至少有三个坐标，并且线段不能接触或交叉。

下图显示了两个数学上无效的线性环示例：

多边形

多边形几何体由一个或多个线性环组成：第一个线性环定义了多边形的轮廓，而额外的线性环定义了多边形内部的空间。例如：

然而，这种多边形的数学理想表示如果内部环重叠、相互接触或接触多边形的外部，就会崩溃。如果发生任何这些情况，您的多边形将变得无效，并且程序可能会崩溃。

多边形

如其名称所示，MultiPolygon 几何体是由两个或多个多边形组成的集合。如果其中两个多边形沿着边相接触，则 MultiPolygon 在数学上无效——在这种情况下，这两个多边形应该被合并成一个更大的多边形，因此 MultiPolygon 被认为是无效的。

修复无效几何体

即使您的几何数据一开始是有效的，在操作它们时几何体也可能变得无效。例如，如果您尝试将多边形分成两部分，或者合并 LineString 几何体，结果有时可能是无效的。

现在您已经了解了几何体可能无效的方式，让我们来看看一些修复它们的小技巧。首先，您可以通过检查 is_valid 属性来询问 Shapely 是否认为几何体是有效的：

if not geometry.is_valid:
    ....

同样，您可以使用 IsValid() 方法来检查 OGR 几何对象是否有效：

if not geometry.IsValid():
    ...

不幸的是，这两个有效性检查并不完美：有时您会发现一个几何体被标识为有效，尽管它实际上不是。当这种情况发生时，您必须在程序中添加 try...except 子句来捕获崩溃，然后尝试在再次尝试之前自己修复几何体。

当一个几何体无效时，您首先应该尝试的将是 buffer(0) 技术。buffer() 操作是一种将几何体扩展以包括原始几何体一定距离内的所有点的操作，例如：

通过调用 buffer(0)，您正在告诉 Shapely（或 OGR）构建一个新副本的几何体，该副本包括与几何体零距离内的所有点。这实际上是从头开始重建几何体，并且通常会将无效的几何体转换回有效的几何体。

很不幸，这并不总是有效的。有时buffer()函数无法在不崩溃的情况下重建复杂几何形状。在这种情况下，你可能需要将几何形状拆分成单个部分，然后逐一检查每个部分，看它是否是崩溃的原因。当你重新构建时，你可以排除表现不佳的部分。以下是一个示例 Python 代码，它尝试使用此技术修复无效的 Shapely 几何形状：

def fix_geometry(geometry):
    buffer_worked = True
    try:
        geometry = geometry.buffer(0)
    except:
        buffer_worked = False

    if buffer_worked:
        return geometry

    polygons = []
    if geometry.geom_type == "Polygon":
        polygons.append(geometry)
    elif geometry.geom_type == "MultiPolygon":
        polygons.extend(geometry.geoms)

    fixed_polygons = []
    for n,polygon in enumerate(polygons):
        if not linear_ring_is_valid(polygon.exterior):
            continue # Unable to fix.

        interiors = []
        for ring in polygon.interiors:
            if linear_ring_is_valid(ring):
                interiors.append(ring)

        fixed_polygon = shapely.geometry.Polygon(polygon.exterior,
                                                 interiors)

        try:
            fixed_polygon = fixed_polygon.buffer(0)
        except:
            continue

        if fixed_polygon.geom_type == "Polygon":
            fixed_polygons.append(fixed_polygon)
        elif fixed_polygon.geom_type == "MultiPolygon":
            fixed_polygons.extend(fixed_polygon.geoms)

    if len(fixed_polygons) > 0:
        return shapely.geometry.MultiPolygon(fixed_polygons)
    else:
        return None # Unable to fix.

def linear_ring_is_valid(ring):
    points = set() # Set of (x,y) tuples.

    for x,y in ring.coords:
        points.add((x,y))

    if len(points) < 3:
        return False
    else:
        return True

注意

记住，此代码与 Shapely 几何形状一起工作。如果你有一个 OGR 几何形状，你可以使用shapely.wkt.loads(ogrGeometry.ExportToWkt())将其转换为 Shapely 几何形状。

摘要

在本章中，我们更详细地研究了用于地理空间分析的数据。我们看到了为什么拥有高质量的地理空间数据很重要，你可能会遇到的地理空间数据的各种类型，以及提供免费高质量地理空间数据的主要网站。然后我们探讨了如何使用 GDAL 和 OGR 读取和写入矢量格式和栅格格式的地理空间数据，并了解了空间参考系统。最后，我们研究了地理空间数据可能变得无效的方式以及如何修复它。

在下一章中，我们将探讨空间数据库以及它们如何作为地理空间分析的有力工具。

第三章：空间数据库

在本章中，我们将探讨如何使用数据库存储、分析和操作地理空间数据。虽然空间数据库可能相当复杂，优化空间查询的过程可能具有挑战性，但它们可以以简单直接的方式使用，而无需太多麻烦，并且是地理空间分析师工具箱的重要组成部分。

在本章中，我们将：

学习在使用空间数据库之前需要了解的重要概念
在您的计算机上安装 PostgreSQL 关系数据库系统
安装 PostGIS 扩展到 PostgreSQL 以支持空间数据库
安装psycopg2数据库适配器，以便您可以从 Python 程序访问 Postgres
学习如何使用 PostGIS 创建空间数据库
发现如何使用 Python 将数据导入您的空间数据库
学习如何使用 Python 代码查询您的空间数据库
看看您如何从 Python 操作您的空间数据
学习如何从空间数据库导出数据

让我们从查看空间数据库的概念及其工作方式开始。

空间数据库概念

如前一章所述，空间数据库是可以存储和查询空间数据的数据库。在空间启用数据库表的每个记录中都有一个或多个几何字段，这些字段将记录定位在地球表面的某个位置。几何字段（们）的使用将取决于您在数据库表中存储的信息类型。例如：

表示送货车辆的记录可能包括一个表示车辆当前位置的点几何形状。
表示道路的记录可能包括一个表示道路形状的 LineString 几何形状。
表示森林火灾的记录可能包括一个表示火灾影响区域的 Polygon 几何形状。

注意

一些空间数据库允许您有多个几何字段，而其他数据库则限制每个记录只有一个。

单独来看，几何字段只是一个可以存储编码几何数据的数据库blob。数据通常以已知二进制（WKB）格式存储。这允许您从数据库中存储和检索几何数据。然而，单独来看，这并不非常有用——定义空间数据库的是使用存储的几何值构建空间索引的能力。

空间索引允许您根据记录在地球表面的位置在数据库中搜索记录。空间索引不会直接索引几何形状。相反，它为每个几何形状计算边界框，然后索引该边界框。以下插图显示了这是如何工作的：

空间索引的常见任务是指定包含给定点的几何形状（或几何形状）。例如，如果用户在地图上点击一个位置，你可能想知道用户点击的是哪个国家（如果有的话）。这可以通过以下空间数据库查询表示：

SELECT * FROM table WHERE ST_Contains(table.geometry, click_point);

注意

ST_Contains函数是一个数据库查询函数的例子。此函数由 PostGIS 空间数据库提供。不同的空间数据库使用不同的名称来命名它们的查询函数；本章中列出的所有查询函数都来自 PostGIS，因为我们将在本书中使用该数据库。

要执行此查询，数据库首先使用空间索引来识别那些包含所需点的边界框的记录。此过程如下所示：

十字准线代表所需点，矩形代表边界框。正如你所见，有两个边界框包含所需点。这些边界框对应于数据库中标记为法国和德国的记录。数据库使用这些信息将每个匹配的几何形状加载到内存中，并逐一检查它们是否包含所需点：

以这种方式，数据库能够确定点击点位于德国境内。

让我们回顾这个过程，因为它是一个非常重要的概念。数据库首先使用存储在空间索引中的边界框来识别可能匹配的记录，然后加载每个可能的几何形状到内存中进行检查。这个两步过程非常高效：通过使用空间索引中的边界框，它立即排除了绝大多数不可能匹配的记录。然后，它只对少数可能的匹配进行相对耗时的加载几何形状到内存的任务，然后逐一检查这些记录。

理解执行空间查询的两步过程非常重要，因为您必须做一些事情来确保它正常工作。特别是：

您必须确保您想要查询的几何形状包含在空间索引中。
您必须仔细措辞您的查询，以便数据库实际上可以使用您设置的索引。例如，如果数据库必须将您的几何形状从一种空间参考系统转换到另一种，或者在进行查询之前对数据进行某种空间操作，那么您的空间索引将被忽略，数据库将回退到对您所有数据进行顺序扫描。这可能非常慢，可能需要数小时甚至数天才能完成单个查询。
例如，如果你有一个具有大边界框的极其复杂的几何形状，比如美国的一个详细轮廓，你可能会发现你的查询仍然需要很长时间才能完成。这是因为边界框覆盖了地球表面如此大的区域，它被包含在许多查询中，而轮廓的复杂性意味着查询仍然需要很长时间来处理。解决这个问题的方法之一是将一个大而复杂的几何形状分成更小的部分，这样数据库只需要处理一小部分而不是整个。

尽管存在这些潜在问题，空间数据库是存储和分析地理空间数据的绝佳工具。当然，空间数据库不仅限于使用 ST_Contains() 搜索记录。它们可以用于各种空间查询，如下表所示：

空间查询函数	描述
`ST_Within`	这匹配完全被给定多边形包围的几何形状的记录。
`ST_Intersects`	这匹配记录的几何形状与给定几何形状相交的记录。
`ST_Crosses`	这匹配记录的几何形状与给定的线或多边形交叉的记录。
`ST_DWithin`	这匹配位于给定位置或几何形状给定距离内的记录。

这些空间查询函数中有些细微之处你需要熟悉——这些在 PostGIS 文档中有详细描述。然而，这个表格应该能让你了解空间数据库的力量，并告诉你如何使用适当的空间索引的空间数据库可以是一个处理地理空间数据的强大工具，尤其是在你需要处理大量记录时。

现在你已经对空间数据库的工作原理有了些了解，让我们在你的电脑上安装一个，然后看看我们如何从你的 Python 程序中访问它。

安装空间数据库

在这本书中，我们将使用最受欢迎和最强大的地理空间数据库之一：PostGIS。PostGIS 是一个免费提供的 PostgreSQL 关系型数据库的扩展。为了在我们的 Python 程序中使用它，我们需要安装三个独立的软件组件：

The PostgreSQL database server itself
PostgreSQL 的 PostGIS 扩展
The psycopg2 database adapter for Python

注意

PostgreSQL 通常简称为 Postgres。在这本书中，我们将经常使用这个更口语化的名称。

让我们逐一了解安装这些软件组件的过程。

安装 PostgreSQL

PostgreSQL (postgresql.org) 是最强大的开源关系型数据库之一。虽然它以设置和使用困难而闻名，但实际上并不太复杂，并且由于每个主要操作系统都有预构建的安装程序，设置过程现在相当直接。

让我们继续在您的计算机上安装 PostgreSQL。您如何操作取决于您正在运行的操作系统：

如果您的计算机运行的是 Microsoft Windows，您可以从 www.enterprisedb.com/products-services-training/pgdownload 下载 PostgreSQL 的安装程序。选择适合您 Windows 版本（32 位或 64 位）的安装程序，并下载安装程序文件。然后只需双击下载的安装程序，并按照说明操作。
如果您正在运行 Mac OS X，您可以从 KyngChaos 网站下载 PostgreSQL 的有效版本，www.kyngchaos.com/software/postgres。只需下载磁盘映像，打开它，然后双击 PostgreSQL.pkg 包文件来在您的计算机上安装 PostgreSQL。
如果您正在使用 Linux 机器，您可以遵循 PostgreSQL 下载页面上的说明，www.postgresql.org/download。选择适合您使用的 Linux 发行版的链接，您将看到相应的安装说明。

安装 PostgreSQL 后，您可以通过在终端或命令行窗口中输入 psql 命令来检查它是否正在运行。如果一切顺利，您应该会看到 Postgres 命令行：

psql (9.3.4)
Type "help" for help.

postgres=#

小贴士

如果 psql 命令在用户身份验证方面出现问题，您可能需要指定在连接到 Postgres 时使用的用户账户。例如：

% psql –U postgres

许多 PostgreSQL 安装都有一个名为 postgres 的用户，当访问数据库时，您需要选择该用户（使用 –U 命令行选项）。或者，如果您正在运行 Microsoft Windows，您可能需要使用 sudo 切换到 root 用户，或者以管理员身份打开命令提示符。

安装 PostGIS

现在我们已经安装了 Postgres 本身，接下来我们需要安装 PostGIS 空间数据库扩展。PostGIS 的主要网站可以在 postgis.net 找到。您应该访问此网站，点击文档选项卡，并下载最新版本 PostGIS 的用户手册。您会发现这本手册非常有帮助，因为它详细解释了 PostGIS，包括您可以执行的各类查询。

您如何安装 PostGIS 取决于您正在运行的操作系统：

如果您的计算机运行的是 MS Windows，您可以从 download.osgeo.org/postgis/windows 下载 PostGIS 的安装程序。
对于 Mac OS X，从 kyngchaos.com/software/postgres 下载并运行 PostGIS 安装程序。

注意

注意，您还需要安装 GDAL 完整包，您应该在处理上一章时已经完成了这一步骤。
如果您使用的是基于 Linux 的操作系统，请遵循 PostGIS 安装页面上的说明：postgis.net/install。

要检查 PostGIS 是否已成功安装，请尝试在终端窗口中输入以下命令序列：

% createdb test_database
% psql -d test_database -c "CREATE EXTENSION postgis;"
% dropdb test_database

小贴士

如果您需要在不同的用户账户下运行 PostgreSQL，则需要为这些命令添加 –U postgres 选项或使用 sudo。

如您可能猜到的，createdb 命令创建一个新的数据库。然后我们使用 psql 命令使用 PostGIS 扩展初始化该数据库，最后使用 dropdb 命令再次删除数据库。如果这一系列命令运行无误，则表示您的 PostGIS 安装（以及 Postgres 本身）已正确设置并运行。

安装 psycopg2

现在我们已经拥有了一个空间数据库，让我们安装 psycopg2 库，以便我们可以使用 Python 访问它。

注意事项

psycopg2 是一个标准的 Python 数据库适配器——也就是说，它是一个符合 PEP 249 中指定的 Python 数据库 API 的库（www.python.org/dev/peps/pep-0249）。我们将探讨如何使用 psycopg2 存储和查询空间数据，但如果您之前没有使用过 Python 数据库适配器，您可能想查看有关此主题的可用教程。有关此主题的一个很好的教程可以在 halfcooked.com/presentations/osdc2006/python_databases.html 找到。

psqcopg2 的网站可以在 initd.org/psycopg 找到。通常，您安装此库的方式取决于您使用的操作系统：

对于 MS Windows 系统，您可以从 www.stickpeople.com/projects/python/win-psycopg 下载可双击安装的程序包。
如果您的计算机运行的是 Mac OS X 系统，可以在 www.kyngchaos.com/software/python 找到可双击安装的程序包。
对于 Linux 机器，您需要从源代码安装 psycopg2。有关如何操作的说明，请参阅 initd.org/psycopg/docs/install.html。

要检查是否安装成功，请启动您的 Python 解释器并输入以下命令：

>>> import psycopg2
>>>

如果 psycopg2 安装正确，您应该会看到 Python 解释器提示符重新出现，且没有错误信息，如本示例所示。如果出现错误信息，您可能需要遵循 psycopg2 网站上的故障排除说明。

从 Python 访问 PostGIS

到目前为止，我们已经将一些工具和库安装到了您的计算机上。现在，是时候使用这些工具和库做一些有趣的事情了。在本章的剩余部分，我们将把世界边界数据集导入到名为 world_borders 的 PostGIS 数据库中，然后使用 Python 对该数据进行各种查询。我们还将了解如何使用 PostGIS 和 Python 操作该数据集。

首先，创建一个名为 world_borders 的新目录并将其放置在方便的位置。您将使用此目录来存储您创建的各种文件。

设置空间数据库

当使用 psycopg2 访问数据库时，我们首先必须指定我们要使用哪个数据库。这意味着在您的 Python 代码可以使用它之前，数据库必须存在。为了设置一切，我们将使用 Postgres 命令行工具。在终端或命令行窗口中输入以下内容：

% createdb world_borders

提示

不要忘记包含 -U postgres 选项，或者如果需要以不同的用户账户访问 Postgres，请使用 sudo。

这创建了数据库本身。接下来，我们想要为我们的数据库启用 PostGIS 空间扩展。为此，请输入以下命令：

% psql -d world_borders -c "CREATE EXTENSION postgis;"

现在我们已经设置了数据库本身，让我们在数据库中创建一个表来存储我们的空间数据。为此，我们将创建一个名为 create_table.py 的 Python 程序。请创建一个名为 world_borders 的目录，并在其中创建此文件，然后输入以下内容到文件中：

import psycopg2

我们现在想要打开到数据库的连接。为此，我们必须告诉 psycopg2 我们要使用哪个数据库以及哪个用户账户（以及可能，哪个密码）来访问该数据库。这是通过向 psycopg2.connect() 函数提供关键字参数来完成的，如下所示：

connection = psycopg2.connect(database="world_borders", user="...", password="...")

提示

您只需要 user 参数，如果您在运行 Postgres 命令行工具时需要提供 -U 命令行参数。您也只需要 password，如果该用户账户受密码保护。

一旦我们建立了数据库连接，我们接下来设置一个 cursor 对象，我们将使用它向数据库发出命令：

cursor = connection.cursor()

下一步可能有点反直觉：我们不会创建数据库表，而是会删除它如果它已经存在。这样做可以让我们多次运行 create_table.py 脚本而不会产生任何错误。以下是相关代码：

cursor.execute("DROP TABLE IF EXISTS borders")

execute() 语句告诉游标运行给定的 SQL 命令。在这种情况下，命令是 DROP TABLE IF EXISTS，它告诉数据库如果表已存在，则删除（删除）该表。

我们现在可以使用以下命令创建我们的数据库表：

cursor.execute("CREATE TABLE borders (" +
                   "id SERIAL PRIMARY KEY," +
                   "name VARCHAR NOT NULL," +
                   "iso_code VARCHAR NOT NULL," +
                   "outline GEOGRAPHY)")

注意，我们将此命令拆分到多行，以便更容易阅读。除了最后一行之外，这是一个标准的 SQL 数据库表定义：我们正在创建一个表，其中每个记录都有一个由数据库自动分配的唯一id值，一个name值和一个iso_code值。在最后一行中，我们创建outline字段，并给它一个GEOGRAPHY类型。地理字段是 PostGIS 特有的；它们是GEOMETRY字段类型的一种变体，旨在与使用未投影纬度和经度坐标的空间数据一起工作。

既然我们已经创建了我们的数据库表，让我们为这些数据设置一个空间索引。正如我们所见，空间索引将大大加快对数据库的查询速度。让我们为我们的outline字段创建一个空间索引：

cursor.execute("CREATE INDEX border_index ON borders USING GIST(outline)")

最后，因为 Postgres 是一个事务型数据库，我们需要提交我们所做的更改，以使它们永久。以下是完成此操作的必要代码：

connection.commit()

这完成了我们的create_table.py程序，它应该看起来像下面这样：

import psycopg2

connection = psycopg2.connect(database="world_borders", user="...", password="...")
cursor = connection.cursor()

cursor.execute("DROP TABLE IF EXISTS borders")

cursor.execute("CREATE TABLE borders (" +
                   "id SERIAL PRIMARY KEY," +
                   "name VARCHAR NOT NULL," +
                   "iso_code VARCHAR NOT NULL," +
                   "outline GEOGRAPHY)")

cursor.execute("CREATE INDEX border_index ON borders USING GIST(outline)")
connection.commit()

如果你运行这个程序，你的数据库表和相关的空间索引将被创建。现在，让我们将世界边界数据集的内容导入我们新创建的表中。

导入空间数据

复制你之前下载的TM_WORLD_BORDERS-0.3目录，并将其放置在你的world_borders目录内。然后创建另一个名为import_data.py的 Python 脚本。这就是你将放置将数据导入数据库的代码的地方。

我们将使用 OGR 库从 shapefile 导入数据，并使用psycopg2将其插入数据库。因此，我们程序的前两行应该看起来像下面这样：

import osgeo.ogr
import psycopg2

接下来，我们需要打开到数据库的连接。执行此操作的代码与我们在create_table.py脚本中使用的代码相同：

connection = psycopg2.connect(database="world_borders", user="...", password="...")
cursor = connection.cursor()

小贴士

不要忘记调整psycopg2.connect()的关键字参数，以匹配你需要连接到 PostgreSQL 的用户账户。

我们现在准备开始从 shapefile 导入数据。不过，首先我们将删除我们数据库表中的现有内容；这将使我们能够多次运行我们的import_data.py程序，在添加新记录之前清除现有记录，这样我们每次都是从一张白纸开始：

cursor.execute("DELETE FROM borders")

我们现在准备从 shapefile 导入数据到数据库中。让我们先打开 shapefile，并逐个提取我们想要的信息：

shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature  = layer.GetFeature(i)
    name     = feature.GetField("NAME")
    iso_code = feature.GetField("ISO3")
    geometry = feature.GetGeometryRef()

这应该对你来说很熟悉，因为我们已经在上一章中使用了 OGR 来读取 shapefile 的内容。现在我们有了几何形状，我们可以将其转换为 WKT 格式，如下所示：

    wkt = geometry.ExportToWkt()

我们现在有了将要素插入数据库所需的所有信息。以下是执行实际插入的代码：

    cursor.execute("INSERT INTO borders (name, iso_code, outline) VALUES (%s, %s, ST_GeogFromText(%s))", (name, iso_code, wkt))

这里有很多事情在进行中，所以让我们更仔细地看看这个命令。我们在这里使用 INSERT，这是一个标准的 SQL 命令。INSERT 命令有以下基本结构：

INSERT INTO table (field, field, ...) VALUES (value, value, ...);

如您所见，我们指定了数据库表名、字段列表以及要存储到这些字段中的值。

作为标准的 Python 数据库适配器，psycopg2 将自动将 Python 值（如整数、浮点数、字符串、datetime 对象等）转换为它们的 SQL 等价物。这就是那些 %s 占位符的作用所在——我们在 SQL 命令字符串的每个想要提供值的地方使用 %s，然后作为 cursor.execute() 命令的第二个参数提供实际的值。例如，考虑以下 Postgres 命令：

cursor.execute("INSERT INTO users (name, age) VALUES (%s, %s)", (user_name, user_age))

这个命令将记录插入到 users 表中，将 name 字段设置为 user_name 变量的值，并将 age 字段设置为 user_age 变量的值。这种将 Python 值转换为 SQL 字面量的转换非常强大，并且是使用数据库适配器的主要好处之一。

在我们将 shapefile 的内容导入到 borders 表的 INSERT 语句中，有一个最后的复杂性：我们使用 ST_GeogFromText() 函数在插入到 outline 字段之前将我们的 WKT 格式字符串转换为地理值。我们必须这样做，因为 OGR 和 Postgres 使用不同的内部表示形式来表示几何数据。WKT 格式字符串是这些两种内部表示之间的 通用语言。

在我们从 shapefile 导入各种特征之后，我们必须提交我们的更改，以便它们被写入数据库：

connection.commit()

将所有这些放在一起，这就是我们的 import_data.py 程序的样子：

import osgeo.ogr
import psycopg2

connection = psycopg2.connect(database="world_borders", user="...", password="...")
cursor = connection.cursor()

cursor.execute("DELETE FROM borders")

shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    iso_code = feature.GetField("ISO3")
    geometry = feature.GetGeometryRef()
    wkt = geometry.ExportToWkt()

    cursor.execute("INSERT INTO borders (name, iso_code, outline) VALUES (%s, %s, ST_GeogFromText(%s))", (name, iso_code, wkt))

connection.commit()

当我们运行这个程序时，World Borders Dataset shapefile 中的所有记录都应导入到数据库中。请注意，这只需要几秒钟就能完成——尽管我们必须将轮廓从 OGR 几何形状转换为 WKT，然后再从 WKT 转换为 PostGIS 地理对象，但这并不需要很长时间。

如果您愿意，可以运行 psql 命令行客户端，并输入如 SELECT id,name,iso_code FROM borders 这样的命令来查看您已导入的数据。但当然，我们真正想要的是使用 Python 来查询我们的空间数据库。现在让我们这样做。

查询空间数据

让我们再写一个 Python 程序来对数据库内容执行各种查询。首先创建另一个名为 query_data.py 的 Python 文件，并将其放置在 world_borders 目录中。我们首先导入 psycopg2 库，打开到我们数据库的连接，并创建一个数据库游标：

import psycopg2
connection = psycopg2.connect(database="world_borders", user="...", password="...")
cursor = connection.cursor()

这一切都应该来自我们之前创建的 create_table.py 程序。

现在让我们执行一个简单的（非空间）数据库查询，只是为了看看它是如何工作的。将以下内容添加到你的程序末尾：

cursor.execute("SELECT id,name FROM borders ORDER BY name")
for row in cursor:
    print row

当你运行你的 query_data.py 程序时，你应该看到来自你的 borders 表的记录 ID 和相关名称的列表：

(1264, 'Afghanistan')
(1237, 'Albania')
(1235, 'Algeria')
...

注意，你使用 cursor.execute() 来执行你的查询，然后遍历游标以获取匹配的行。每行的值是一个包含你请求的字段的元组。

注意

当然，你也可以使用 %s 在你的查询中包含 Python 值，例如：

cursor.execute("SELECT id FROM borders WHERE name=%s", (country_name,))

到目前为止，我们一直在使用 PostgreSQL 的非空间方面。现在让我们针对这些数据执行一个空间查询。我们将要求数据库识别所有距离巴黎 1,000 公里的国家。使用 PostGIS 中的 GEOGRAPHY 数据类型，这很容易做到：

lat    = 48.8567
long   = 2.3508
radius = 1000000

cursor.execute("SELECT name FROM borders WHERE ST_DWITHIN(" +
                   "ST_MakePoint(%s, %s), outline, %s)", (long, lat, radius))
for row in cursor:
    print row[0]

ST_DWithin 命令识别距离指定点 radius 米内的国家；运行程序应该返回一个列表，列出距离巴黎 1,000 公里的国家：

San Marino
Denmark
Ireland
Austria
...
Switzerland
United Kingdom

这让你了解 PostGIS 的强大功能以及你可以使用 psycopg2 数据库适配器执行的查询类型。确保你研究 PostGIS 手册中的 PostGIS 参考部分，以了解你可以执行的各种类型的空间查询。

操作空间数据

在你的空间分析中，你不仅限于使用静态数据。你还可以在 PostGIS 数据库中直接创建新几何形状并操作现有几何形状。虽然使用我们之前使用的 ST_GeogFromText() 函数等函数创建全新的几何形状很容易，但你也可以使用 PostGIS 几何编辑和几何处理函数从旧几何形状派生新的地理值。

注意

当你使用 PostGIS 的 GEOGRAPHY 类型时，有一些函数可用性的限制。PostGIS 最初只支持 GEOMETRY 数据类型，该数据类型设计为仅与投影到平坦笛卡尔平面的空间数据一起工作。当使用 GEOGRAPHY 字段时，请查看 PostGIS 手册以了解哪些函数受支持。

为了了解我们如何根据现有数据计算新的空间值，让我们编写一个程序来缓冲我们的轮廓，并将它们存储到我们数据库表中的新 GEOGRAPHY 列中。

我们在上一章中看到了 buffer() 操作，我们看到了它通常可以用来修复无效的几何形状。如果你还记得，buffer() 操作构建了一个新的几何形状，它包括现有几何形状一定距离内的所有点。以下图像显示了英国的轮廓，以及缓冲后的相同轮廓：

让我们编写一个程序来计算这些缓冲轮廓。在你的 world_borders 目录中创建一个新的 Python 脚本，并将其命名为 buffer.py。将以下内容输入到这个文件中：

import psycopg2
connection = psycopg2.connect(database="world_borders", user="...", password="...")
cursor = connection.cursor()

我们现在想要创建一个新的字段来存储缓冲的大纲。为此，请将以下内容添加到文件末尾：

try:
    cursor.execute("ALTER TABLE borders ADD COLUMN buffered_outline GEOGRAPHY")
except psycopg2.ProgrammingError:
    connection.rollback()

ALTER TABLE命令是一个标准的 Postgres 命令，用于更改数据库表的结构；在这种情况下，我们添加一个名为buffered_outline的新GEOGRAPHY列。

注意

注意，我们将我们的ALTER TABLE命令包裹在try...except语句中。这是因为如果列已经存在，psycopg2将引发ProgrammingError。通过捕获这个错误，我们可以多次运行我们的buffer.py程序而不会失败，因为buffered_outline字段已经添加到表中。

由于psycopg2中异常的事务问题，当发生异常时，我们必须调用connection.rollback()。这允许程序即使在抛出异常的情况下也能继续运行。

我们下一个任务是计算缓冲的大纲。使用 PostGIS，这非常简单：

cursor.execute("UPDATE borders SET buffered_outline=ST_Buffer(outline, 1000)")

在这个 SQL 语句中，我们将buffered_outline字段的值设置为ST_Buffer()命令的结果。ST_Buffer()命令接受一个地理值和一个以米为单位的距离；它返回一个包含所有在给定距离内的现有地理点的新的地理值。

我们最后的任务是提交我们对数据库所做的更改：

connection.commit()

这实际上完成了我们的buffer.py程序，如果我们运行它，我们将得到存储在buffered_outline字段中的每个大纲的缓冲版本。然而，因为这个程序没有显示任何内容，所以无法知道它是否真的工作。为了解决这个问题，让我们添加一个最终的空问查询来计算并显示每个大纲的面积。

我们查询的基本结构将如下所示：

cursor.execute("SELECT name, ST_Area(outline), ST_Area(buffered_outline) FROM borders ORDER BY name")
for name, area1, area2 in cursor:
    ...

ST_Area()函数的结果是以平方米为单位的地理面积。因为这些数字可能非常大，我们希望将它们转换为平方千米进行显示。然而，有一个小问题：当我们缓冲几何形状时，它有时可能变得无效，因为缓冲的几何形状超出了有效的纬度和经度值范围。尽管我们只缓冲了千米，但任何接近北极或南极，或接近经度-180 或+180 度极限的地理形状都将有一个无效的缓冲轮廓。当我们尝试计算这些无效轮廓的面积时，结果将是一个NaN（不是一个数字）值。

让我们添加一些代码来检查无效的面积并优雅地处理它们；将之前代码列表中的...行替换为以下内容：

    if not math.isnan(area1):
        area1 = int(area1/1000000)
    else:
        area1 = "n/a"
    if not math.isnan(area2):
        area2 = int(area2/1000000)
    else:
        area2 = "n/a"
    print name, area1, area2

您还需要在程序顶部添加一个import math语句。

运行这个程序需要大约一分钟的时间来计算所有缓冲区，之后将显示计算出的面积：

Afghanistan 641915 646985
Albania 28676 29647
Algeria 2317478 2324740
American Samoa 229 363
...
Zimbabwe 389856 392705
Åland Islands 817 1144

如您所见，缓冲区域比原始区域略大，这是您所期望的。

导出空间数据

我们对空间数据库的介绍几乎完成；唯一剩下要检查的是如何再次从 PostGIS 中获取空间数据，例如将其保存回 shapefile。要从 GEOGRAPHY 字段中提取空间值，请使用 ST_AsText() 函数。例如：

cursor.execute("SELECT name,ST_AsText(outline) FROM borders")
for name,wkt in cursor:
    geometry = osgeo.ogr.CreateGeometryFromWkt(wkt)
    ...

您可以使用 OGR 几何对象将空间数据写入 shapefile，或者对它进行任何您希望执行的操作。

摘要

在本章中，我们探讨了空间数据库如何成为地理空间数据分析的强大工具。我们介绍了空间数据库背后的重要概念，并在您的计算机上安装了 PostgreSQL、PostGIS 和 psycopg2。然后，我们通过创建空间数据库、将数据导入该数据库、执行空间查询、使用 PostGIS 操作空间数据以及使用 Python 代码从空间数据库导出数据，亲自动手实践。

在下一章中，我们将探讨如何使用 Mapnik 库根据我们的地理空间数据制作出外观精美的地图。

第四章：创建地图

在本章中，我们将探讨 Python 程序如何使用 Mapnik 库创建外观精美的地图。你将在计算机上安装 Mapnik，学习 Mapnik 库的基础知识，并了解如何使用它生成简单的地图。然后，我们将探索 Mapnik 的更多高级特性，并了解它如何用于产生各种复杂的视觉效果。最后，我们将创建一个有用的程序，该程序将任何形状文件的内容以地图的形式显示出来。

介绍 Mapnik

如果不能可视化地理空间数据，很难理解其含义。通常，通过绘制地图来使空间数据可见——实际上，地图只不过是由空间数据创建的图像。Mapnik (mapnik.org) 是一个强大的工具，可以将原始地理空间数据转换为地图图像。

Mapnik 本身是用 C++ 编写的，但它附带了一些绑定，允许你从 Python 访问它。使用 Python 代码，你可以定义组成地图的各个图层，指定包含要显示的数据的数据源，然后设置样式，以控制如何绘制各种特征。

当你刚开始使用 Mapnik 时，它可能会让你感到有些难以接近，所以让我们直接动手，立即开始。首先，让我们在你的计算机上安装 Mapnik，并使用它生成一个简单的地图，然后再深入探讨如何使用 Mapnik 库构建和样式化地图。

安装 Mapnik

要安装 Mapnik，请访问 Mapnik 主网站上的下载页面 (mapnik.org/pages/downloads.html)，并选择适合你操作系统的安装程序。对于 Mac OS X 和 MS Windows，都有预构建的包可用。对于 Linux 计算机，你需要从源代码编译程序，或者使用包管理器下载、编译和安装 Mapnik 及其各种依赖项；有关如何操作的完整说明可在 Mapnik 下载页面上找到。

Mapnik 的初体验

我们将通过编写一个简单的程序来开始对 Mapnik 的探索，该程序使用我们之前下载的世界边界数据集生成地图。将 TM_WORLD_BORDERS-0.3 目录的内容复制到方便的位置，然后在同一目录中创建一个新的 Python 程序。将你的新程序命名为 mapnik_example.py。此程序将根据世界边界数据集形状文件的内容生成一个 PNG 格式的图像文件。

将以下内容输入到你的 mapnik_example.py 文件中：

import mapnik

map = mapnik.Map(1200, 600)
map.background = mapnik.Color("#e0e0ff")

layer = mapnik.Layer("countries")
layer.datasource = mapnik.Shapefile(file="TM_WORLD_BORDERS-0.3.shp")
layer.styles.append("country_style")
map.layers.append(layer)

fill_symbol = mapnik.PolygonSymbolizer(mapnik.Color("#60a060"))
line_symbol = mapnik.LineSymbolizer(mapnik.Color("black"), 0.5)

rule = mapnik.Rule()
rule.symbols.append(fill_symbol)
rule.symbols.append(line_symbol)

style = mapnik.Style()
style.rules.append(rule)

map.append_style("country_style", style)

map.zoom_all()
mapnik.render_to_file(map, "map.png", "png")

当你运行此程序时，应在同一目录中创建一个名为 map.png 的新文件。打开此文件将显示生成的地图：

现在我们已经看到了我们的示例程序是如何工作的，让我们仔细看看它，并依次检查每个部分。让我们从程序的开始部分开始：

import mapnik

map = mapnik.Map(1200, 600)
map.background = mapnik.Color("#e0e0ff")

在这里，我们简单地导入 Mapnik 库，然后创建并初始化一个新的地图对象。地图图像宽度为 1,200 像素，高度为 600 像素，地图背景为浅蓝色，由十六进制颜色值 #e0e0ff 定义。

一个地图由一个或多个 地图层 组成。在我们的程序中，我们只有一个地图层，我们将其设置为访问 TM_WORLD_BORDERS shapefile：

layer = mapnik.Layer("countries")
layer.datasource = mapnik.Shapefile(file="TM_WORLD_BORDERS-0.3.shp")
layer.styles.append("country_style")
map.layers.append(layer)

在这个层定义中，有几个需要注意的地方：

每个地图层都有一个名称，该名称在地图中唯一标识了层；在我们的程序中，我们称我们的地图层为 countries。
每个层都有一个 数据源，它告诉 Mapnik 数据应该从哪里获取。在这种情况下，我们使用 mapnik.Shapefile 类从 shapefile 加载数据，尽管可以使用许多不同类型的数据源。例如，您可以直接从空间数据库加载数据，甚至可以使用 Python 数据源以编程方式创建和显示要素。
layer.styles.append("country_style") 这一行告诉 Mapnik 应该使用哪个风格来绘制层的数据。Mapnik 风格通过名称引用，并且可以与每个层关联任意数量的风格。

注意

Mapnik 层还可以与一个空间参考系统相关联。如果您没有指定空间参考系统，Mapnik 将假设数据位于标准的 EPSG 4326 空间参考系统中。

接下来，我们想要定义 country_style 风格，该风格将绘制我们的地图层的内容。一个风格由任意数量的规则组成，其中每个规则都有一个可选的 过滤器，用于标识数据源中哪些要素应该使用此规则绘制，以及一个用于将匹配的要素绘制到地图上的 符号化器 列表。

我们首先创建两个符号化器：一个用于用淡绿色填充每个多边形的内部，另一个用于用细黑线绘制每个多边形的轮廓：

fill_symbol = mapnik.PolygonSymbolizer(mapnik.Color("#60a060"))
line_symbol = mapnik.LineSymbolizer(mapnik.Color("black"), 0.5)

对于填充符号，我们再次使用十六进制颜色代码来定义用于绘制多边形内部的颜色，而对于线符号，我们使用一个命名颜色。请注意，0.5 值定义了绘制每个多边形轮廓的宽度，单位为像素。

现在我们已经有了我们的两个符号化器，接下来我们定义一个规则，使用它们来绘制 shapefile 的内容：

rule = mapnik.Rule()
rule.symbols.append(fill_symbol)
rule.symbols.append(line_symbol)

注意，这个规则没有过滤器；过滤器的缺失告诉 Mapnik 应该使用这两个符号化器绘制层数据源中的每个要素。

为了完成定义我们的 country_style 风格，我们初始化 Style 对象本身，将我们唯一的规则添加到风格中，然后将风格对象添加到我们的映射中：

style = mapnik.Style()
style.rules.append(rule)

map.append_style("country_style", style)

注意，我们在将风格添加到地图对象时给出了一个名称；因为这个名称用于标识地图层使用的风格，所以当我们添加风格到地图以及引用地图层中的风格时，使用完全相同的名称是很重要的。

我们最后的任务是告诉地图要显示世界的哪个区域，以及如何将地图的可见部分渲染成图像文件：

map.zoom_all()
mapnik.render_to_file(map, "map.png", "png")

在我们的示例中，我们放大以显示地图层中的所有数据，并将结果保存为名为map.png的 PNG 格式图像文件。

这就完成了我们使用 Mapnik 生成地图图像的示例 Python 程序。你可以使用 Mapnik 做很多更复杂的事情，但这将给你一个了解它是如何工作以及你可以用它做什么的思路。

构建地图

现在我们已经看到了 Mapnik 的一个示例，让我们更仔细地看看 Mapnik 库背后的某些理念。例如，考虑以下美国西海岸的地图：

这张地图实际上由四个不同的地图层组成：

正如你所见，地图层是逐层叠加绘制的，如图表右侧的箭头所示。为了达到正确的视觉效果，层需要以相反的顺序添加，这样添加的每一层都会出现在地图中已有的层之前。也就是说，基础层应该首先添加，然后是城市区域层，依此类推。层添加到地图中的顺序非常重要；如果你顺序错误，某些层可能会被遮挡。

地图样式化

正如我们之前看到的，地图的样式是通过创建mapnik.Style对象并将其添加到地图中定义的，每个对象都有一个独特的名称：

style = mapnik.Style()
# ...setup style
map.append_style("style_name", style)

我们通过将样式名称添加到层的styles列表中，告诉每个地图层我们希望该层使用哪些样式：

layer.styles.append("style_name")

你可能会认为直接将样式定义添加到地图层会更简单，但通过名称引用样式的这个过程是有意为之的：它将显示的内容与显示的方式分离开来。这种方法允许你在多个地图层中使用相同的样式集，或者只需通过交换样式名称就可以完全改变地图的外观。

注意

定义你的地图样式有另一种方式。而不是创建自己的mapnik.Style对象并逐个将其添加到地图中，你可以使用 XML 格式的样式表一次性定义所有样式。虽然这非常强大，但 XML 样式表很难阅读，而且非常不符合 Python 风格。出于这些原因，我们在这本书中不会使用 XML 样式表。

在我们的示例程序中，我们创建了一个单一的mapnik.Style对象，它只包含一条规则。这条规则有两个符号化器与之相关联，告诉 Mapnik 如何绘制每个国家多边形的内部和外部。然而，规则可以更加复杂，包括在规则被使用之前必须满足的各种条件。例如，考虑以下 Python 代码片段：

symbol1 = ...
symbol2 = ...
symbol3 = ...

rule1 = mapnik.Rule()
rule1.filter = mapnik.Filter("[POPULATION] < 500000")
rule1.symbols.append(symbol1)

rule2 = mapnik.Rule()
rule.filter = mapnik.Filter("[POPULATION] < 1000000")
rule2.symbols.append(symbol2)

rule3 = mapnik.Rule()
rule3.set_else(True)
rule3.symbols.append(symbol3)

style.rules.append(rule1)
style.rules.append(rule2)
style.rules.append(rule3)

由于样式的规则是依次评估的，因此当特征的POPULATION属性值小于 500,000 时，此样式将使用symbol1绘制特征；如果特征的POPULATION属性值在 500,000 到 1,000,000 之间，则使用symbol2绘制特征；如果特征的POPULATION属性值是 1,000,000 或更多，则使用symbol3绘制特征。

注意

除了具有过滤器外，规则还可以有一个最小和最大比例因子，在这个比例因子下规则将生效。例如，当地图被放大到最大时，可以使用此功能隐藏较小的特征。

由于规则中可以有多个符号，特征的绘制方式也可以变得相当复杂。例如，您可以定义一个规则，该规则使用三个不同的符号化器来绘制 LineString 几何形状作为街道：

如您所想象，结合符号化、规则、过滤器和样式将为您在选择哪些特征应在地图中显示以及如何绘制这些特征方面提供极大的灵活性。

学习 Mapnik

现在您已经看到了 Mapnik 能做什么，并对 Mapnik 的工作方式有了一些了解，让我们更深入地了解 Mapnik 库的其他方面。在本章的这一节中，我们将介绍数据源、符号化和地图渲染。

数据源

每个地图层都与一个数据源（mapnik.Datasource的子类）相关联，该数据源提供要在地图上显示的数据。各种类型的数据源通过 C++插件提供，这些插件在编译 Mapnik 时启用或禁用。要检查是否安装了特定类型的数据源，请检查相关的插件是否已安装到您的 Mapnik 副本中。您可以通过在 Python 命令提示符中输入以下内容来查看已安装的插件列表（因此，支持的源数据）：

import mapnik
print list(mapnik.DatasourceCache.plugin_names())

Mapnik 当前支持以下数据源插件：

csv：此插件提供mapnik.CSV数据源，可以从文本文件或字符串中读取表格数据。默认情况下，数据以CSV（逗号分隔值）格式存储，尽管也支持其他类似格式。

CSV 数据源将自动根据包含“lat”、“latitude”、“lon”、“long”和“longitude”等名称的列头识别点几何形状。如果列头命名为“geojson”或“wkt”，数据源还将检测 GeoJSON 和 WKT 格式的几何形状。有关csv插件的文档，请参阅github.com/mapnik/mapnik/wiki/CSV-Plugin。
gdal: 此插件提供了 mapnik.Gdal 数据源。此数据源使用 GDAL 库读取栅格格式数据并将其提供给地图层。要在地图层中使用此数据源，您需要向地图层添加一个样式，该样式包含一个 mapnik.RasterSymbolizer 以在地图上绘制栅格数据。有关 gdal 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/GDAL。
ogr: 此插件实现了 mapnik.Ogr 数据源。此数据源使用 OGR 库读取矢量格式数据。有关 ogr 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/OGR。
osm: osm 插件提供了 mapnik.Osm 数据源。此数据源读取 OpenStreetMap XML 格式的数据。有关 osm 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/OsmPlugin。
postgis: 此插件提供了 mapnik.PostGIS 数据源。此数据源连接到 PostGIS 数据库，并从指定的数据库表中读取空间数据。在创建 PostGIS 数据源时，您使用 host、dbname、user 和 password 参数来告诉 Mapnik 如何连接到特定的 PostGIS 数据库，而 table 参数指定从哪个表中读取数据。

对于特殊目的，还提供了额外的参数，例如限制显示数据的范围，或使用 SQL 子查询仅包括数据库表中的某些记录。有关 postgis 插件的完整文档，请参阅 github.com/mapnik/mapnik/wiki/PostGIS。
python: 此插件提供了 mapnik.Python 数据源。这允许您通过编写一个自定义的 Python 类来访问要显示的数据，从而实现自己的数据源。要编写自定义 Python 数据源，您通常需要创建 mapnik.PythonDatasource 的子类，然后在调用 mapnik.Python() 函数实例化数据源时使用您自定义类的名称作为 factory 参数。然后，在您的类中实现必要的函数以提供对数据的访问。有关 python 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/Python-Plugin。
栅格：此插件实现了 mapnik.Raster 数据源，用于显示 TIFF 或 GeoTIFF 格式的栅格图像文件的内容。虽然您也可以使用 gdal 插件读取栅格格式数据，但在读取这些类型的文件时，raster 插件速度更快。要在地图层中使用此数据源，您需要向地图层添加一个样式，该样式包括一个 RasterSymbolizer 以在地图上绘制图像文件的内容。有关 raster 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/Raster。
形状：此插件提供了 mapnik.Shapefile 数据源，允许您读取形状文件。虽然 ogr 数据源也能读取形状文件，但通常使用 mapnik.Shapefile 数据源更方便。有关 shape 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/ShapeFile。
sqlite：此插件提供了 mapnik.SQLite 数据源。此数据源从 SQLite 数据库读取空间数据。数据库可以是包含 WKB 格式几何数据的普通 SQLite 数据库，也可以是使用 Spatialite 数据库扩展的空间数据库。有关 sqlite 插件的文档，请参阅 github.com/mapnik/mapnik/wiki/SQLite。

符号化器

符号化器负责在地图上绘制特征。通常使用多个符号化器来绘制单个特征——我们之前在同时使用 PolygonSymbolizer 绘制多边形内部和 LineSymbolizer 绘制多边形轮廓时看到了这一点。

Mapnik 内部有许多不同类型的符号化器可用，许多符号化器都与复杂的选项相关联。我们不会详尽地列出所有符号化器和它们的选项，而是仅查看一些更常见的符号化器类型及其使用方法。

点符号化

PointSymbolizer 类用于在点几何上绘制图像。默认情况下，每个点都显示为一个 4 x 4 像素的黑色方块：

要使用不同的图像，您必须创建一个 mapnik.PathExpression 对象来表示所需图像文件的路径，然后在实例化时将其传递给 PointSymbolizer 对象：

path = mapnik.PathExpression("/path/to/image.png")
point_symbol = PointSymbolizer(path)

注意，PointSymbolizer 在期望的点周围绘制图像。要使用前面示例中显示的推针图像，您需要添加额外的透明空白区域，以便推针的尖端位于图像的中间，如下所示：

你可以通过设置符号化的 opacity 属性来控制绘制图像的不透明度。你还可以通过设置 allow_overlap 属性为 True 来控制是否在图像上绘制标签。最后，你可以通过将 transform 属性设置为包含标准 SVG 变换表达式的字符串来对图像应用 SVG 变换，例如 point_symbol.transform = "rotate(45)"。

PointSymbolizer 的文档可以在 github.com/mapnik/mapnik/wiki/PointSymbolizer 找到。

LineSymbolizer

mapnik.LineSymbolizer 用于绘制 LineString 几何形状和 Polygon 几何形状的轮廓。当你创建一个新的 LineSymbolizer 时，通常会使用两个参数来配置它：用于绘制线的颜色，作为一个 mapnik.Color 对象，以及线的宽度，以像素为单位。例如：

line_symbol = mapnik.LineSymbolizer(mapnik.Color("black"), 0.5)

注意，你可以使用分数线宽；由于 Mapnik 使用抗锯齿，如果你在绘制很多紧密排列的线时，宽度小于 1 像素的线通常比整数宽度的线看起来更好。

除了颜色和宽度之外，你还可以通过设置 opacity 属性使线半透明。这应该设置为介于 0.0 和 1.0 之间的数字，其中 0.0 表示线将完全透明，1.0 表示线将完全不透明。

你还可以使用 stroke 属性来访问（或替换）由线符号化器使用的线对象。线对象是一个 mapnik.Stroke 的实例，可以用于更复杂的视觉效果。例如，你可以通过调用线的 add_dash() 方法来创建虚线效果：

line_symbol.stroke.add_dash(5, 7)

这两个数字都是以像素为单位的；第一个数字是虚线段的长度，而第二个数字是虚线之间的间隙长度。

提示

注意，你可以通过多次调用 add_dash() 来创建交替的虚线模式。

你还可以通过设置 stroke 的 line_cap 属性来控制线端应该如何绘制，以及通过设置 stroke 的 line_join 属性来控制当 LineString 改变方向时，各个线段之间的连接应该如何绘制。line_cap 属性可以设置为以下值之一：

mapnik.line_cap.BUTT_CAP
mapnik.line_cap.ROUND_CAP
mapnik.line_cap.SQUARE_CAP

line_join 属性可以设置为以下值之一：

mapnik.line_join.MITER_JOIN
mapnik.line_join.ROUND_JOIN
mapnik.line_join.BEVEL_JOIN

LineSymbolizer 类的文档可以在 github.com/mapnik/mapnik/wiki/LineSymbolizer 找到。

PolygonSymbolizer

mapnik.PolygonSymbolizer 类用于用给定的颜色填充 Polygon 几何形状的内部。当你创建一个新的 PolygonSymbolizer 时，通常会传递给它一个参数：用于填充多边形的 mapnik.Color 对象。你也可以通过设置 fill_opacity 属性来改变符号化的不透明度，例如：

fill_symbol.fill_opacity = 0.8

一次又一次，不透明度是从 0.0（完全透明）到 1.0（完全不透明）测量的。

除了 PolygonSymbolizer 属性外，还有一个你可能觉得有用的属性：gamma。gamma值可以设置为 0.0 到 1.0 之间的数字。gamma值控制绘制多边形边缘时使用的反走样量；默认的gamma值为1.0，多边形的边缘将完全反走样。虽然这通常是好事，但如果你尝试用相同的颜色绘制相邻的多边形，反走样会导致多边形的边缘可见，而不是将它们组合成更大的单一区域。通过稍微降低 gamma 值（例如，fill_symbol.gamma = 0.6），相邻多边形之间的边缘将消失。

PolygonSymbolizer类的文档可以在github.com/mapnik/mapnik/wiki/PolygonSymbolizer找到。

文本符号化

TextSymbolizer类用于在地图上绘制文本标签。此类符号化器可用于点、LineString 和 Polygon 几何形状。以下示例显示了如何使用 TextSymbolizer：

text_symbol = mapnik.TextSymbolizer(mapnik.Expresion("[label]"), "DejaVu Sans Book", 10, mapnik.Color("black"))

如你所见，通常将四个参数传递给 TextSymbolizer 的初始化器：

定义要显示的文本的mapnik.Expression对象。在这种情况下，要显示的文本将来自数据源中的label属性。
用于绘制文本的字体名称。要查看可用的字体，请在 Python 命令行中输入以下内容：
```
import mapnik
for font in mapnik.FontEngine.face_names():
    print font
```
字体大小，以像素为单位。
用于绘制文本的颜色。

默认情况下，文本将在几何形状的中心绘制；例如：

这种标签定位被称为点放置。TextSymbolizer 允许你将其更改为使用所谓的线放置，其中标签将沿着线条绘制：

text_symbol.label_placement = mapnik.label_placement.LINE_PLACEMENT

如你所见，这会导致标签沿着 LineString 几何形状的长度或多边形的周长绘制。对于点几何形状，不会绘制文本，因为点内没有线条。

TextSymbolizer 通常会只绘制一次标签，但你可以通过指定每个标签之间的像素间隔来告诉符号化器重复标签：

text_symbol.label_spacing = 30

默认情况下，Mapnik 足够智能，可以阻止标签相互重叠。如果可能，它会稍微移动标签以避免重叠，如果仍然会重叠，则会完全隐藏标签。例如：

你可以通过设置allow_overlap属性来更改此设置：

text_symbol.allow_overlap = True

最后，你可以设置一个光环效果，在文本周围绘制较浅颜色的边框，使其即使在深色背景上也能可见。例如，

text_symbol.halo_fill = mapnik.Color("white")
text_symbol.halo_radius = 1

有许多其他标签选项，所有这些都在TextSymbolizer类的文档中详细描述。这可以在github.com/mapnik/mapnik/wiki/TextSymbolizer找到。

RasterSymbolizer

RasterSymbolizer类用于在地图上绘制栅格格式数据。这种符号化类型通常与栅格或 GDAL 数据源一起使用。要创建一个新的栅格符号化，你需要实例化一个新的mapnik.RasterSymbolizer对象：

raster_symbol = mapnik.RasterSymbolizer()

栅格符号化将自动绘制地图层数据源提供的任何栅格格式数据。这通常用于绘制基图，在基图上显示矢量数据；例如：

虽然有一些高级选项可以控制栅格数据的显示方式，但在大多数情况下，你可能感兴趣的唯一选项是opacity属性。像往常一样，这设置了显示图像的不透明度，允许你将半透明的栅格图像一层层叠加。

RasterSymbolizer的文档可以在github.com/mapnik/mapnik/wiki/RasterSymbolizer找到。

地图渲染

我们现在已经详细检查了生成地图的许多构建块：图层、数据源、样式、规则、过滤器以及符号化。使用你所学的知识，你应该能够构建和设计你自己的地图。但是，一旦设置了mapnik.Map对象，你能用它做什么呢？

在本章开头我们检查的示例程序中，我们使用了mapnik.render_to_file()函数将生成的地图保存到图像文件中。在渲染地图时，你首先必须设置地图的范围——即定义地图可见部分的矩形：

只有地图的可见范围将被包含在生成的图像中；其他所有内容都将被忽略。

在我们的示例程序中，我们使用了map.zoom_all()来设置地图的可见范围，包括所有地图层中的所有要素。当然，有时你可能只想显示整体地图的一部分。为此，你可以使用map.zoomToBox()方法设置地图的可见范围。例如：

map.zoomToBox(mapnik.Box2d(-124.5, 32.0, -114.0, 43.0))

这四个数字分别代表最小经度、最小纬度、最大经度和最大纬度。如果你使用这些纬度和经度值执行此语句，地图的可见范围将覆盖大约美国的加利福尼亚州。

注意，你不仅限于只渲染一次地图。如果你想，你可以从单个mapnik.Map对象创建多个图像，改变可见范围，然后每次调用mapnik.render_to_file()将新可见的地图部分保存到不同的文件中。

一个工作示例

让我们将我们所学的所有内容结合起来，编写一个可以显示 shapefile 内容的程序。这是一个非常有用的程序，因为你可以操作或生成一些空间数据，将结果保存到 shapefile 中，然后运行此程序以显示生成的地图图像。

我们将我们的程序命名为 shapeToMap.py。创建此文件，并将以下 Python 代码开始输入到其中：

import mapnik

LAYERS = [
    {'shapefile'  : "TM_WORLD_BORDERS-0.3.shp",
     'lineColor'  : "black",
     'lineWidth'  : 0.4,
     'fillColor'  : "#709070",
     'labelField' : "NAME",
     'labelSize'  : 12,
     'labelColor' : "black"
    }
]

BACKGROUND_COLOR = "#a0c0ff"

BOUNDS_MIN_LAT  = 35.26
BOUNDS_MAX_LAT  = 71.39
BOUNDS_MIN_LONG = -10.90
BOUNDS_MAX_LONG = 41.13

MAX_WIDTH  = 1600
MAX_HEIGHT = 800

注意，我们在这里定义的各种常量将被用于配置我们将要生成的地图：

LAYERS: 这是一个要在地图上显示的地图层列表。列表中的每个项目应是一个字典，包含以下条目之一或全部：
- shapefile: 目标 shapefile 的名称和路径
- lineColor: 用于绘制特征外部的十六进制颜色代码（如果有的话）
- lineWidth: 用于绘制特征外部的线条宽度，以像素为单位
- fillColor: 用于绘制特征内部的十六进制颜色代码（如果有的话）
- labelField: 如果有的话，用于标记每个特征的源文件中的属性名称
- labelSize: 标记特征时使用的字体大小，以像素为单位
- labelColor: 用于绘制标签的十六进制颜色代码
BACKGROUND_COLOR: 这是用于绘制地图背景的十六进制颜色代码。
BOUNDS_MIN_LAT, BOUNDS_MIN_LONG, BOUNDS_MAX_LAT, 和 BOUNDS_MAX_LONG: 这些定义了你想要生成的地图的可视范围。
MAX_WIDTH 和 MAX_HEIGHT: 这些指定生成的地图图像的最大尺寸。请注意，生成的图像实际上可能小于这些值，具体取决于边界矩形的纵横比。

每次你想使用此程序生成地图时，你都需要编辑这些常量以适应你的需求。

我们接下来需要计算用于地图的高度和宽度。因为可视范围可以是任何形状，我们计算实际宽度和高度尽可能大，同时匹配可视范围的纵横比。我们通过首先计算地图的宽度和高度来实现这一点，使宽度为最大允许宽度，高度为匹配可视范围纵横比所需的任何值。为此，请将以下代码添加到程序末尾：

extent = mapnik.Envelope(BOUNDS_MIN_LONG, BOUNDS_MIN_LAT,  BOUNDS_MAX_LONG, BOUNDS_MAX_LAT)
aspectRatio = extent.width() / extent.height()

mapWidth = MAX_WIDTH
mapHeight = int(mapWidth / aspectRatio)

我们接下来查看计算出的高度是否太大，如果是，则缩小地图，使高度不超过允许的最大值：

if mapHeight > MAX_HEIGHT:
    scaleFactor = float(MAX_HEIGHT) / float(mapHeight)
    mapWidth = int(mapWidth * scaleFactor)
    mapHeight = int(mapHeight * scaleFactor)

这确保生成的地图尽可能大，同时确保地图具有与可视范围相同的纵横比。

现在我们知道了我们的地图有多大，我们可以创建并初始化我们的 mapnik.Map 对象：

map = mapnik.Map(mapWidth, mapHeight)
map.background = mapnik.Color(BACKGROUND_COLOR)

我们接下来需要定义我们各种地图样式，为每个地图层使用单个样式和规则。注意，我们使用 LAYERS 列表中的各种字典条目来为每个层定义地图样式：

for i,src in enumerate(LAYERS):
    style = mapnik.Style()
    rule = mapnik.Rule()

    if src['fillColor'] != None:
        symbol = mapnik.PolygonSymbolizer(
                    mapnik.Color(src['fillColor']))
        rule.symbols.append(symbol)
    if src['lineColor'] != None:
        symbol = mapnik.LineSymbolizer(
                    mapnik.Color(src['lineColor']),
                    src['lineWidth'])
        rule.symbols.append(symbol)
    if src['labelField'] != None:
        symbol = mapnik.TextSymbolizer(
                     mapnik.Expression(
                            "[" + src['labelField'] + "]"),
                            "DejaVu Sans Bold",
                            src['labelSize'],
                            mapnik.Color(src['labelColor']))
        symbol.allow_overlap = True
        rule.symbols.append(symbol)

    style.rules.append(rule)

    map.append_style("style-"+str(i+1), style)

我们现在需要定义我们地图的各种层：

for i,src in enumerate(LAYERS):
    layer = mapnik.Layer("layer-"+str(i+1))
    layer.datasource = mapnik.Shapefile(file=src['shapefile'])
    layer.styles.append("style-"+str(i+1))
    map.layers.append(layer)

最后，我们渲染地图图像：

map.zoom_to_box(extent)
mapnik.render_to_file(map, "map.png", "png")

由于你已经学习了本章“学习 Mapnik”部分的各种类和方法，你应该希望能够理解所有这些代码的作用。如果有什么不清楚的地方，请返回并复习该部分的相应部分。本章节的源代码中可以下载此程序的完整副本。

使用我们之前定义的各种常数，你应该能够使用这个程序来绘制世界边界数据集的内容。只需将TM_WORLD_BORDERS-0.3目录放置在与shapeToMap.py程序相同的文件夹中，然后尝试运行程序。如果一切顺利，程序应该生成一个map.png图像，该图像显示西欧和中欧的世界边界数据集内容：

如果你仔细观察这张图片，你会注意到一些标签被其他多边形遮挡。这是因为我们告诉程序在同一地图层中绘制多边形及其标签。为了解决这个问题，将你的LAYERS定义替换为以下内容：

LAYERS = [
    {'shapeFile'  : "TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp",
     'lineColor'  : "black",
     'lineWidth'  : 0.4,
     'fillColor'  : "#709070",
     'labelField' : None,
     'labelSize'  : None,
     'labelColor' : None,
    },	
    {'shapeFile'  : "TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp",
     'lineColor'  : None,
     'lineWidth'  : None,
     'fillColor'  : None,
     'labelField' : "NAME",
     'labelSize'  : 12,
     'labelColor' : "black"
    }
]

如你所见，我们现在正在两个独立的地图层中显示 shapefile，一个用于绘制国家多边形，另一个地图层用于在多边形前面绘制所有标签。如果你再次运行程序，你会看到标签问题已经得到解决。

这应该会让你了解shapeToMap.py程序有多有用。只需更改程序顶部的常数，你就可以快速查看任何 shapefile 的内容。实际上，本书中的许多插图都是使用此程序的修改版本生成的。

下一步

虽然shapeToMap.py程序被故意保持简单，以便更容易理解，但还有很多可以做的事情来改进这个程序并使其更有用。你可能喜欢通过实现以下新功能来挑战自己：

向每个地图层添加一个可选的labelHalo条目，用于在标签文本周围绘制光环。
向地图层添加一个labelPlacement条目，以便轻松控制标签放置选项。
向地图层添加一个labelAllowOverlap条目，用于控制标签是否允许重叠。
向地图层添加一个filter条目，用于构建一个mapnik.Filter()表达式，以限制在地图层中显示的特征集。
添加一个选项，根据每个特征的边界框动态计算地图的可视范围。这将允许你在事先计算边界之前生成地图。
添加一个调用os.system("open map.png")（适用于 Mac OS X）或os.startfile("map.png")（适用于 MS Windows）的调用，以便在生成图像后自动显示该图像。
添加对使用除默认 EPSG 4326 以外的投影的 shapefile 的支持。
从单独的模块中加载配置常量，这样你就不必每次想要更改要显示的数据时都编辑 Python 源文件。

作为本章源代码的一部分，提供了一个名为generateMap.py的shapeToMap.py的更复杂版本。generateMap.py程序实现了前面提到的所有建议。

摘要

在本章中，我们介绍了 Mapnik 地图生成库以及如何在 Python 程序中使用它来创建外观精美的地图。你安装了 Mapnik，查看了一个简单的示例，了解它如何使用，然后开始学习构建和样式化地图的过程。

我们随后对 Mapnik 进行了更详细的考察，研究了你可以用来将空间数据加载到地图层的各种数据源类型。我们还考察了可以用来显示空间特征的符号化工具，如何使用可见范围来控制要显示的地图部分，以及如何将地图渲染为图像文件。然后我们创建了一个有用的 Python 程序，名为shapeToMap.py，它可以用来从存储在 shapefiles 中的任何空间数据生成地图，并最终探讨了如何改进shapeToMap.py使其更加有用。

在下一章中，我们将探讨分析空间数据的各种工具和技术，包括如何使用 Python 解决各种有趣的地理空间问题。

第五章：分析地理空间数据

在本章中，我们将探讨分析地理空间数据的过程。有时，你的地理空间分析结果可能是一或多个数字，例如，有多少个国家位于赤道以南？ 或者 澳大利亚海岸线的长度是多少？ 在其他时候，你的分析结果可能是一个几何对象，例如，阿拉斯加的最北点在哪里？ 或者 德克萨斯州中哪些部分位于新墨西哥州以东？ 在其他时候，你的分析结果可能是一个列表，例如，哪些国家位于新几内亚 1,000 公里范围内？ 在所有这些情况下，你需要熟悉可用于计算所需结果的工具和技术。

为了帮助你学习这些工具和技术，我们将检查以下内容：

如何安装和使用两个强大的 Python 库来解决地理空间问题
如何计算和操作位置
如何计算 LineString 几何形状在现实世界单位中的长度
如何使用现实世界单位计算多边形的面积
如何使用包含道路的 shapefile 构建连接的 LineString 几何形状的抽象模型，然后使用该模型计算两点之间的最短路径。

让我们从查看一些你可以用于地理空间分析的 Python 库开始。

空间分析库

你已经有一些对分析地理空间数据有用的库：OGR 库包括比较和操作几何形状的方法，Shapely 是一个用于处理和分析几何数据的优秀库。然而，还有两个其他库你可能想熟悉：PyProj，这是一个用于在地球表面上计算距离和位置的强大库，以及NetworkX，它可以从地理空间数据中构建抽象的数学模型，然后分析这些模型以解决各种问题。

让我们更详细地研究这两个库，并将它们都安装到你的计算机上。

PyProj

PyProj (pypi.python.org/pypi/pyproj) 是一个使用 Python 处理空间参考系统的强大工具。PyProj 本身只是 PROJ.4 地图投影库的 Python 接口，该库是用 C 编写的。因此，要安装 PyProj，你通常需要安装 PROJ.4 库，然后安装或构建 PyProj 本身。

在我们深入了解安装 PyProj（和 PROJ.4）的细节之前，让我们看看这个库的功能以及它如何有用。如果你还记得第二章中的内容，“地理空间数据”，空间参考系统是使用坐标表示地球表面位置的一种方式。未投影坐标，如纬度和经度值，通过从地球中心到所需点的线条追踪，然后测量该线条在东西和南北方向上的角度，直接表示地球表面的一个位置：

另一方面，投影坐标表示为二维笛卡尔平面上的位置：

空间参考系统，也称为地图投影，是将地球表面的点转换为二维平面上的点的一种方式。PROJ.4（以及 PyProj）是处理这些投影的工具。

了解用于生成你的地理空间数据的地图投影至关重要。使用错误的投影将毁掉你所有的计算和地图可视化。地图投影也很重要，因为使用未投影坐标系统（如 EPSG 4326）中的数据进行空间计算几乎是不可能的。例如，想象一下你想计算以下多边形的面积，它代表苏格兰尼斯湖的轮廓：

这个多边形的坐标是 EPSG 4326——也就是说，它们是纬度和经度值。如果你愿意，你可以将这个多边形加载到 Shapely 中，并要求它计算面积：

>>> wkt = "POLYGON ((-4.335556 57.373056,-4.676389 57.138611,-4.447778 57.324722,-4.349167 57.386944,-4.334444 57.379444,-4.335556 57.373056))"
>>> polygon = shapely.wkt.loads(wkt)
>>> print polygon.area
0.006077434151

然而，结果是按度数计算的面积，这是一个没有意义的数字。这是因为 Shapely 不了解空间参考系统。它天真地对待纬度和经度值作为 (x,y) 坐标，这意味着这种空间计算无法产生有用或准确的结果。

你真正想要的是以有意义的单位（如平方米或平方英里）测量的面积。这正是 PyProj 发挥作用的地方。PyProj 允许你使用任何空间参考系统进行计算和转换。PyProj 承担了所有的繁重数学运算，因此你不必亲自去做。

现在，让我们安装 PyProj 并看看它是如何工作的。你安装 PyProj（以及底层的 PROJ.4 库）的方式取决于你使用的操作系统：

对于 MS Windows，你可以从 www.lfd.uci.edu/~gohlke/pythonlibs/#pyproj 安装预构建的 PyProj 版本。此安装程序包括 PROJ.4，因此你不需要单独安装它。
对于 Mac OS X，你需要执行以下操作：
1. 下载并安装 PROJ.4 库。可以从 www.kyngchaos.com/software/frameworks 下载 PROJ.4 的 Mac 安装程序。
2. 如果你电脑上还没有安装 XCode，请前往 Mac App Store 下载最新版本。XCode 是苹果的开发系统，可以免费下载。
3. 如果你使用的是低于 10.9（Yosemite）版本的 Mac OS X，你需要单独安装命令行工具。为此，启动 XCode 并从 XCode 菜单中选择 Preferences… 命令。在 Downloads 选项卡中，将有一个安装命令行工具的选项；启用此选项并等待所需工具安装完成。
  
  提示
  
  如果你使用的是 Mac OS X 10.9（Yosemite）或更高版本，你可以跳过此步骤。
4. 从 pypi.python.org 下载 PyProj 的源代码。
5. 使用终端，cd 到你下载的 PyProj 目录，并输入以下命令：
```
python setup.py build
sudo python.setup.py install
```
6. 最后，启动 Python 并尝试输入以下命令：
```
import pyproj
```
  Python 提示符应该会重新出现，没有任何错误信息。
对于 Linux，你可以使用你喜欢的包管理器安装 PROJ.4 和 PyProj，或者你可以按照 trac.osgeo.org/proj 和 github.com/jswhit/pyproj 上提供的说明从源代码构建它们。

现在你已经安装了 PyProj，让我们看看这个库是如何使用的。PyProj 提供了两个主要的类：

Proj：这个类代表一个空间投影，并允许你将坐标在投影之间进行单次或批量转换
Geod：这是一个 大地测量 计算类，允许你根据使用给定空间参考系统的坐标执行各种计算

让我们看看如何使用 PyProj 来计算两点之间的距离。打开一个终端窗口，启动你的 Python 解释器，并输入以下代码：

import pyproj
geod = pyproj.Geod(ellps="WGS84")
lat1 = -38.137
long1 = 176.349
lat2 = -41.286
long2 = 174.776
heading1,heading2,distance = geod.inv(long1, lat1, long2, lat2)
print int(distance)

这两个坐标代表新西兰罗托鲁阿和惠灵顿两座城市的地理位置：

如果一切顺利，你应该会看到打印出数字 374,729，这是这两座城市之间的直线距离（大圆距离），单位是米。

注意

注意，我们使用 ellps="WGS84" 来设置 Geod 对象的空间参考系统。此值设置了用于大地测量计算器地球形状的数学模型——WGS84 是 EPSG 4326 空间参考系统使用的椭球体的名称，因此我们实际上是在告诉 PyProj 坐标是以纬度和经度测量的。

PyProj 还可以用于在坐标系之间进行转换。我们将在看到它如何与 Shapely 结合以执行精确的空间计算时简要介绍这一点。

NetworkX

NetworkX 是一个用于定义和分析数学图的 Python 库。在数学术语中，图是由顶点通过边连接在一起的一组：

每条边通常都会分配一个值，称为权重，它可以用来对图进行查询。每条边可以是有向的——也就是说，你只能沿着一个方向跟随边——或者它可以是无向的，允许你沿着任意方向跟随边。

虽然这些图是一个有趣的数学概念，但它们在地理空间分析中也非常有用，因为你可以使用图来表示地点以及在这些地点之间移动的方式。例如，顶点可能代表城镇，而边可能代表连接这些城镇的道路。当以这种方式使用时，边通常根据道路的长度进行加权，这样较长的道路就有更大的权重，整个图可以用来计算两点之间的最短路径。

NetworkX 是一个用于处理数学图非常强大的库。更好的是，它包括读取 shapefile 并将其转换为图的能力。这允许你非常简单地转换地理空间数据为抽象图表示，然后你可以用它以各种有用的方式分析地点之间的关系。

注意

在 NetworkX 中，图中的顶点被称为节点。

让我们继续安装 NetworkX 库。NetworkX 的主网站可以在networkx.github.io找到，你可以直接从pypi.python.org/pypi/networkx下载库。

由于 NetworkX 是用纯 Python 编写的，你可以简单地下载源代码，然后输入python setup.py install将其安装到你的site-packages目录中，或者如果你更喜欢，你可以通过输入pip install networkx使用 pip 安装它。

对于更多详细信息，请参阅 NetworkX 的安装说明，可以在networkx.github.io/documentation/latest/install.html找到。

一旦你安装了 NetworkX，可以通过在 Python 命令提示符中输入以下内容来检查它是否工作：

import networkx

graph = networkx.Graph()
graph.add_edge("New York", "San Francisco", weight=2908)
graph.add_edge("San Francisco", "Los Angeles", weight=382)
graph.add_edge("Los Angeles", "New York", weight=2776)
graph.add_edge("San Francisco", "Portland", weight=635)

print networkx.shortest_path(graph, "New York", "Portland")

这个简单的程序构建了一个 NetworkX 图，其中的节点代表城市，而边代表连接这些城市的道路。对于每条边，权重代表这两个城市之间的驾驶距离（英里）。使用这个简单的图，我们然后要求 NetworkX 显示纽约和俄勒冈州波特兰之间的最短路径。

如果一切顺利，运行前面的代码将告诉你最短路径是从纽约到旧金山，然后从那里到波特兰：

['New York', 'San Francisco', 'Portland']

显然，你可以用 NetworkX 做更多的事情，我们将在本章后面使用这个库。然而，目前来说，了解你可以使用 NetworkX 从你的空间数据中构建一个抽象图，然后使用 NetworkX 的分析函数根据你的图计算有用的信息就足够了。

空间分析食谱

现在我们已经拥有了一套完整的地理空间分析库，让我们看看我们如何使用它们来解决一些实际问题。我们将探讨如何计算和操作位置、长度和面积，以及如何使用 NetworkX 来计算两点之间的最短路径。

计算和比较坐标

如我们之前所见，PyProj 可以用来计算两个位置之间的实际距离。它还可以用来测量两点之间线的角度，并基于起点、距离和航向计算新点。

让我们使用 PyProj 来计算两点之间的距离。然后我们将用它来计算从给定点出发的一定距离和航向的位置。

首先创建一个新的 Python 程序，命名为 coord_analysis.py。将以下代码输入到这个程序的开始部分：

import pyproj
geod = pyproj.Geod(ellps="WGS84")

到目前为止，这与我们之前看到的代码相同：我们只是导入 PyProj 库并基于 WGS84 椭球体创建一个大地测量计算对象。如果你还记得，这是标准 EPSG 4326 空间参考系统所使用的地球表面的数学模型。

现在我们将添加一些代码来提示用户输入所需的坐标。这完全是标准的 Python 代码，不需要任何进一步的解释：

def get_coord(prompt):
    while True:	
        s = raw_input(prompt + " (lat,long): ")
        if "," not in s: continue
        s1,s2 = s.split(",", 1)
        try:	
            latitude = float(s1.strip())
        except ValueError:
            continue
        try:
            longitude = float(s2.strip())
        except ValueError:
            continue
        return latitude,longitude
lat1,long1 = get_coord("Starting coordinate")
lat2,long2 = get_coord("Ending coordinate")

现在我们有了两组纬度和经度值，我们可以使用 PyProj 来计算这两点之间的实际距离：

heading1,heading2,distance = geod.inv(long1, lat1, long2, lat2)

这正是我们之前看到的相同代码。geod.inv() 方法接受两个坐标，并返回航向（从第一个点到第二个点的线角度，从正北顺时针测量），逆航向（从第二个点回到第一个点的线角度，同样从正北顺时针测量），以及两点之间的距离（以米为单位）。

注意，调用 geod.inv() 需要我们在纬度值之前提供经度值。这是因为 PyProj 与任何坐标系统一起工作，而经度代表 x（从左到右）值，而纬度代表 y（从下到上）值。在处理可能位于任何空间参考系统中的通用坐标时，x 值始终列在第一位。

现在我们已经计算了这三个数字，让我们显示它们，以便用户可以看到我们计算的结果：

print "Heading = %0.2f degrees" % heading1
print "Inverse heading = %0.2f degrees" % heading2
print "Distance = %0.2f kilometers" % (distance/1000)

要查看结果，请保存你的程序并在终端窗口中运行它。然后输入以下坐标到你的程序中：

Starting coordinate: 37.774929, -122.419416
Ending coordinate: 34.052234, -118.243685

这两个坐标代表旧金山和洛杉矶的位置。假设你已正确输入程序，以下结果应该会显示：

Heading = 136.38 degrees
Inverse heading = -41.17 degrees
Distance = 559.04 kilometers

这告诉我们，如果我们直接在旧金山市中心上空的飞机上飞行 559 公里，航向为 136.38 度（从正北顺时针测量），你将最终到达洛杉矶市中心。同样，如果你在洛杉矶，航向为-41.17 度（再次从正北顺时针测量），你将最终到达旧金山：

现在，让我们添加一些代码来计算从现有位置出发，一定距离和航向的点的坐标。注释掉get_coord()函数之后的全部内容，然后添加以下内容到程序末尾：

def get_num(prompt):
    while True:
        s = raw_input(prompt + ": ")
        try:
            value = float(s)
        except ValueError:
            continue
        return value

这是一个简单的实用函数，用于提示用户输入数值。我们将使用这个（以及我们之前编写的get_coord()函数）来提示用户输入我们需要的信息。现在将以下内容添加到程序末尾：

sLat,sLong = get_coord("Starting coordinate")
distance = get_num("Distance in kilometers") * 1000
heading = get_num("Heading")

注意，我们将以公里为单位的距离转换为以米为单位的距离——PyProj 始终以米为单位工作，因此我们必须以米为单位提供距离。

我们现在可以计算终点坐标了。使用 PyProj，这很简单：

eLong,eLat,iHeading = geod.fwd(sLong, sLat, heading, distance)

geod.fwd()方法返回所需的坐标（X 值列在前面），以及反向航向。我们的最后一个任务是向用户显示这些结果：

print "End point = (%0.4f,%0.4f)" % (eLat, eLong)
print "Inverse heading = %0.2f degrees" % iHeading

如果你运行这个程序，你可以尝试输入起点、航向和距离，程序将显示终点。例如：

Starting coordinate (lat,long): 37.774929, -122.419416
Distance in kilometers: 559.04
Heading: 136.38
End point = (34.0521,-118.2440)
Inverse heading = -41.17 degrees

注意

计算出的终点坐标与之前看到的不完全相同。这是因为距离和航向只指定到小数点后两位的精度。

当然，对于我们的程序，我们注释掉了第一次计算，以便我们能够专注于第二次计算。一个明显的改进是添加一个简单的文本提示，询问用户要执行哪个计算。但你可以看到 PyProj 是如何用来计算和比较地球表面上的点的——当你使用经纬度值时，这不容易做到。

计算长度

现在我们知道了如何计算两点之间的距离（以米为单位），让我们应用这项技术来计算任何 LineString 几何形状的真实长度。

要计算 LineString 几何形状的总长度，我们需要将 LineString 分割成单个线段，计算每个线段的长度，然后将结果相加以获得整个 LineString 的总长度：

要看到这个动作，我们需要一些 LineString 几何来工作。在这个例子中，我们将使用代表加利福尼亚州主要和次要道路的 LineStrings。这些数据可以从美国人口普查局的网站www.census.gov/geo/maps-data/data/tiger-line.html下载。向下滚动到标记为2014 TIGER/Line Shapefiles的部分，然后点击下载选项，然后点击Web interface。在下载页面，从Select a layer type下拉菜单中选择Roads，然后点击Submit按钮。最后，从Primary and Secondary Roads下拉菜单中选择California，然后点击Download按钮以下载所需数据。

生成的 shapefile 将是一个名为tl_2014_06_prisecroads.zip的压缩 ZIP 存档。解压这个存档，并将生成的 shapefile 放置在方便的位置。在 shapefile 所在的目录中创建一个名为calc_lengths.py的新文件，然后将以下代码输入到这个文件中：

import osgeo.ogr
import shapely.wkt
import pyproj

geod = pyproj.Geod(ellps="WGS84")

shapefile = osgeo.ogr.Open("tl_2014_06_prisecroads.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = shapely.wkt.loads(
        feature.GetGeometryRef().ExportToWkt())

这一切对你来说都应该非常熟悉——我们只是导入我们将要使用的各种库，创建一个新的pyproj.Geod对象用于我们的长度计算，然后逐个遍历 shapefile 的内容。正如你所见，我们使用shapely.wkt.loads()函数将特征转换为 Shapely 几何对象。

现在我们有了 Shapely 几何，我们的下一个任务是将其分割成单独的线段并计算每个线段的长度。让我们来做这件事：

    tot_length = 0
    prev_long,prev_lat = geometry.coords[0]
    for cur_long,cur_lat in geometry.coords[1:]:
        heading1,heading2,distance = geod.inv(
            prev_long, prev_lat, cur_long, cur_lat)
        tot_length = tot_length + distance
        prev_long,prev_lat = cur_long,cur_lat

由于 Shapely 几何是 LineString，我们可以使用geometry.coords访问构成 LineString 的各个坐标。然后我们依次处理每一对坐标，使用我们之前学到的技术来计算两个坐标之间以米为单位的距离。我们跟踪所有坐标对的总计算长度，从而得到 LineString 几何的总长度。

小贴士

如果我们的道路数据是在一个保持距离的投影坐标系中，我们可以简单地要求 Shapely 通过访问geometry.length属性来计算每个 LineString 的总长度。然而，这对于 EPSG 4326 数据是不行的，因为结果又将是度数长度。

我们最后一个任务是处理计算出的长度。让我们简单地打印出来，包括道路的名称：

    print feature.GetField("FULLNAME"), int(tot_length)

理论上，我们的程序现在应该可以工作了，让我们尝试运行它：

$ python calc_lengths.py
N Wheeler Ridge Rd 1616
N Main St 1595
...

到目前为止一切顺利；正如你所见，每条道路的总长度（以米为单位）正在被计算。不幸的是，如果我们再等几秒钟，我们的程序将引发 Python 异常并停止：

Traceback (most recent call last):
 File "calc_lengths.py", line 23, in <module>
 prev_long,prev_lat = geometry.coords[0]
 File "/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7/shapely/geometry/base.py", line 634, in coords
 "Multi-part geometries do not provide a coordinate sequence")
NotImplementedError: Multi-part geometries do not provide a coordinate sequence

这里发生了什么？看起来geometry.coords不可用，因为几何形状不是一个普通的 LineString。实际上，如果您还记得第二章中的内容，shapefiles 在简单几何形状和这些几何形状的集合之间没有区别，因此 shapefile 中的 LineString 实际上可能是一系列 LineString 的集合。这正是本例中发生的情况——如果您使用交互式 Python 命令提示符将受影响的要素加载到内存中，您可以打印出几何类型以查看问题所在：

>>> geometry = shapely.wkt.loads(feature.GetGeometryRef().ExportToWkt())
>>> print geometry.geom_type
MultiLineString

因此，我们用 MultiLineString 几何形状而不是 LineString 几何形状来表示道路。幸运的是，很容易将 MultiLineString 拆分并逐个处理单个 LineString。以下是添加此功能后我们整个程序的外观：

import osgeo.ogr
import shapely.wkt
import pyproj

geod = pyproj.Geod(ellps="WGS84")

shapefile = osgeo.ogr.Open("tl_2014_06_prisecroads.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = shapely.wkt.loads(
        feature.GetGeometryRef().ExportToWkt())

    lineStrings = []
    if geometry.geom_type == "LineString":
        lineStrings.append(geometry)
    elif geometry.geom_type == "MultiLineString":
        for lineString in geometry:
            lineStrings.append(lineString)

    tot_length = 0

    for lineString in lineStrings:
        prev_long,prev_lat = lineString.coords[0]
        for cur_long,cur_lat in lineString.coords[1:]:
            heading1,heading2,distance = geod.inv(
                prev_long, prev_lat, cur_long, cur_lat)
            tot_length = tot_length + distance
            prev_long,prev_lat = cur_long,cur_lat

    print feature.GetField("FULLNAME"), int(tot_length)

使用这种技术，我们可以计算任何线性几何形状（如 LineString 或 MultiLineString）的确切长度，单位为米。实际上，我们甚至可以使用相同的技巧来计算多边形几何形状的周长，通过访问多边形的外部线性环，然后将其处理成 LineString：

lineString = polygon.exterior

计算面积

在第二章中，我们了解到如何使用 OGR 和世界莫勒韦德投影（EPSG 54009）来计算多边形的面积。世界莫勒韦德是一种等面积地图投影，在全球范围内具有较高的准确性，因此适用于计算平方米面积。

不使用 OGR，让我们使用 Shapely 应用相同的技巧。这样做的好处是，我们将能够访问 Shapely 的所有功能，允许我们以各种有用的方式操作和测量几何形状。为此，我们将使用一个方便的 Shapely 函数shapely.ops.transform()。这个函数允许您将转换函数应用于几何形状中的每个坐标。转换可以是任何您想要的内容（如果您想的话，您可以在 Python 中编写自己的转换函数），但最重要的是，您可以使用 PyProj 实现一个转换函数，将 EPSG 4326 转换为 EPSG 54009，以便您能够准确计算任何几何形状的面积。

让我们看看它是如何工作的。将您之前下载的TM_WORLD_BORDERS-0.3 shapefile 的副本放入一个方便的目录中，并在同一目录中创建一个名为calc_areas.py的新文件。然后，将以下代码输入到这个新文件中：

import osgeo.ogr
import shapely.wkt
import shapely.ops
import pyproj

shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)

src_proj = pyproj.Proj(proj="longlat", ellps="WGS84", datum="WGS84")
dst_proj = pyproj.Proj(proj="moll", lon_0=0, x_0=0, y_0=0, ellps="WGS84", datum="WGS84", units="m")

def latlong_to_mollweide(longitude, latitude):
    return pyproj.transform(src_proj, dst_proj,
                            longitude, latitude)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    wkt = feature.getGeometryRef().ExportToWkt()
    geometry = shapely.wkt.loads(wkt)

    transformed = shapely.ops.transform(latlong_to_mollweide,
                                        geometry)
    area = int(transformed.area/1000000)

    print feature.GetField("NAME"), area

如您所见，我们定义了一个转换函数latlong_to_mollweide()，它将给定的纬度和经度值转换为莫勒韦德坐标系中的(x,y)坐标。莫勒韦德基于米，因此 Shapely 可以针对该几何形状进行计算，并以米为单位返回结果。

当你运行 calc_areas.py 程序时，你应该会看到一个国家列表以及与 World Borders Dataset 中相关多边形的面积，以平方公里为单位：

$ python calc_areas.py
Antigua and Barbuda 546
Algeria 2326137
Azerbaijan 86014
Albania 28702
...

使用 shapely.ops.transform() 的好处是，你可以对生成的几何形状使用 Shapely 的所有计算和几何操作功能。例如，新西兰有一个专属经济区，从海岸线向外延伸 200 英里。使用 buffer() 方法，你可以通过将新西兰的轮廓扩展到包括海岸线 200 英里内的所有点来计算这个专属经济区的形状：

注意

缓冲区是一个极其强大的操作。例如，你可以使用 buffer() 与 intersects() 方法一起识别给定起始几何形状一定距离内的所有国家。例如：

buffered_area = test_country['geom'].buffer(1000000)
for country in countries:
    if country['geom'].intersects(buffered_area):
        print "%s is within 1000 km of %s" % (country['name'], test_country['name'])

计算最短路径

在我们的最终示例中，我们将使用包含道路数据的 shapefile 来计算两点之间的最短路径。这是一个相当复杂的示例，使用了各种分析和操作几何数据的技术。它还使用了 NetworkX 库在道路网络的抽象表示上执行最短路径计算。

让我们从看看 NetworkX 如何将包含 LineString 几何形状的 shapefile 转换为抽象网络开始。如果你查看 tl_2014_06_prisecroads shapefile 的一小部分，你会看到看起来像是一系列相连的道路，例如：

然而，道路实际上并没有在它们相交的地方停止——道路特征只是继续延伸，根据需要与其他道路重叠。在地图上，这些可能看起来像是交叉点，但实际上在两条道路相遇或交叉的地方没有真正的交叉点：

这很重要，因为 NetworkX 将 LineString 几何形状转换为抽象图的方式——NetworkX 会认为只有当两个 LineString 具有相同的起始或结束点时它们才会相交；简单的交叉并不意味着两条道路相交。在先前的例子中，道路 2 和 道路 4 将是唯一被认为是相连的道路——即使 道路 2 看起来与 道路 1 和 道路 3 相交，由于没有匹配的端点，这些道路将被排除在图之外。

为了允许 NetworkX 将道路 shapefile 转换为网络，我们需要在道路相交的点处分割道路：

从数学的角度讲，这被称为从重叠道路网络中生成平面图。这个过程并不完美——它忽略了桥梁和立交桥，以及禁止进入标志阻止旅行者采取特定转弯的地方。然而，将道路形状文件转换为平面图是一个好的起点，如果你有一份桥梁和其他禁止进入点的列表，你总是可以排除特定的交叉口来计算。

让我们继续将我们的道路形状文件转换为平面图。为此，创建一个新的 Python 程序，命名为split_roads.py，并将以下代码输入到这个文件中：

import os
import os.path
import shutil
import osgeo.ogr
import osgeo.osr
import shapely.wkt

SRC_SHAPEFILE = "tl_2014_06_prisecroads.shp"

all_roads = []
shapefile = osgeo.ogr.Open(SRC_SHAPEFILE)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    wkt = feature.GetGeometryRef().ExportToWkt()
    geometry = shapely.wkt.loads(wkt)
    all_roads.append(geometry)

除了额外的import语句（我们很快就会需要它们）之外，这段代码应该相当清晰：我们只是将tl_2014_06_prisecroads.shp形状文件中的 LineString 几何形状加载到all_roads列表中。

提示

如果你的形状文件在不同的目录中，编辑SRC_SHAPEFILE常量，以便程序可以找到形状文件。

我们接下来的任务是分割交叉口处的道路。幸运的是，Shapely 使这变得相当简单。将以下内容添加到程序的末尾：

split_roads = []

for i in range(len(all_roads)):
    cur_road = all_roads[i]
    crossroads = []
    for j in range(len(all_roads)):
        if i == j: continue
        other_road = all_roads[j]
        if cur_road.crosses(other_road):
            crossroads.append(other_road)
    if len(crossroads) > 0:
        for other_road in crossroads:
            cur_road = cur_road.difference(other_road)
        if cur_road.geom_type == "MultiLineString":
            for split_road in cur_road.geoms:
                split_roads.append(split_road)
        elif cur_road.geom_type == "LineString":
            split_roads.append(cur_road)
    else:
        split_roads.append(cur_road)

正如你所见，我们通过调用 Shapely 的crosses()方法来识别与当前道路交叉的任何道路。然后，我们使用difference()方法从当前道路中移除每个交叉口；这会在其他道路与之交叉的点将道路分割开来：

最后，我们想要将分割的道路保存回形状文件。为此，将以下代码添加到程序的末尾：

driver = osgeo.ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists("split_roads"):
    shutil.rmtree("split_roads")
os.mkdir("split_roads")
dstFile = driver.CreateDataSource("split_roads/split_roads.shp")

spatialReference = osgeo.osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("WGS84")

layer = dstFile.CreateLayer("Layer", spatialReference)

for road in split_roads:
    wkt = shapely.wkt.dumps(road)
    linestring = osgeo.ogr.CreateGeometryFromWkt(wkt)

    feature = osgeo.ogr.Feature(layer.GetLayerDefn())
    feature.SetGeometry(linestring)

    layer.CreateFeature(feature)
    feature.Destroy()

dstFile.Destroy()

这段代码你应该很熟悉，因为我们已经在第二章中使用了相同的技巧，当我们查看如何将矢量数据写入形状文件时。这里唯一的新事物是我们使用shutil.rmtree()后跟os.mkdir()来删除并重新创建存储形状文件的目录；这允许我们再次运行程序，而无需每次都记得删除形状文件。

这就完成了我们的split_roads.py程序。将tl_2014_06_prisecroads形状文件中的 8,000 条道路分割开来可能需要几分钟时间，所以当你开始编写下一个程序时，只需让它运行即可。

一旦我们得到了分割道路的集合，我们就会想要另一个程序来使用它们来计算两点之间的最短路径。现在让我们开始编写这个程序。我们将把这个程序命名为calc_shortest_path.py。创建这个文件，并将以下代码输入到其中：

import shapely.wkt
import pyproj
import networkx

现在，我们将编写一些实用函数，这些函数是我们进行最短路径计算所需的。首先，我们将使用我们之前看到的技巧来计算两点之间以米为单位的距离：

def calc_distance(lat1, long1, lat2, long2):
    geod = pyproj.Geod(ellps="WGS84")
    heading1,heading2,distance = geod.inv(long1, lat1, long2, lat2)
    return distance

我们将使用这个函数来编写另一个函数，用于计算 LineString 几何形状的总长度：

def calc_length(linestring):
    tot_length = 0
    prev_long,prev_lat = linestring.coords[0]
    for cur_long,cur_lat in linestring.coords[1:]:
        distance = calc_distance(prev_lat, prev_long,
                                 cur_lat, cur_long)
        tot_length = tot_length + distance
        prev_long,prev_lat = cur_long,cur_lat
    return int(tot_length)

接下来，我们需要复制我们之前编写的get_coord()函数：

def get_coord(prompt):
    while True:
        s = raw_input(prompt + " (lat,long): ")
        if "," not in s: continue
        s1,s2 = s.split(",", 1)
        try:
            latitude = float(s1.strip())
        except ValueError:
            continue
        try:
            longitude = float(s2.strip())
        except ValueError:
            continue
        return latitude,longitude

我们还需要编写一个额外的函数：find_closest_node。这个函数将在 NetworkX 图中找到最接近给定纬度和经度值的节点。我们需要这个函数来识别最短路径计算中的起始和结束节点。

下面是find_closest_node函数的代码，你应该将其添加到程序的末尾：

def find_closest_node(graph, latitude, longitude):
    closest_node = None
    min_distance = None
    for node in graph.nodes():
        distance = calc_distance(node[1], node[0],
                                 latitude, longitude)
        if closest_node == None:
            closest_node = node
            min_distance = distance
        elif distance < min_distance:
            closest_node = node
            min_distance = distance
    return closest_node

要找到最近的节点，我们只需遍历图中的所有节点（顶点），并计算节点与所需坐标之间的距离（以米为单位）。然后我们选择并返回距离最小的节点。

现在我们已经准备好开始编写程序的主要部分。第一步是让 NetworkX 读取split_roads形状文件，并从道路数据中创建一个图：

graph = networkx.read_shp("split_roads/split_roads.shp")

这将读取形状文件并生成一个 NetworkX 图，其中每条边代表一条道路，每个节点代表道路的端点。因为 NetworkX 无法知道如何识别道路或端点，它使用纬度和经度来识别每个端点（即每个节点），以及起始和结束的纬度和经度来识别每条道路（即每条边）。因此，生成的图将包含类似于以下图示的节点和边：

结果图将会相当大，因为我们几乎需要从我们的 shapefile 中导入 10,000 条道路。

我们接下来的任务可能看起来有些奇怪：因为没有保证每条道路都能从其他每条道路到达，我们需要将图缩减到仅包含可到达的道路集合。如果我们不这样做，我们的最短路径计算很可能会失败。为了移除不可到达的道路，我们使用connected_component_subgraphs()函数来识别包含最多连接道路的图的一部分，并使用这个子图进行最短路径计算。以下是必要的代码：

graph = networkx.connected_component_subgraphs(graph.to_undirected()).next()

注意，因为connected_component_subgraphs()函数需要一个无向图，而read_shp()函数返回一个有向图，所以我们必须使用to_undirected()方法使图变为无向。

提示

如果你遇到'list' object has no attribute 'next'错误，你可能在使用 NetworkX 的不同版本。在这种情况下，将此行替换为graph = networkx.connected_component_subgraphs(graph.to_undirected())[0]。

现在我们已经有了可用的道路集合，我们的下一个任务是计算每条道路的长度。这个长度值将作为最短路径计算的基础。幸运的是，长度计算相当直接：

for node1,node2 in graph.edges():
    wkt = graph[node1][node2]['Wkt']
    linestring = shapely.wkt.loads(wkt)
    length = calc_length(linestring)
    graph.edge[node1][node2]['length'] = length

如您所见，我们提取了每个边的原始 LineString 几何形状，以 WKT 格式表示，然后使用它来创建一个 Shapely 几何对象。然后我们使用我们的calc_length()函数来计算道路的总长度，并将结果值作为length属性存储到边中。运行此代码将计算并存储图中每条道路的长度。

完成这些后，我们终于可以计算最短路径了。我们首先要求用户输入所需起始点和终点的大致纬度和经度值：

start_lat, start_long = get_coord("Starting Coordinate")
end_lat, end_long = get_coord("Ending Coordinate")

用户输入的值定义了两个坐标；我们需要使用这些坐标来识别起始节点和结束节点。我们可以使用我们之前编写的find_closest_node()函数来完成这项工作：

start_node = find_closest_node(graph, start_lat, start_long)
end_node   = find_closest_node(graph, end_lat, end_long)

现在我们可以根据之前计算出的长度值，获取两个节点之间的最短路径：

path = networkx.shortest_path(graph, start_node, end_node, "length")

返回的path值是由构成最短路径的节点组成的列表。让我们通过打印出这条路径的详细信息来完成我们的程序：

tot_length = 0
prev_node = path[0]
for cur_node in path[1:]:
    edge = graph.edge[prev_node][cur_node]
    print (str(prev_node) + " -> " + str(cur_node) +
           ", length = " + str(edge['length']))
    tot_length = tot_length + edge['length']
    prev_node = cur_node
print "Total length = " + str(tot_length)

现在我们已经完成了程序，让我们来测试一下。运行calc_shortest_path.py脚本。程序将首先将道路网络加载到内存中，并计算每条道路的长度：

$ python calc_shortest_path.py 
Loading road network into memory...
graph has 7976 nodes and 9709 edges
Calculating road lengths...

在计算长度之后，程序将提示您输入所需的起始和结束坐标。让我们输入奥克兰和圣路易斯奥比斯波这两个位于加利福尼亚州的城市坐标：

Starting Coordinate (lat,long): 37.794189, -122.276469
Ending Coordinate (lat,long): 35.281107, -120.661211

程序接下来将计算最近的匹配节点，以及这两点之间的最短路径：

start node = (-122.272515, 37.797457)
end node = (-120.66285, 35.285892)
(-122.272515, 37.797457) -> (-122.176834, 37.719054), length = 12528
(-122.176834, 37.719054) -> (-122.176734, 37.718964), length = 13
...
(-120.663604, 35.286751) -> (-120.663466, 35.286594), length = 21
(-120.663466, 35.286594) -> (-120.66285, 35.285892), length = 95
Total length = 358838

当然，像这样打印出每个终点的纬度和经度并不特别有用——如果我们在地图上显示最短路径，对用户来说会更好。如果您愿意，可以将计算出的路径保存到 shapefile 中，然后使用 Mapnik 将 shapefile 的内容作为地图的一部分显示出来。但您可以看到最短路径计算是如何工作的，以及将道路数据转换为 NetworkX 可以处理格式的所需条件。

摘要

在本章中，您学习了两个用于分析地理空间数据的新有用库。然后我们探讨了各种操作和分析空间数据的技术，包括准确计算距离、长度、位置和面积的方法。

接下来，我们探讨了如何将交叉道路转换为平面图，并将其存储为 shapefile，以便基于道路数据进行最短路径计算。最后，我们编写了一个程序来计算任意两点之间的最短路径。在我们解决这些各种问题时，我们学习了许多操作和分析地理空间数据的技术——这些技术您在编写自己的地理空间分析程序时将经常使用。

在下一章中，我们将把您所学的一切结合起来，使用 Python 实现一个完整的地理空间分析系统。

第六章.构建一个完整的地理空间分析系统

在本章中，我们将把在前几章中学到的技能应用到构建一套解决复杂地理空间问题的程序中。在这个过程中，我们将学习：

地图匹配是什么，以及它是如何工作的
如何使用地图匹配根据 GPS 记录生成道路热图
如何下载公路地图数据并将其转换为道路网络
如何将公路网络存储在数据库中
如何使用 GPS 跟踪设备生成自己的旅程记录
如何实现地图匹配算法，将 GPS 记录与现有的公路网络匹配，并使用结果计算每段公路被使用的频率
如何使用这些计算数据生成一个看起来很棒的 GPS 热图

让我们从考察地图匹配的概念开始，看看它如何有助于解决各种地理空间问题。

将 GPS 数据与地图匹配

GPS 记录设备在一段时间内捕捉一系列纬度和经度坐标。随着设备被某人携带从一个地方移动到另一个地方，GPS 坐标记录了人的移动。以下地图显示了典型的 GPS 记录：

GPS 设备相对便宜且非常准确，允许你记录步行、骑自行车、开车或开卡车所经历的旅程。然而，仅凭 GPS 设备本身，它所做的只是记录一系列坐标——GPS 设备不知道你在旅程中跟随了哪些道路。

地图匹配是将 GPS 记录与公路数据库进行匹配的过程，以确定在该旅程中使用了哪组公路。为了使地图匹配成功，你需要三样东西：

一个准确的旅程 GPS 记录，包括足够的 GPS 坐标以确定所走的道路。
一个准确的公路数据库。
一个合适的算法来匹配 GPS 坐标与公路数据库。

一旦你知道了哪些道路被走了，你可以用这些信息做各种用途。例如，一个逐段导航系统将使用其对公路数据库和已走的路径的了解来建议下一步要采取的转弯。

你也可以出于历史目的使用地图匹配：跟踪旅行者在旅程中采取的路线，可能为了优化未来的旅程，或者记录哪些道路被最频繁地使用。

在本章中，我们将实现一个完整的地图匹配系统，使用适当的公路数据库和一种复杂的算法来匹配 GPS 记录与这些道路。我们将使用这些信息来生成一系列历史 GPS 记录中常用公路的热图。

GPS 热图系统的概述

我们正在实施的系统将被命名为GPS 热图。我们将下载一组道路数据，并将这些数据转换为有向道路段网络。然后，我们将生成或下载来自各种旅程的 GPS 记录，我们将使用这些记录来识别常行驶的道路。最后，我们将根据道路被行驶的频率生成热图，提供 GPS 设备捕获的最常行驶道路的视觉摘要。生成的输出将类似于以下内容：

注意

热图通常使用从蓝色到红色的颜色渐变来绘制道路，蓝色用于较少行驶的道路，红色用于最常行驶的道路。然而，在这本书的印刷版中，热图将以黑白形式出现，因此我们选择了一种单色的蓝色，以便在打印时热图仍然有意义。

为了使我们的 GPS 热图程序尽可能高效地运行，我们将使用 PostGIS 数据库来存储底层道路数据。我们将生成非相交道路段的平面图，使用这些数据来构建连接的道路段网络，并将这个网络存储在数据库中以供快速访问。然后，我们将使用地图匹配算法为每个道路段计算一个计数。这些计数然后成为生成 GPS 热图图像的基础，每个道路段将根据该段计算的计数选择合适的颜色。

由于我们的 GPS 热图程序由许多部分组成，我们将将其实现为一套单独的 Python 程序：

init_db.py: 这个脚本将初始化 PostGIS 数据库，为我们提供一个在加载和处理数据时存储数据的地方。
import_roads.py: 这个脚本将道路数据从 shapefile 导入到 PostGIS 数据库中。
split_roads.py: 这个脚本会将导入的道路转换为一系列非重叠的道路段，通过基于原始道路数据计算平面图来实现。
calc_directed_network.py: 这个脚本将使用道路段生成一个连接的道路段有向网络。这告诉我们各种道路段是如何连接的，以及哪些段从一个给定点出发。
map_matcher.py: 这个脚本执行实际的地图匹配过程，从 GPX 格式文件中读取原始 GPS 数据，使用这些数据来识别使用该 GPS 设备的用户所经过的道路段，并在使用时为每个道路段增加计数。
generate_heatmap.py: 这个脚本将使用计算出的计数数据生成用于显示的热图图像。

在我们可以开始实施这些各种程序之前，我们需要获取底层数据。现在让我们看看如何做这件事。

获取必要的数据

您的第一步是决定使用哪个 GPS 数据集。如果您愿意，您可以使用 GPS 设备捕获自己的 GPS 记录，或者您可以使用本章示例代码中提供的 GPS 记录。一旦您决定了要使用哪组 GPS 数据，您将需要下载该区域的道路数据。

让我们先决定您希望使用哪套 GPS 数据。

获取 GPS 数据

如果您有自己的 GPS 记录设备，您可能想出去并在您在当地地区往返时捕获自己的 GPS 记录。如果您没有自己的 GPS 记录设备，或者您不想捕获自己的 GPS 数据，本章提供的示例代码中包含了一些示例 GPS 记录。

这些示例记录是在作者位于新西兰罗托鲁阿市及其周边地区拍摄的。这些记录使用 Garmin Edge 500 自行车电脑捕获，并以GPX格式导出。如果您想捕获自己的 GPS 记录，请确保记录至少 20 次不同的旅程，以便您有一套良好的数据可以工作。

注意

并非每个 GPS 记录都是可用的。有时，GPS 数据可能过于不准确，或者可能错过您的旅程的部分，导致匹配错误。此外，地图匹配算法的限制意味着任何在道路上掉头或使用相同道路段两次的旅程都无法匹配。因此，您可能会发现一些记录不起作用。

无论您使用示例记录还是创建自己的，将生成的GPX格式文件放入名为gps-data的目录中。

下载道路数据

一旦您有了 GPS 记录，接下来您需要一组匹配的道路数据。如果您的 GPS 记录是在美国拍摄的，您可以使用美国人口普查局的 TIGER 数据。我们在上一章中使用了这个网站来下载加利福尼亚州的所有主要和次要道路。在这种情况下，您将想要下载您所在地区的所有道路。为此：

访问www.census.gov/geo/maps-data/data/tiger-line.html，滚动到标记为2014 TIGER/Line Shapefiles的部分，点击下载选项，然后点击Web 界面。
在下载页面，从选择图层类型下拉菜单中选择Roads，然后点击提交按钮。
从所有道路下拉菜单中选择您的州，然后点击提交。
然后，您需要选择您所在的县，最后您可以点击下载按钮以获取所需的道路数据。

如果您不在美国，您需要找到一个合适的替代方案。OpenStreetMap (openstreetmap.org) 是数据的一个可能来源，尽管您可能需要寻找您可以使用格式的道路数据。另外，koordinates.com 网站可能也有您可以使用的数据。

如果你想要使用本章提供的示例 GPS 录音，请按照以下步骤操作：

前往 koordinates.com。
点击页面右上角的登录链接。
通过点击注册链接注册一个免费账户。
登录后，你将看到一个以你的位置为中心的地图。将地图平移到新西兰，并放大到找到北岛中心的城市罗托瓦。
点击页面左上角的搜索框，输入 road centrelines 并按 Enter。

注意

注意新西兰拼写单词 "centrelines"；如果你输入 centerlines，你将找不到你要找的数据。
你要查找的数据集名为 Improved NZ Road Centrelines (August 2011)。点击此数据集旁边的 + 图标，道路数据将显示在你的地图上。
接下来，进一步放大以显示罗托瓦市及其周边地区，点击页面右上角的裁剪工具 ()，并拖出一个大约如下截图所示的矩形：
选择完成后，点击窗口右上角的 下载或订购 链接。默认选项（WGS 84 地图投影和 shapefile 格式）正是你需要的，所以只需点击 接受条款并创建下载 按钮。
大约一分钟后，道路数据将可供你下载。生成的文件将被命名为 kx-improved-nz-road-centrelines-august-2011-SHP.zip。解压此 ZIP 存档，将生成的目录重命名为 roads，并将此目录放置在方便的位置。

实现 GPS 热图系统

现在我们已经获得了必要的数据，我们准备开始实现我们的 GPS 热图系统。创建一个名为 gps-heatmap 的目录来存放程序及其相关数据文件，然后将你之前创建的两个数据目录（gps-data 和 roads）放入此目录中。

现在我们已经准备好开始编码了。让我们先实现 init_db.py 程序来初始化我们的 PostGIS 数据库。

初始化数据库

当你完成第三章 空间数据库 的工作时，你应该已经安装了 Postgres。我们将使用 Postgres 创建并初始化一个数据库来存储所有处理过的道路数据。第一步是创建数据库本身，你可以在终端窗口中输入以下内容来完成：

% createdb gps_heatmap

这将创建一个名为 gps_heatmap 的数据库。如果你遇到认证错误，你需要输入密码或使用 -U postgres 命令行选项，以便 createdb 命令可以运行。

现在您已经创建了数据库本身，下一步是将它转换为空间数据库，以便我们可以用它来存储几何数据。为此，请在终端窗口中输入以下命令：

% psql -d gps_heatmap -c "CREATE EXTENSION postgis;"

小贴士

如果需要，不要忘记添加-U postgres命令行选项。

您现在已经为 Python 代码创建了一个空间数据库。我们现在将编写init_db.py脚本，该脚本初始化数据库中的各种表和索引。请创建主gps-heatmap目录中的init_db.py文件，并将以下代码输入到该文件中：

import psycopg2

connection = psycopg2.connect(database="gps_heatmap",
                              user="postgres")
cursor = connection.cursor()

cursor.execute("DROP TABLE IF EXISTS roads")
cursor.execute("CREATE TABLE roads (" +
                   "id SERIAL PRIMARY KEY," +
                   "name VARCHAR," + 
                   "centerline GEOMETRY)")
cursor.execute("CREATE INDEX ON roads USING GIST(centerline)")

 connection.commit()

如您所见，我们正在使用psycopg2库来访问我们的 PostGIS 数据库。我们创建了一个数据库连接和一个相关的cursor对象。然后我们创建了roads数据库表，首先删除它，这样我们就可以在需要更改数据库结构时再次运行此脚本。

小贴士

如果需要使用不同的用户名或密码访问您计算机上的 PostGIS，请勿忘记更改psycopg2.connect()语句的参数。

roads表将存储从我们刚刚下载的 shapefile 中导入的原始道路数据。如您所见，这个表将有三组不同的字段：

id：这是数据库中这条道路的唯一 ID
name：这是道路的名称
centerline：这是一个 LineString 几何对象，表示这条道路的形状

这不是我们需要的唯一数据库表，但足以让我们开始。随着我们的进行，我们将在init_db.py程序中添加更多的表定义。不过，现在您应该能够运行此程序来创建roads表，这是我们将在下一节中导入下载的道路数据到数据库时需要的。

导入道路数据

我们现在已准备好将下载的 shapefile 中的道路数据导入数据库。将执行此操作的程序称为import_roads.py。请创建此文件，并将以下 Python 代码输入到其中：

import psycopg2
from osgeo import ogr

connection = psycopg2.connect(database="gps_heatmap",
                              user="postgres")
cursor = connection.cursor()

cursor.execute("DELETE FROM roads")

到目前为止，我们所做的一切只是打开数据库连接并删除roads表中的现有内容。接下来，我们需要从我们下载的 shapefile 中导入道路数据。当然，我们如何做将取决于道路数据来自哪里。对于Improved NZ Road Centrelines数据，我们将使用以下代码：

shapefile = ogr.Open("roads/improved-nz-road-centrelines-august-2011.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = feature.GetGeometryRef()

    if feature.GetField("descr") != None:
        name = feature.GetField("descr")
    elif feature.GetField("label") != None:
        name = feature.GetField("label")
    else:
        name = None

    centerline_wkt = geometry.ExportToWkt()

    cursor.execute("INSERT INTO roads (name, centerline) " +
                   "VALUES (%s, ST_GeomFromText(%s))",
                   (name, centerline_wkt))

connection.commit()

如您所见，我们使用 OGR 从 shapefile 导入道路数据，并使用descr或label字段作为道路名称。这对应于NZ Road Centrelines数据定义的方式，有时道路名称在descr字段中，有时在label字段中。有些道路没有名称（例如，道路穿过环岛的地方），因此在这种情况下，道路名称将被设置为None。

如果您使用的是不同的道路数据源，您需要修改此代码以适应您的道路数据组织方式。只需确保将中心线和（如有适用）道路名称存储到roads表中。

分割道路数据成段

如我们在上一章所见，道路接触或交叉的点并不是自动被认为是构建道路网络连接点的。我们首先需要从交叉的道路中创建一个平面图，就像我们在上一章所做的那样。因此，我们的下一个任务是将道路分割成段，从原始道路数据中形成一个平面图的道路段。

split_roads.py程序将分割原始道路数据成段。然而，在我们能够编写这个程序之前，我们需要在数据库中添加一个表来存储道路段。为此，将以下代码添加到您的init_db.py程序中，在调用connection.commit()之前立即执行：

cursor.execute("DROP TABLE IF EXISTS road_segments")
cursor.execute("CREATE TABLE road_segments (" +
               "id SERIAL PRIMARY KEY," +
               "name VARCHAR," + 
               "centerline GEOMETRY," +
               "tally INTEGER)")
cursor.execute("CREATE INDEX ON road_segments USING GIST(centerline)")

如您所见，我们创建了一个名为road_segments的新表来存储平面图中的各个段。每个道路段将包含以下字段：

id：这是该道路段的唯一 ID。
name：这是该段所属的道路名称。
centerline：这是一个表示该道路段形状的 LineString 几何形状。
tally：这是 GPS 记录使用该道路段次数的数量。这是我们在本章后面将要实现的地图匹配算法的输出。

现在我们已经创建了road_segments表，我们可以开始实现split_roads.py程序。创建此文件并将以下代码添加到其中：

import psycopg2
import shapely.wkt
import shapely

connection = psycopg2.connect(database="gps_heatmap",
                              user="postgres")
cursor = connection.cursor()

cursor.execute("DELETE FROM road_segments")

到目前为止，我们只是打开了一个数据库连接并删除了任何现有的road_segments记录。像往常一样，这给了我们一个空白表，我们可以存储我们的计算出的道路段，移除我们可能之前存储的任何段。这允许我们按需运行程序。

接下来，我们想要将roads表的内容转换成一系列连接的道路段。在上一章中，我们使用 Shapely 在内存中存储的道路数据执行了这个任务。这次，我们将使用 PostGIS 实现相同的过程。首先，我们将加载所有road记录 ID 的主列表到内存中：

all_road_ids = []
cursor.execute("SELECT id FROM roads")
for row in cursor:
    all_road_ids.append(row[0])

我们将依次处理每条道路。对于每条道路，我们首先将道路的名称和几何形状加载到内存中：

for road_id in all_road_ids:
    cursor.execute("SELECT name,ST_AsText(centerline) " +
                   "FROM roads WHERE id=%s", (road_id,))
    name,wkt = cursor.fetchone()
    cur_road = shapely.wkt.loads(wkt)

现在我们已经得到了道路的 LineString 几何形状，我们想要在每个它接触或交叉其他道路的点处将其分割。为此，我们将为这条道路构建一个十字路口列表：

    crossroads = []
    cursor.execute("SELECT ST_AsText(centerline) FROM ROADS" +
                   "WHERE ST_Touches(roads.centerline, " +
                   "ST_GeomFromText(%s)) OR ST_Crosses(" +
                   "roads.centerline, ST_GeomFromText(%s))",
                   (wkt, wkt))
    for row in cursor:
        crossroad = shapely.wkt.loads(row[0])
        crossroads.append(crossroad)

然后，我们使用十字路口将当前道路的几何形状分割成一个或多个段：

    for crossroad in crossroads:
        cur_road = cur_road.difference(crossroad)

接下来，我们需要处理生成的道路，将其分割成每个道路段的一个单独的 LineString：

    segments = []
    if cur_road.geom_type == "MultiLineString":
        for segment in cur_road.geoms:
            segments.append(segment)
    elif cur_road.geom_type == "LineString":
        segments.append(cur_road)

然后，我们将计算出的段保存到road_segments表中：

    for segment in segments:
        centerline_wkt = shapely.wkt.dumps(segment)
        cursor.execute("INSERT INTO road_segments (name, " +
                       "centerline, tally) VALUES (%s, " +
                       "ST_GeomFromText(%s), %s)",
                       (name, centerline_wkt, 0))

最后（在for road_id in all_road_ids循环外部），我们将更改提交到数据库：

connection.commit()

这完成了我们的split_roads.py程序。如果你按顺序运行程序，然后使用psql命令行客户端访问数据库，你可以看到程序确实从原始道路数据中生成了许多路段：

% python init_db.py
% python import_roads.py
% python split_roads.py
% psql gps_heatmap
# SELECT count(*) from roads;
1556
# SELECT count(*) from road_segments;
3240

如你所预期，路段的数量比道路的数量要多得多。

注意

如果你想要查看路段，你可以轻松地编写一个使用 Mapnik 的程序来显示road_segments表的内容。这个程序的版本名为preview_segments.py，包含在本章的示例代码中。

构建有向路段的网络

当我们尝试将记录的 GPS 数据与路段数据库匹配时，我们需要回答的一个重要问题是，“从这条路段出发还有哪些其他路段？”为了回答这个问题，我们需要构建一个有向网络的路段。

在上一章中，我们做了类似的事情，在那里我们使用 NetworkX 来计算两点之间的最短路径。然而，在这种情况下，我们将把网络存储在数据库中以供以后使用。

为了使地图匹配算法更容易实现，我们的路段网络将是有向的——也就是说，我们road_segments表中的每个路段实际上将由两个独立的路段表示：

注意

当然，并非每条道路都是双向的，但我们忽略这一点以简化问题。

如你所见，每个路段都被转换成两个有向路段：一个从点 A 到点 B 运行，另一个从点 B 返回到点 A。

我们将使用一个名为directed_segments的新表来存储有向路段。每个有向路段记录将包含以下字段：

id: 这是此路段段落的唯一 ID
road_segment_id: 这是此有向路段派生出的路段记录的 ID
centerline: 这是此有向路段的 LineString 几何形状

注意，有向路段的 LineString 与路段本身的方向一致——也就是说，有向路段的起点在centerline.coords[0]，终点在centerline.coords[-1]。

第二个表，endpoints，将包含各种有向路段的端点坐标。每个端点记录将包含以下字段：

id: 这是此端点的唯一 ID
endpoint: 这是一个包含此端点坐标的 Point 几何形状

注意

我们在这里使用 Point 几何形状，以便可以对这张表进行空间查询。

最后，我们需要一个表来标识从给定端点出发的有向路段。这个表，endpoint_segments，将包含以下字段：

id: 这是此endpoint_segments记录的唯一 ID
directed_segment_id：这是有向道路段的记录 ID
endpoint_id：这是此有向道路段离开的端点的记录 ID

这三个表将用于存储有向道路段网络。让我们修改我们的init_db.py程序以创建这三个新表。为此，将以下代码添加到文件末尾，紧接在调用connection.commit()之前：

cursor.execute("DROP TABLE IF EXISTS directed_segments")
cursor.execute("CREATE TABLE directed_segments (" +
               "id SERIAL PRIMARY KEY," +
               "road_segment_id INTEGER," +
               "centerline GEOMETRY)")
cursor.execute("CREATE INDEX ON directed_segments USING GIST(centerline)")

cursor.execute("DROP TABLE IF EXISTS endpoints")
cursor.execute("CREATE TABLE endpoints (" +
               "id SERIAL PRIMARY KEY," +
               "endpoint GEOMETRY)")
cursor.execute("CREATE INDEX ON endpoints USING GIST(endpoint)")

cursor.execute("DROP TABLE IF EXISTS endpoint_segments")
cursor.execute("CREATE TABLE endpoint_segments (" +
               "id SERIAL PRIMARY KEY," +
               "directed_segment_id INTEGER," +
               "endpoint_id INTEGER)")
cursor.execute("CREATE INDEX ON endpoint_segments(directed_segment_id)")
cursor.execute("CREATE INDEX ON endpoint_segments(endpoint_id)")

这是我们需要对数据库结构进行的最后一次更改，所以请继续重新创建数据库表，导入道路，并再次分割它们：

% python init_db.py
% python import_roads.py
% python split_roads.py

注意

如果我们使用了数据库迁移，我们就不需要每次都重新运行我们的程序，但我们正在尽可能简化数据库逻辑。幸运的是，这是我们最后一次需要这样做的时候。

我们现在已准备好计算有向道路网络并将其存储到数据库中。我们将创建的程序称为calc_directed_network.py；创建此文件，并将以下代码输入其中：

import networkx
import psycopg2
import shapely.wkt
import shapely.geometry

connection = psycopg2.connect(database="gps_heatmap", user="postgres")
cursor = connection.cursor()

我们现在已准备好创建表示道路网络的 NetworkX 图。在上一章中我们这样做时，我们使用了networkx.read_shp()函数，从 shapefile 的内容创建一个 NetworkX DiGraph对象。不幸的是，没有等效函数可以用于从 PostGIS 数据库的内容创建图；然而，由于 NetworkX 是用 Python 实现的，因此很容易修改read_shp()函数的源代码以实现我们想要的功能。将以下代码添加到calc_directed_network.py程序末尾：

network = networkx.Graph()

cursor.execute("SELECT id,ST_AsText(centerline) FROM road_segments")
for row in cursor:
    road_segment_id,wkt = row
    linestring = shapely.wkt.loads(wkt)

    first_pt = linestring.coords[0]
    last_pt  = linestring.coords[-1]

    network.add_edge(first_pt, last_pt,
                     {'road_segment_id' : road_segment_id})

NetworkX 图中的节点是一个(long,lat)元组，用于标识每个道路段端点的唯一标识符，边代表有向道路段。请注意，我们在图中存储道路段的记录 ID 作为属性，以便以后引用。

现在我们有了 NetworkX 图，让我们准备使用它来生成连接道路段的有向网络。为此，我们必须从我们构建的图中计算最大的连接子图，就像我们在上一章中所做的那样。以下是必要的代码：

sub_graphs = list(networkx.connected_component_subgraphs(network))
largest = sub_graphs[0]

我们现在有一个包含所有连接道路段的图。我们可以现在使用它将道路段端点存储到数据库中：

cursor.execute("DELETE FROM endpoints")

endpoint_ids = {}
for node in largest.nodes():
    point = shapely.geometry.Point(node)
    wkt = shapely.wkt.dumps(point)

    cursor.execute("INSERT INTO endpoints (endpoint) " +
                   "VALUES (ST_GeomFromText(%s)) RETURNING id",
                   (wkt,))
    endpoint_id = cursor.fetchone()[0]

    endpoint_ids[node] = endpoint_id

注意，endpoint_ids字典将一个(long,lat)坐标映射到数据库中端点的记录 ID。我们将使用它将有向道路段与其端点链接起来。

我们的最后任务是存储有向道路段，以及该段开始于的端点。我们将首先删除数据库中现有的记录，然后遍历我们图中的道路段：

cursor.execute("DELETE FROM directed_segments")
cursor.execute("DELETE FROM endpoint_segments")

for node1,node2 in largest.edges():
    endpoint_id_1 = endpoint_ids[node1]
    endpoint_id_2 = endpoint_ids[node2]
    road_segment_id = largest.get_edge_data(node1, node2)['road_segment_id']

    cursor.execute("SELECT ST_AsText(centerline) " +
                   "FROM road_segments WHERE id=%s",
                   (road_segment_id,))
    wkt = cursor.fetchone()[0]
    linestring = shapely.wkt.loads(wkt)

我们现在有了段端点的记录 ID 和定义此道路段的 LineString。我们现在需要将此段转换为两个有向段，每个方向一个：

    reversed_coords = list(reversed(linestring.coords))
    if node1 == linestring.coords[0]:
        forward_linestring = linestring
        reverse_linestring = shapely.geometry.LineString(reversed_coords)
    else:
        reverse_linestring = linestring
        forward_linestring = shapely.geometry.LineString(reversed_coords)

这给我们提供了两个 LineString 几何形状，一个从第一个端点到第二个端点，另一个从第二个端点回到第一个端点。现在我们可以将我们计算的信息保存到数据库中：

    cursor.execute("INSERT INTO directed_segments " +
                   "(road_segment_id,centerline) VALUES " +
                   "(%s, ST_GeomFromText(%s)) RETURNING id",
                   (road_segment_id, forward_linestring.wkt))
    forward_segment_id = cursor.fetchone()[0]

    cursor.execute("INSERT INTO directed_segments " +
                   "(road_segment_id,centerline) VALUES " +
                   "(%s, ST_GeomFromText(%s)) RETURNING id",
                   (road_segment_id, reverse_linestring.wkt))
    reverse_segment_id = cursor.fetchone()[0]

    cursor.execute("INSERT INTO endpoint_segments " + 
                       "(directed_segment_id, endpoint_id) " +

                       "VALUES (%s, %s)", 

                       (forward_segment_id, endpoint_id_1))

        cursor.execute("INSERT INTO endpoint_segments " + 
                       "(directed_segment_id, endpoint_id) " +

                       "VALUES (%s, %s)", 

                       (reverse_segment_id, endpoint_id_2))

为了完成我们的程序，我们必须提交我们所做的更改：

connection.commit()

你现在可以通过运行这个程序来创建有向道路段网络：

% python calc_directed_network.py

实现地图匹配算法

我们终于准备好实现我们的地图匹配算法了。我们将使用的算法来源于论文：《使用多假设技术（MHT）在高清导航网络上进行 GPS 轨迹的地图匹配》，由 Nadine Schuessler 和 Kay Axhausen 为瑞士交通规划与系统研究所撰写。

如果你对阅读这篇论文感兴趣，可以在www.ivt.ethz.ch/vpl/publications/reports/ab568.pdf找到。这个算法基于路线候选的概念，即旅行者在记录 GPS 点时可能采取的路径。每个路线候选都有一个包含方向道路段和得分的列表，该得分标识 GPS 点与这些道路段的匹配程度。

通过逐个跟随 GPS 点来重现旅程。在任何特定时刻，都有一个可能匹配迄今为止已处理的 GPS 坐标的路线候选列表。随着每个新的 GPS 点被处理，路线候选会逐个更新，通过比较 GPS 坐标与路线候选的最终道路段。

如果 GPS 点被认为仍然在最终段落的某个地方，那么 GPS 点就简单地添加到该段落，并更新路线候选的得分。另一方面，如果 GPS 点超出了路线候选的最终道路段，那么我们会查看道路网络，看看哪些其他道路段从这个点出发。然后，我们为每个从该端点出发的道路段创建新的路线候选：

一旦处理完所有 GPS 点，我们就选择得分最低的路线候选，认为它最有可能是这次旅行的使用路径。

为了实现这个算法，我们将使用 Python 字典来表示路线段——即在路线候选的旅程中的单个段落。每个路线段字典将包含以下条目：

directed_segment_id：这是跟随该路线段的道路段directed_segment的记录 ID
linestring：这是道路段中心线，作为一个 Shapely LineString 对象
gps_points：这是一个定义分配给该路线段的 GPS 点的(long,lat)坐标的列表
gps_distances：这是一个包含每个 GPS 点与段 LineString 之间计算出的最小距离的列表

每个路由候选将用一个包含以下条目的 Python 字典来表示：

segments: 这是一个由构成此路由候选的路由段组成的列表。
directed_segment_ids: 这是一个包含此路由所使用的每个定向段记录 ID 的列表。我们使用这个列表来快速丢弃另一个路由正在使用相同道路段序列的新路由候选。
score: 这是为此路由候选计算出的分数。分数是每个路由段中 GPS 距离的总和——换句话说，分数越低，GPS 点越接近此路由。

在了解这些信息的基础上，让我们开始实现地图匹配器。创建一个新的 Python 程序，命名为map_matcher.py，并将以下内容输入到该文件中：

import os
import osgeo.ogr
import shapely.geometry
import shapely.wkt
import psycopg2
import pyproj

gps_tracks = []
for fName in os.listdir("gps-data"):
    if fName.endswith(".gpx"):
        srcFile = osgeo.ogr.Open("gps-data/" + fName)
        layer = srcFile.GetLayerByName("tracks")

        for feature_num in range(layer.GetFeatureCount()):
            feature = layer.GetFeature(feature_num)
            geometry = feature.GetGeometryRef()

            if geometry.GetGeometryName() == "MULTILINESTRING":
                for geom_num in range(geometry.GetGeometryCount()):
                    wkt = geometry.GetGeometryRef(geom_num).ExportToWkt()
                    gps_tracks.append((fName, wkt))
            elif geometry.GetGeometryName() == "LINESTRING":
                wkt = geometry.ExportToWkt()
                gps_tracks.append((fName, wkt))

connection = psycopg2.connect(database="gps_heatmap", user="postgres")
cursor = connection.cursor()

如您所见，我们导入了所需的各个库，将记录的 GPS 数据加载到内存中，并打开到我们的数据库的连接。接下来，我们想要重置道路段的总计值：

cursor.execute("UPDATE road_segments SET tally=0")
connection.commit()

我们现在准备好开始处理记录的 GPS 数据。虽然基于此算法的论文使用了一种复杂的过程，将 GPS 数据分割成行程段，处理每个行程段，然后将结果路由连接起来，但我们将使用一种更简单的方法；我们假设每个 GPS 记录没有间隙，但它可能开始或结束在道路上。为了处理这个问题，我们修剪起点和终点坐标，直到我们找到一个距离道路 10 米以内的点。

将以下内容添加到程序末尾：

for fName,track_wkt in gps_tracks:
    print "Processing " + fName

    gps_track  = shapely.wkt.loads(track_wkt)
    gps_points = list(gps_track.coords)

    while len(gps_points) > 0:
        circle = calc_circle_with_radius(gps_points[0], 10)
        cursor.execute("SELECT count(*) FROM road_segments " +
                "WHERE ST_Intersects(ST_GeomFromText(%s)," +
                       "centerline)", (circle.wkt,))
        if cursor.fetchone()[0] == 0:
            del gps_points[0]
        else:
            break

    while len(gps_points) > 0:
        circle = calc_circle_with_radius(gps_points[-1], 10)
        cursor.execute("SELECT count(*) FROM road_segments " +
                "WHERE ST_Intersects(ST_GeomFromText(%s)," +
                       "centerline)", (circle.wkt,))
        if cursor.fetchone()[0] == 0:
            del gps_points[-1]
        else:
            break

我们只是依次处理每个 GPS 轨迹，从起点和终点修剪点，直到找到一个距离道路段 10 米以内的 GPS 点。请注意，我们正在使用一个名为calc_circle_with_radius()的函数来创建一个描述 GPS 坐标周围 10 米圆的 Shapely 多边形，然后要求数据库在该圆内查找任何道路段。让我们继续实现这个calc_circle_with_radius()函数；将其放置在程序顶部，紧接在import语句之后：

def calc_circle_with_radius(center_point, radius):
    geod = pyproj.Geod(ellps="WGS84")
    sLong,sLat = center_point
    eLong,eLat,iHeading = geod.fwd(sLong, sLat, 0, radius)
    lat_delta = abs(sLat - eLat)
    return shapely.geometry.Point(sLong, sLat).buffer(lat_delta)

现在我们已经得到了每个记录中相关的 GPS 点集合，我们准备好开始地图匹配。第一步是根据起始 GPS 点构建一组初始路由候选。我们通过识别距离 GPS 点 750 米以内的所有道路端点来完成此操作，并为从这些端点出发的每个道路段创建一个（单段）路由候选。遵循 Schuessler 和 Axhausen 的论文，我们确保至少有 25 个路由候选，如果没有，我们将搜索区域扩大 100 米并再次尝试。

将以下代码添加到你的程序末尾：

    search_distance = 750
    while True:
        circle = calc_circle_with_radius(gps_points[0],
                                         search_distance)

        cursor.execute("SELECT id FROM endpoints " +
                       "WHERE ST_Contains(ST_GeomFromText(%s)," +
                       "endpoint)", (circle.wkt,))
        possible_endpoints = []
        for row in cursor:
            possible_endpoints.append(row[0])

        possible_road_segments = []
        for endpoint_id in possible_endpoints:
            cursor.execute("SELECT directed_segment_id " +
                           "FROM endpoint_segments " +
                           "WHERE endpoint_id=%s", (endpoint_id,))
            for row in cursor:
                directed_segment_id = row[0]
                possible_road_segments.append(
                    (directed_segment_id, endpoint_id))

        route_candidates = []
        for directed_segment_id,endpoint_id in possible_road_segments:
            cursor.execute("SELECT ST_AsText(centerline) " +
                           "FROM directed_segments WHERE id=%s",
                           (directed_segment_id,))
            wkt = cursor.fetchone()[0]
            linestring = shapely.wkt.loads(wkt)
            gps_distance = calc_distance(gps_points[0],
                                         linestring)

            segment = {
                'directed_segment_id' : directed_segment_id,
                'linestring' : linestring,
                'gps_points': [gps_points[0]],
                'gps_distances': [gps_distance]}
            route_segments = [segment]

            candidate = {
                'segments': route_segments,
                'directed_segment_ids' : [directed_segment_id],
                'score': calc_score(route_segments)}
            route_candidates.append(candidate)

        if len(route_candidates) >= 25:
            break
        else:
            search_distance = search_distance + 100
            continue

如您所见，我们为每个可能的路线候选者创建一个单独的路线段字典和一个路线候选者字典，并将结果存储在route_candidates列表中。这里我们还需要两个额外的函数：一个用于计算从给定点到给定 Shapely 几何体内部最近点的距离，另一个用于计算给定路线候选者的得分。请将这些两个函数添加到程序顶部：

def calc_distance(point, geometry):
    return shapely.geometry.Point(point).distance(geometry)

def calc_score(route_segments):
    total = 0
    for segment in route_segments:
        total = total + sum(segment['gps_distances'])
    return total

现在我们已经有一组初始路线候选者，我们必须通过将每个连续的 GPS 点添加到每个路线候选者中来开发它们，在到达每个路段的末尾时生成新的路线候选者。同时，我们修剪路线候选者列表，以防止其变得过大。

我们的大部分工作将在名为develop_route()的函数中完成。这个函数将接受一个路线候选者和一个 GPS 点（以及一些其他参数），并返回一个包含新或更新路线候选者的列表，以便在下一轮迭代中进行处理。让我们编写使用此函数的代码；将以下内容添加到程序末尾：

    for next_point in gps_points[1:]:
        num_routes_to_process = len(route_candidates)
        for i in range(num_routes_to_process):
            route = route_candidates.pop(0)
            new_candidates = develop_route(next_point, route, route_candidates, cursor)
            route_candidates.extend(new_candidates)

我们对每个路线候选者只处理一次，首先使用route_candidates.pop(0)将其从候选者列表中删除，然后将候选者传递给develop_route()函数。然后我们将新或更新的路线候选者添加到route_candidates列表的末尾。当我们的for i in range(num_routes_to_process)循环结束时，我们将已经处理了每个路线候选者一次，要么将 GPS 点纳入该路线候选者，要么用一组新的路线候选者替换它。

在我们开始处理下一个 GPS 点之前，我们需要修剪路线候选者列表。根据 Schuessler 和 Axhausen 的说法，一种非常有效的方法是反复移除得分最高的路线候选者，直到剩余的候选者不超过 40 个。我们现在就来做这件事：

        while len(route_candidates) > 40:
            highest = None
            for index,route in enumerate(route_candidates):
                if highest == None:
                    highest = index
                elif route['score'] > route_candidates[highest]['score']:
                    highest = index
            del route_candidates[highest]

小贴士

确保将此代码放在for next_point in...循环内。

在我们实现develop_route()函数之前，让我们完成程序的主要部分。我们现在已经处理了所有 GPS 点，因此我们可以检查每个剩余路线候选者的得分，并选择得分最低的候选者（排除任何路线段少于两个的候选者）。这是 GPS 点最有可能采取的路线候选者。然后我们增加该路线使用的每个路段的计数。以下是相关代码：

    best_route = None
    for route in route_candidates:
        if len(route['segments']) >= 2:
            if best_route == None:
                best_route = route
            elif route['score'] < best_route['score']:
                best_route = route

    if best_route == None: continue

    for segment in best_route['segments']:
        cursor.execute("SELECT road_segment_id " +
                       "FROM directed_segments WHERE id=%s",
                       (segment['directed_segment_id'],))
        road_segment_id = cursor.fetchone()[0]
        cursor.execute("UPDATE road_segments SET tally=tally+1" +
                       "WHERE id=%s", (road_segment_id,))

connection.commit()

我们现在需要实现develop_route()函数。这个函数将使用 Schuessler 和 Axhausen 论文中的以下逻辑：

如果路线候选者只有一个路段，检查 GPS 点是否到达该路段 LineString 的起点。如果发生这种情况，GPS 记录必须是在错误方向上跟随该有向路段，因此我们丢弃该路线候选者。
检查 GPS 点是否仍然位于路线候选者的最终段内。如果是这样，将该 GPS 点添加到该最终段，重新计算候选者的得分，并返回它以供进一步处理。
如果 GPS 点超出了路线候选者最终段落的末尾，识别我们已达到的端点，并为从该端点出发的每个有向道路段创建一个新的路线候选者。我们检查每个新路线候选者的有效性，并返回有效候选者以供进一步处理。

让我们开始实现这个函数：

def develop_route(next_point, route, route_candidates, cursor):
    if len(route['segments']) == 1:
        if point_at_start_of_segment(next_point,
                                     route['segments'][0]):
            return []

    last_segment = route['segments'][-1]

    if point_in_route_segment(next_point, last_segment):
        next_distance = calc_distance(next_point,
                                      last_segment['linestring'])
        last_segment['gps_points'].append(next_point)
        last_segment['gps_distances'].append(next_distance)
        route['score'] = calc_score(route['segments'])
        return [route]

这实现了之前描述的前两个步骤。请注意，我们使用了两个新的函数，point_at_start_of_segment() 和 point_in_route_segment() 来完成所有繁重的工作。我们将在不久后实现这些函数，但首先让我们通过创建一组新的路线候选集来处理 GPS 点经过最后一个路线段末尾的过程。

这个过程的第一个步骤是确定我们已达到的当前端点。将以下内容添加到 develop_route() 函数的末尾：

    last_point = last_segment['linestring'].coords[-1]
    endpoint = shapely.geometry.Point(last_point)

    cursor.execute("SELECT id FROM endpoints " +
                   "WHERE endpoint=ST_GeomFromText(%s)",
                   (endpoint.wkt,))
    endpoint_id = cursor.fetchone()[0]

接下来，我们将构建一个包含所有从该端点出发的有向道路段列表：

    possible_segment_ids = []
    cursor.execute("SELECT directed_segment_id " +
                   "FROM endpoint_segments " +
                   "WHERE endpoint_id=%s", (endpoint_id,))
    for row in cursor:
        possible_segment_ids.append(row[0])

我们现在需要为每个可能的路段创建一个新的路线候选者。对于每个路线候选者，我们使用该有向路段创建一个单独的路线段：

    new_candidates = []
    for directed_segment_id in possible_segment_ids:
        cursor.execute("SELECT road_segment_id," +
                       "ST_AsText(centerline) " +
                       "FROM directed_segments " +
                       "WHERE id=%s", (directed_segment_id,))
        road_segment_id,wkt = cursor.fetchone()
        linestring = shapely.wkt.loads(wkt)

        next_distance = calc_distance(next_point, linestring)

        new_segment = {}
        new_segment['directed_segment_id'] = directed_segment_id
        new_segment['linestring'] = linestring
        new_segment['gps_points'] = [next_point]
        new_segment['gps_distances'] = [next_distance]

        new_candidate = {}
        new_candidate['segments'] = []
        new_candidate['segments'].extend(route['segments'])
        new_candidate['segments'].append(new_segment)
        new_candidate['directed_segment_ids'] = []
        new_candidate['directed_segment_ids'].extend(
                            route['directed_segment_ids'])
        new_candidate['directed_segment_ids'].append(directed_segment_id)

        if not route_is_valid(new_candidate, route_candidates,
                              new_candidates):
            continue

        new_candidate['score'] = calc_score(new_candidate['segments'])
        new_candidates.append(new_candidate)
    return new_candidates

注意，我们使用另一个函数 route_is_valid() 来检查新路线的有效性。我们还将必须实现这个函数。

这完成了 develop_route() 函数本身。现在让我们编写 point_at_start_of_segment() 函数，该函数确定 GPS 轨迹是否沿着有向道路段错误地运行：

def point_at_start_of_segment(next_point, segment):
    num_points = len(segment['gps_points'])
    if num_points > 0:
        average_distance = sum(segment['gps_distances']) / num_points

        startpoint_coord = segment['linestring'].coords[0]
        startpoint = shapely.geometry.Point(startpoint_coord)
        endpoint_coord = segment['linestring'].coords[-1]
        endpoint = shapely.geometry.Point(endpoint_coord)

        distance_to_start = calc_distance(next_point, startpoint)
        distance_to_end   = calc_distance(next_point, endpoint)

        if distance_to_start < 2 * average_distance:
            if distance_to_end > 2 * average_distance:
                return True
    return False

这段代码有点混乱，比较当前点到路段起点和终点的距离，但对我们来说已经足够好了。

接下来，我们需要实现 point_in_route_segment() 函数。我们将使用两个单独的测试来查看点是否已达到段落的端点。首先，我们知道如果 GPS 点到段落的 LineString 上最近点的距离等于点到该 LineString 末端的距离，那么我们已经到达了端点：

这是 point_in_route_segment() 函数的第一部分，它实现了这个测试：

def point_in_route_segment(point, segment):
    endpoint = shapely.geometry.Point(segment['linestring'].coords[-1])

    distance_to_linestring = calc_distance(point,
                                           segment['linestring'])
    distance_to_endpoint = calc_distance(point, endpoint)

    if distance_to_linestring == distance_to_endpoint:
        return False

第二个测试涉及比较最终路线段与由分配给该路线段的 GPS 点构建的 LineString 的长度。如果 GPS LineString 比路段长，那么我们必须已经到达了该段落的末尾：

    gps_coords = []
    gps_coords.extend(segment['gps_points'])
    gps_coords.append(point)

    gps_length = shapely.geometry.LineString(gps_coords).length
    segment_length = segment['linestring'].length

    if gps_length > segment_length:
        return False

最后，如果 GPS 点未通过这两个测试，那么它仍然位于当前路线段内：

    return True

注意

Schuessler 和 Axhausen 的论文提出了一种第三种测试方法，即比较 GPS 轨迹的方向与路段方向。然而，由于路段是复杂的 LineStrings 而不是直线段，因此不清楚如何实现这种测试，所以我们不会在我们的地图匹配算法实现中使用这种测试。

这完成了 point_in_route_segment() 函数的实现。我们需要实现的最后一个函数是 route_is_valid()。如果一条路线候选被认为是有效的，那么：

它是唯一的；也就是说，没有其他路线候选具有完全相同的路段序列
其最终路段不会回到前一个路段的起点；也就是说，路线不会自我回环。
路线不包括相同的定向路段两次

为了计算唯一性，route_is_valid() 函数不仅需要一个当前路线候选的列表，还需要一个由 develop_route() 函数创建的新候选列表。因此，route_is_valid() 函数接受当前路线候选列表和正在创建的新候选列表。

这是该函数实现的第一部分，包括唯一性检查：

def route_is_valid(route, route_candidates, new_candidates):
    route_roads = route['directed_segment_ids']

    for other_route in route_candidates:
        if route_roads == other_route['directed_segment_ids']:
            return False

    for other_route in new_candidates:
        if route_roads == other_route['directed_segment_ids']:
            return False

以下代码检查一条路线是否不会自我回环：

    if len(route['segments']) >= 2:
        last_segment = route['segments'][-1]
        prev_segment = route['segments'][-2]

        last_segment_end   = last_segment['linestring'].coords[-1]
        prev_segment_start = prev_segment['linestring'].coords[0]

        if last_segment_end == prev_segment_start:
            return False

最后，我们确保相同的定向路段不会被使用两次：

    directed_segment_ids = set()
    for segment in route['segments']:
        directed_segment_id = segment['directed_segment_id']
        if directed_segment_id in directed_segment_ids:
            return False
        else:
            directed_segment_ids.add(directed_segment_id)

如果路线通过所有三个检查，则被认为是有效的：

    return True

这完成了 route_is_valid() 函数的实现，实际上也完成了整个 map_matcher.py 程序的实现。您应该能够从命令行运行它，并依次看到每个 GPS 记录被处理：

% python map_matcher.py
Processing ride_2015_01_08.gpx
Processing ride_2015_01_11.gpx
Processing ride_2015_01_23.gpx
...

由于每个 GPS 记录中都有数千个点，程序处理每个文件将需要几分钟。一旦完成，road_segments 表中的 tally 字段就会更新，以显示每个路段被使用的次数。您可以使用 Postgres 命令行客户端进行检查：

% psql gps_heatmap
# SELECT name,tally FROM road_segments WHERE tally > 0 ORDER BY tally DESC;
 3560 | otonga rd                        |    42
 6344 | wychwood cres                    |    42
 3561 | otonga rd                        |    42
 3557 | otonga rd                        |    42
 3558 | otonga rd                        |    42
 3559 | otonga rd                        |    42
 6343 | wychwood cres                    |    41
 6246 | springfield rd                   |    19
 6300 | old taupo rd                     |    19

如您所见，地图匹配是一个相当复杂的过程，但这个程序实际上工作得相当好。现在我们已经计算了计数，我们可以编写我们 GPS 热力图系统的最后一部分：基于计算计数值显示热力图的程序。

生成 GPS 热力图

我们将使用 Mapnik 生成热力图，为每个独特的计数值创建一个单独的 mapnik.Rule，这样每个路段使用的颜色就会根据其计数值而变化。这个程序将被命名为 generate_heatmap.py；创建这个程序并将以下代码输入其中：

import mapnik
import psycopg2

MAX_WIDTH = 1200
MAX_HEIGHT = 800
MIN_TALLY = 3

connection = psycopg2.connect(database="gps_heatmap",
                              user="postgres")
cursor = connection.cursor()

在导入所需的库和定义一些常量之后，我们打开数据库连接，以便我们可以计算最高的计数值和计算热力图的边界。现在让我们来做这件事：

cursor.execute("SELECT max(tally) FROM road_segments")
max_tally = cursor.fetchone()[0]

cursor.execute("SELECT ST_XMIN(ST_EXTENT(centerline)), " +
               "ST_YMIN(ST_EXTENT(centerline)), " +
               "ST_XMAX(ST_EXTENT(centerline)), " +
               "ST_YMAX(ST_EXTENT(centerline)) " +
               "FROM road_segments WHERE tally >= %s" % MIN_TALLY)
min_long,min_lat,max_long,max_lat = cursor.fetchone()

如你所见，我们使用MIN_TALLY常量来放大热图中更受欢迎的部分。如果你想改变这个值，将其设置为1将显示所有被 GPS 轨迹覆盖的道路段，而将其设置为更高的值将聚焦于地图上最常用的部分。

现在我们知道了热图覆盖的地球区域，我们可以计算地图图像的尺寸。我们希望使用指定的最大尺寸，同时保持地图的宽高比：

extent = mapnik.Envelope(min_long, min_lat,  max_long, max_lat)
aspectRatio = extent.width() / extent.height()

mapWidth = MAX_WIDTH
mapHeight = int(mapWidth / aspectRatio)
if mapHeight > MAX_HEIGHT:
    scaleFactor = float(MAX_HEIGHT) / float(mapHeight)
    mapWidth = int(mapWidth * scaleFactor)
    mapHeight = int(mapHeight * scaleFactor)

接下来，我们初始化地图本身：

map = mapnik.Map(mapWidth, mapHeight)
map.background = mapnik.Color("white")

尽管只有部分道路段会被 GPS 记录使用，但我们仍然想显示所有未使用的道路段作为热图的背景。为此，我们将创建一个unused_roads图层和相应的 Mapnik 样式：

layer = mapnik.Layer("unused_roads")
layer.datasource = mapnik.PostGIS(host='localhost',
                                  user='postgres',
                                  password='',
                                  dbname='gps_heatmap',
                                  table='road_segments')
layer.styles.append("unused_road_style")
map.layers.append(layer)

line_symbol = mapnik.LineSymbolizer(mapnik.Color("#c0c0c0"), 1.0)

rule = mapnik.Rule()
rule.filter = mapnik.Filter("[tally] = 0")
rule.symbols.append(line_symbol)

style = mapnik.Style()
style.rules.append(rule)
map.append_style("unused_road_style", style)

注意，我们使用mapnik.PostGIS()数据源，这样地图图层就可以直接从我们的 PostGIS 数据库中获取数据。

接下来，我们需要定义一个用于已使用道路（即具有tally值为1或更多的道路）的地图图层。这个我们将称之为used_roads的图层将为每个独特的tally值有一个单独的mapnik.Rule()。这允许我们为每个独特的tally值分配不同的颜色，以便每个道路段的颜色根据该段tally值而变化。为了实现这一点，我们需要一个函数来计算给定tally值要使用的mapnik.Stroke()。以下是这个函数，你应该将其放置在程序顶部附近：

def calc_stroke(value, max_value):
    fraction = float(value) / float(max_value)

    def interpolate(start_value, end_value, fraction):
        return start_value + (end_value - start_value) * fraction

    r = interpolate(0.7, 0.0, fraction)
    g = interpolate(0.7, 0.0, fraction)
    b = interpolate(1.0, 0.4, fraction)

    color = mapnik.Color(int(r*255), int(g*255), int(b*255))
    width = max(4.0 * fraction, 1.5)

    return mapnik.Stroke(color, width)

interpolate()辅助函数用于计算从浅蓝色到深蓝色的颜色范围。我们还根据tally调整显示的道路段宽度，以便更常用的道路在地图上用更宽的线条绘制。

小贴士

如果你想，你可以更改起始和结束颜色，使热图更加多彩。如前所述，我们只是使用了蓝色调，这样当以黑白打印时热图才有意义。

实现了这个功能后，我们可以将used_roads图层添加到我们的地图中。为此，请将以下代码添加到你的程序末尾：

layer = mapnik.Layer("used_roads")
layer.datasource = mapnik.PostGIS(host='localhost',
                                  user='postgres',
                                  password='',
                                  dbname='gps_heatmap',
                                  table='road_segments')
layer.styles.append("used_road_style")
map.layers.append(layer)

style = mapnik.Style()
for tally in range(1, max_tally+1):
    line_symbol = mapnik.LineSymbolizer()
    line_symbol.stroke = calc_stroke(tally, max_tally)

    rule = mapnik.Rule()
    rule.filter = mapnik.Filter("[tally] = %d" % tally)
    rule.symbols.append(line_symbol)

    style.rules.append(rule)
map.append_style("used_road_style", style)

最后，我们可以将地图渲染出来，并将结果保存到磁盘上的图像文件中：

map.zoom_to_box(extent)
mapnik.render_to_file(map, "heatmap.png", "png")

运行此程序后，你应该得到一个包含生成的热图的heatmap.png文件：

恭喜！这个程序绝非平凡，在生成这张图像的过程中解决了许多地理空间问题。当然，你可以使用这个程序来匹配你自己的 GPS 记录与道路网络，但我们真正展示的是如何一步一步地解决复杂的地理空间问题，使用本书中描述的各种技术。

进一步改进

虽然 GPS 热图系统工作得相当不错，但它并不完美。没有任何程序是完美的。如果你有兴趣，你可能想考虑以下方面：

使路段分割算法更复杂，以支持单行道路，以及两条道路相交但不连接的点（例如，在高速公路上跨线桥上）。
改进路线开发过程，使其能够捕捉包括掉头和重复路段的路线。
将原始 GPS 数据分割成连续的段，依次处理每个段，然后将处理过的段重新连接起来。这将允许算法处理包含数据缺失的 GPS 记录。
将地图匹配算法与最短路径计算相结合，编写你自己的逐点导航系统。
研究提高地图匹配算法速度的方法。例如，如果两个不同的路线候选者使用相同的路段，这两个候选者应该能够共享点到该路段之间的计算距离。这将避免需要计算两次相同的距离。肯定还有其他方法可以优化地图匹配器，使其运行得更快。
将一个漂亮的栅格底图图像添加到生成的热图中。

摘要

恭喜！你已经完成了实现一系列程序，这些程序使用一系列地理空间分析技术将记录的 GPS 数据与现有的道路网络匹配。在创建 GPS 热图系统的过程中，你学习了如何将现有的道路数据转换为网络，如何在数据库中表示道路网络，以及如何使用这些数据来实现复杂的地图匹配算法。然后，该算法被用来计算记录的 GPS 数据中每个路段被使用的频率，然后使用这些计数来生成一个显示最常用道路的漂亮热图。

即使你对捕捉自己的 GPS 数据并将其与地图匹配不感兴趣，本章中使用的技巧也会给你自己的开发工作提供许多想法。使用 Python、GDAL 和 OGR、Shapely、PyProj、PostGIS 和 Mapnik 的组合，你现在拥有了一套处理、分析和显示地理空间数据的强大工具。要了解更多信息，请查看以下参考资料：

GDAL 和 OGR 库的主网站
GDAL 和 OGR 的 Python 接口的整体设计描述
使用 Python 的 GDAL/OGR 时需要注意的一些问题列表
Shapely 库的主网站
Shapely 用户手册可以在这里找到
Proj.4 地图投影库的描述
jswhit.github.io/pyproj 是 Proj 库的 Python 接口的主要网站。
www.postgresql.org 是 PostgreSQL 数据库的主要网站。
postgis.net 是 PostGIS 扩展的网站。
mapnik.org 是 Mapnik 库的主要网站。
最后，github.com/mapnik/mapnik/wiki 提供了大量有用的信息，以帮助你使用 Mapnik。

这本书现在已经完成，但我希望我已经让你对可用于地理空间分析的工具和技术有了更深入的理解，以及 Python 如何作为你自己的地理空间开发工作的基础。其余的取决于你。

posted @ 2025-10-23 15:18 绝不原创的飞龙阅读(25) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈