Python-地理空间分析学习手册-三-Python 地理空间分析学习手册（三）七、Python 和高程数据高程数据

Python 地理空间分析学习手册（三）

原文：Learning Geospatial Analysis with Python

协议：CC BY-NC-SA 4.0

七、Python 和高程数据

高程数据是最迷人的地理空间数据类型之一。它代表许多不同类型的数据源和格式。它可以显示向量和栅格数据的属性，从而产生独特的数据产品。高程数据可用于地形可视化、土地覆盖分类、水文建模、运输路线、特征提取和许多其他目的。

栅格数据和向量数据不能执行所有这些选项，但是由于高程数据是三维的，由于包含 x 、 y 和 z 坐标，您通常可以从这些数据中获得比任何其他类型更多的信息。

在本章中，我们将涵盖以下主题:

使用 ASCII 网格高程数据文件进行简单的高程处理
创建阴影浮雕图像
创建高程等高线
激光雷达数据网格化
创建三维网格

在本章中，您将学习如何读写栅格和向量格式的高程数据。我们还将创造一些衍生产品。

访问 ASCII 网格文件

在本章的大部分内容中，我们将使用 ASCII 网格文件。这些文件是一种栅格数据，通常与高程数据相关联。这种网格格式将数据以文本形式存储在大小相等的方形行和列中，并带有简单的标题。行/列中的每个单元格存储一个数值，该数值可以表示地形的某些特征，如海拔、坡度或流向。简单性使其成为一种易于使用且独立于平台的栅格格式。这种格式在第二章学习地理空间数据的 ASCII 网格一节中有描述。

在本书中，我们一直依赖 GDAL，在某种程度上甚至依赖 PIL 来读写地理空间栅格数据，包括gdalnumeric模块，这样我们就可以将栅格数据加载到 NumPy 数组中。ASCII Grid 允许我们仅使用 Python 甚至 NumPy 读写栅格，因为它是简单的纯文本。

As a reminder, some elevation datasets use image formats to store elevation data. Most image formats only support 8-bit values ranging from between 0 to 255; however, some formats, including TIFF, can store larger values.

Geospatial software can typically display these datasets; however, traditional image software and libraries usually don't. For simplicity, in this chapter, we'll mostly stick to the ASCII Grid format for data, which is both human and machine-readable, as well as widely supported.

阅读网格

NumPy 能够使用其loadtxt()方法直接读取 ASCII Grid 格式，该方法旨在从文本文件中读取数组。前六行由标头组成，标头不是数组的一部分。以下几行是网格标题的示例:

ncols 250
nrows 250
xllcorner 277750.0
yllcorner 6122250.0
cellsize 1.0
NODATA_value -9999

让我们看看前面代码中的每一行包含什么:

第 1 行包含网格中的列数，与 x 轴同义。
第 2 行将 y 轴表示为行数。
第 3 行表示左下角的 x 坐标，是以米为单位的最小 x 值。
第 4 行是网格左下角对应的最小 y 值。
第 5 行是栅格的像元大小或分辨率。由于像元是正方形的，因此只需要一个大小值，而大多数地理空间栅格中的分辨率值是分开的 x 和 y 。
第 6 行是NODATA_value，这是一个分配给任何没有提供值的单元格的数字。

地理空间软件在计算时会忽略这些单元格，并且通常允许对其进行特殊的显示设置，例如将其设为黑色或透明。-9999值是一个常见的无数据占位符值，在行业中使用，在软件中很容易检测到，但可以任意选择。例如，负值的高程(即水深测量)可能在-9999米处具有有效数据，并且可能选择9999或其他值。只要在标题中定义了这个值，大多数软件都不会有问题。在一些例子中，我们将使用数字零；然而，零通常也可以是有效的数据值。

numpy.loadtxt()方法包括一个名为skiprows的参数，它允许您在读取数组值之前指定文件中要跳过的行数。

To try out this technique, you can download a sample grid file called myGrid.asc from git.io/vYapU.

因此，对于myGrid.asc，我们将使用以下代码:

myArray  = numpy.loadtxt("myGrid.asc", skiprows=6)

这一行导致myArray变量包含一个从 ascigridmyGrid.asc文件派生的numpy数组。ASC 文件扩展名由 ASCIIGRID 格式使用。这段代码很好用，但是有一个问题。NumPy 允许我们跳过标题，但不能保留它。我们需要保持这一点，以便我们的数据有一个空间参考。我们还将使用它来保存这个网格或创建一个新的网格。

为了解决这个问题，我们将使用 Python 内置的linecache模块来抓取头部。我们可以打开文件，遍历这些行，将每一行存储在一个变量中，然后关闭文件。但是，linecache将解决方案简化为一行。下面一行将文件的第一行读入一个名为line1的变量:

import linecache
line1 = linecache.getline("myGrid.asc", 1)

在本章的示例中，我们将使用这种技术来创建一个简单的头处理器，它可以用几行代码将这些头解析成 Python 变量。现在我们知道如何阅读网格，让我们学习如何编写它们。

书写网格

用 NumPy 写网格和读网格一样简单。我们使用相应的numpy.savetxt()函数将一个网格保存到一个文本文件中。唯一的问题是，在将数组转储到文件之前，我们必须构建并添加六行标题信息。对于不同版本的 NumPy，此过程略有不同。在这两种情况下，首先将头构建为字符串。如果您使用的是 NumPy 1.7 或更高版本，savetext()方法有一个名为 header 的可选参数，允许您指定一个字符串作为参数。您可以使用以下命令从命令行快速检查您的 NumPy 版本:

python -c "import numpy;print(numpy.__version__)"
1.8.2

向后兼容的方法是打开一个文件，写入头，然后转储数组。下面是将名为myArray的数组保存到名为myGrid.asc的 ASCIIGRID 文件的 1.7 版本方法示例:

header = "ncols {}\n".format(myArray.shape[1])
header += "nrows {}\n".format(myArray.shape[0])
header += "xllcorner 277750.0\n"
header += "yllcorner 6122250.0\n"
header += "cellsize 1.0\n"
header += "NODATA_value -9999"
numpy.savetxt("myGrid.asc", myArray, header=header, fmt="%1.2f")

我们使用 Python 格式字符串，这允许您在字符串中放置占位符来格式化要插入的 Python 对象。{}格式变量将您引用的对象转换为字符串。在这种情况下，我们引用数组中的列数和行数。

在 NumPy 中，数组有两个属性:

大小:返回一个整数作为数组中值的个数。
Shape:它返回一个元组，分别包含行数和列数。

因此，在前面的示例中，我们使用 shape 属性 tuple 将行数和列数添加到 ASCII 网格的标题中。请注意，我们还为每一行添加了一个尾随换行符(\n)。没有理由更改x和y值、单元格大小或无数据值，除非我们在脚本中更改它们。

savetxt()方法还有一个fmt参数，它允许您使用 Python 格式字符串来指定数组值的写入方式。在这种情况下，%1.2f值指定至少有一个数字且不超过两位小数的浮点数。1.6 之前的 NumPy 向后兼容版本以相同的方式构建头字符串，但首先创建文件句柄:

with open("myGrid.asc", "w") as f:
 f.write(header)
 numpy.savetxt(f, str(myArray), fmt="%1.2f")

正如您将在接下来的示例中看到的，这种仅使用 NumPy 生成有效地理空间数据文件的能力非常强大。在接下来的几个例子中，我们将使用加拿大不列颠哥伦比亚省温哥华附近山区的 ascigrid数字高程模型 ( DEM )。

You can download this sample as a ZIP file at the following URL: git.io/vYwUX.

下图是使用 QGIS 着色的原始数字高程模型，其色带使较低的高程值为深蓝色，较高的高程值为鲜红色:

虽然我们可以用这种方式从概念上理解数据，但这不是一种直观的数据可视化方式。让我们看看是否可以通过创建阴影浮雕做得更好。

创建阴影浮雕

阴影地形图以这样一种方式对高程进行着色，即地形看起来像是在低角度光线下投射的，这会产生亮点和阴影。这种美学风格创造了一种近乎摄影的错觉，很容易掌握，这样我们就可以理解地形的变化。需要注意的是，这种风格确实是一种错觉，因为从太阳角度来看，光线通常在物理上是不准确的，并且仰角通常被夸大以增加对比度。

在本例中，我们将使用之前引用的 ASCII 数字高程模型来创建另一个网格，该网格表示 NumPy 中地形的阴影浮雕版本。这个地形相当动态，所以我们不需要夸大海拔；但是，该脚本有一个名为z的变量，可以从 1.0 开始增加，以放大高程。

在我们定义了包括输入和输出文件名在内的所有变量之后，我们将看到基于linecache模块的标题解析器，它也使用 Python 列表理解来循环和解析行，然后将这些行从列表中拆分成六个变量。我们还创建了一个名为ycell的y细胞大小，它与常规细胞大小正好相反。如果我们不这样做，得到的网格将被转置。

Note that we define filenames for slope and aspect grids, which are two intermediate products that are combined to create the final product. These intermediate grids are output as well. They can also serve as inputs to other types of products.

该脚本使用三乘三的窗口方法来扫描图像，并平滑这些小网格中的中心值，以有效地处理图像。它是在计算机内存限制的情况下完成的。然而，因为我们使用的是 NumPy，所以我们可以通过矩阵一次处理整个数组，而不是使用一系列冗长的嵌套循环。这项技术是基于一位名叫 Michal Migurski 的开发人员的出色工作，他实现了马修·派瑞 C++实现的智能 NumPy 版本，该版本是 GDAL 套件中 DEM 工具的基础。

计算完坡度和坡向后，它们将用于输出着色地形。坡度是小山或山的陡度，而坡向是网格单元面向的方向，指定为 0 到 360 度之间的角度。最后，所有内容都从 NumPy 保存到磁盘。在savetxt()方法中，我们指定一个四整数格式的字符串，因为峰值海拔是几千米:

首先，我们将导入linecache模块来解析头部，导入numpy模块来进行处理:

from linecache import getline
import numpy as np

接下来，我们将设置定义如何处理着色浮雕的所有变量名:

# File name of ASCII digital elevation model
source = "dem.asc"

# File name of the slope grid
slopegrid = "slope.asc"

# File name of the aspect grid
aspectgrid = "aspect.asc"

# Output file name for shaded relief
shadegrid = "relief.asc"

# Shaded elevation parameters
# Sun direction
azimuth = 315.0

# Sun angle
altitude = 45.0

# Elevation exageration
z = 1.0

# Resolution
scale = 1.0

# No data value for output
NODATA = -9999

# Needed for numpy conversions
deg2rad = 3.141592653589793 / 180.0
rad2deg = 180.0 / 3.141592653589793

现在我们的变量已经设置好了，我们可以解析标题:

# Parse the header using a loop and
# the built-in linecache module
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(h.split(" ")[-1].strip()) for h in hdr]
cols, rows, lx, ly, cell, nd = values
xres = cell
yres = cell * -1

接下来，我们可以通过跳过标题部分使用numpy加载实际数据:

# Load the dem into a numpy array
arr = np.loadtxt(source, skiprows=6)

我们将逐行、逐列地循环数据，来处理它。但是，请注意，我们将跳过包含 nodata 值的外边缘。我们将把数据分成 3×3 像素的小网格，因为对于每个网格单元，我们需要看到它周围的单元:

# Exclude 2 pixels around the edges which are usually NODATA.
# Also set up structure for 3 x 3 windows to process the slope
# throughout the grid
window = []
for row in range(3):
 for col in range(3):
 window.append(arr[row:(row + arr.shape[0] - 2),
 col:(col + arr.shape[1] - 2)])

# Process each 3x3 window in both the x and y directions
x = ((z * window[0] + z * window[3] + z * window[3] + z * 
 window[6]) -
 (z * window[2] + z * window[5] + z * window[5] + z * 
 window[8])) / \
 (8.0 * xres * scale)
y = ((z * window[6] + z * window[7] + z * window[7] + z * 
 window[8]) -
 (z * window[0] + z * window[1] + z * window[1] + z * 
 window[2])) / \
 (8.0 * yres * scale)

对于每个 3×3 的迷你窗口，我们将计算slope、aspect，然后计算shaded浮雕值:

# Calculate slope
slope = 90.0 - np.arctan(np.sqrt(x * x + y * y)) * rad2deg

# Calculate aspect
aspect = np.arctan2(x, y)

# Calculate the shaded relief
shaded = np.sin(altitude * deg2rad) * np.sin(slope * deg2rad) + \
 np.cos(altitude * deg2rad) * np.cos(slope * deg2rad) * \
 np.cos((azimuth - 90.0) * deg2rad - aspect)

接下来，我们需要在 0-255 之间缩放每个值，以便可以将其视为图像:

# Scale values from 0-1 to 0-255
shaded = shaded * 255

现在，我们必须重建头部，因为我们忽略了 nodata 值的外边缘，并且数据集较小:

# Rebuild the new header
header = "ncols {}\n".format(shaded.shape[1])
header += "nrows {}\n".format(shaded.shape[0])
header += "xllcorner {}\n".format(lx + (cell * (cols - 
 shaded.shape[1])))
header += "yllcorner {}\n".format(ly + (cell * (rows - 
 shaded.shape[0])))
header += "cellsize {}\n".format(cell)
header += "NODATA_value {}\n".format(NODATA)

接下来，我们将把任何 nodata 值设置为我们在开始时在变量中设置的所选 nodata 值:


# Set no-data values
for pane in window:
 slope[pane == nd] = NODATA
 aspect[pane == nd] = NODATA
 shaded[pane == nd] = NODATA

我们将分别保存坡度网格和坡向网格，以便以后查看它们，并了解如何创建着色浮雕:

# Open the output file, add the header, save the slope grid
with open(slopegrid, "wb") as f:
 f.write(bytes(header, "UTF-8")
 np.savetxt(f, slope, fmt="%4i")

# Open the output file, add the header, save the aspectgrid
with open(aspectgrid, "wb") as f:
 f.write(bytes(header, "UTF-8")
 np.savetxt(f, aspect, fmt="%4i")

# Open the output file, add the header, save the relief grid
with open(shadegrid, "wb") as f:
 f.write(bytes(header, 'UTF-8'))
 np.savetxt(f, shaded, fmt="%4i")

如果我们将输出着色浮雕网格加载到 QGIS 并指定样式以将图像拉伸到最小值和最大值，我们将看到以下图像:

如果 QGIS 问你要投影，数据是 EPSG:3157。您也可以在我们在第 4 章、地理空间 Python 工具箱的安装 GDAL 部分讨论的 FWTools OpenEV 应用程序中打开图像，该应用程序将自动拉伸图像以获得最佳查看效果。

如您所见，前面的图像比我们最初检查的伪彩色表示更容易理解。接下来，让我们看看用于创建着色浮雕的坡度栅格:

斜率显示了数据集所有方向上从高点到低点的高程逐渐下降。对于许多类型的水文模型，坡度是一个特别有用的输入:

坡向显示了从一个像元到其相邻像元的最大下坡变化率。如果将纵横比图像与着色浮雕图像进行比较，您将看到纵横比图像的红色和灰色值对应于着色浮雕中的阴影。因此，坡度主要负责将数字高程模型转换为地形起伏，而坡向负责着色。

现在我们可以用一种有用的方式显示数据，让我们看看我们是否也可以从中创建其他数据。

创建高程等高线

等高线是数据集中沿同一高程的等值线。等高线通常以一定的间隔步进，以创建一种直观的方式来表示高程数据，包括视觉和数字，使用资源高效的向量数据集。现在，让我们看看另一种使用等高线更好地可视化高程的方法。

输入用于在我们的数字高程模型中生成等高线，输出是一个形状文件。用于生成轮廓的算法(行进的正方形:en.wikipedia.org/wiki/Marchi…)相当复杂，并且很难使用 NumPy 的线性代数来实现。在这种情况下，我们的解决方案是依靠 GDAL 库，它有一个通过 Python 应用编程接口可用的轮廓绘制方法。事实上，这个脚本的大部分只是设置输出 shapefile 所需的 OGR 库代码。实际的轮廓绘制是一个名为gdal.ContourGenerate()的单一方法调用。就在这个调用之前，有定义方法参数的注释。最重要的如下:

contourInterval:这是等高线之间的数据集单位距离。
contourBase:这是等高线的起始高程。
fixedLevelCount:相对于距离，这指定了固定数量的轮廓。
idField:这是一个必需的 shapefile dbf字段的名称，通常只叫 ID。
elevField:这是高程值所需的 shapefile dbf字段的名称，对于地图中的标注非常有用。

你应该从第四章地理空间 Python 工具箱的安装 GDAL 部分安装 GDAL 和 OGR。我们将实施以下步骤:

首先，我们将定义输入的数字高程模型文件名。
然后，我们将输出 shapefile 的名称。
接下来，我们将使用 OGR 创建 shapefile 数据源。
然后，我们将得到 OGR 层。
接下来，我们将打开数字高程模型。
最后，我们将在 OGR 图层上生成等高线。

让我们看一下前面步骤的代码表示:

首先，我们加载gdal和ogr库来处理数据:

import gdal
import ogr

然后我们将为文件名设置一个变量:

# Elevation DEM
source = "dem.asc"

接下来，我们将使用 OGR 创建输出 shapefile 的开头:

# Output shapefile
target = "contour"
ogr_driver = ogr.GetDriverByName("ESRI Shapefile")
ogr_ds = ogr_driver.CreateDataSource(target + ".shp")
ogr_lyr = ogr_ds.CreateLayer(target, 
# wkbLineString25D is the type code for geometry with a z 
# elevation value.
geom_type=ogr.wkbLineString25D)
field_defn = ogr.FieldDefn("ID" ogr.OFTInteger)
ogr_lyr.CreateField(field_defn)
field_defn = ogr.FieldDefn("ELEV" ogr.OFTReal)
ogr_lyr.CreateField(field_defn)

然后，我们将创建一些轮廓:

# gdal.ContourGenerate() arguments
# Band srcBand,
# double contourInterval,
# double contourBase,
# double[] fixedLevelCount,
# int useNoData,
# double noDataValue,
# Layer dstLayer,
# int idField,
# int elevField
ds = gdal.Open(source)

# EPGS:3157
gdal.ContourGenerate(ds.GetRasterBand(1), 400, 10, [], 0, 0, ogr_lyr, 0, 1))
ogr_ds = None

现在，让我们使用pngcanvas绘制刚刚创建的轮廓形状文件，这是我们在第 4 章地理空间 Python 工具箱的 PNGCanvas 部分中介绍的:

import shapefile
import pngcanvas

# Open the contours
r = shapefile.Reader("contour.shp")

# Setup the world to pixels conversion
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 800
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
contours = []

# Loop through all shapes
for shape in r.shapes():
 # Loop through all parts
 for i in range(len(shape.parts)):
   pixels = []
   pt = None
   if i < len(shape.parts) - 1:
     pt = shape.points[shape.parts[i]:shape.parts[i+1]]
   else:
     pt = shape.points[shape.parts[i]:]
   for x, y in pt:
     px = int(iwidth - ((r.bbox[2] - x) * xratio))
     py = int((r.bbox[3] - y) * yratio)
     pixels.append([px, py])
     contours.append(pixels)

# Set up the output canvas
canvas = pngcanvas.PNGCanvas(iwidth, iheight)

# PNGCanvas accepts rgba byte arrays for colors
red = [0xff, 0, 0, 0xff]
canvas.color = red

# Loop through the polygons and draw them
for c in contours:
 canvas.polyline(c)

# Save the image
with open("contours.png", "wb") as f:
 f.write(canvas.dump())

我们将以下面的图像结束:

如果我们将我们的着色浮雕 ASCIIGRID 和 shapefile 带入一个 GIS，比如 QGIS，我们可以创建一个简单的地形图，如下所示。您可以使用在脚本中指定的高程(即ELEV ) dbf字段，用高程标记等高线:

这些数字网格示例中使用的技术为各种高程产品提供了构建模块。接下来，我们将使用最复杂的高程数据类型之一:LIDAR 数据。

使用激光雷达数据

LIDAR 代表光探测和测距。它类似于基于雷达的图像，但使用每秒钟撞击地面数十万次的有限激光束来收集大量非常精细的( x 、 y 、 z )位置以及时间和强度。强度值是激光雷达与其他数据类型的真正区别。例如，建筑物的沥青屋顶可能与附近的树顶高度相同，但强度会有所不同。就像遥感一样，多光谱卫星图像中的辐射值允许我们建立分类库。激光雷达数据的强度值允许我们对激光雷达数据进行分类和着色。

激光雷达的高体积和高精度实际上使其难以使用。激光雷达数据集被称为点云，因为数据集的形状通常是不规则的，因为数据是带有外围点的三维数据。没有多少软件包可以有效地可视化点云。

此外，不规则形状的有限点集合很难交互，即使我们使用适当的软件也是如此。

由于这些原因，对激光雷达数据最常见的操作之一是投影数据并将其重新采样到常规网格中。我们将使用一个小型激光雷达数据集来完成这项工作。该数据集大约 7 MB 未压缩，包含超过 600，000 个点。这些数据捕捉了一些容易识别的特征，如建筑物、树木和停车场中的汽车。您可以从git.io/vOERW下载压缩数据集。

文件格式是激光雷达特有的一种非常常见的二进制格式，称为 LAS ，是激光的缩写。将此文件解压缩到您的工作目录。为了阅读这种格式，我们将使用一个名为laspy的纯 Python 库。您可以使用以下命令安装 Python 版:

pip install http://git.io/vOER9

安装laspy后，我们准备从 LIDAR 创建一个网格。

从激光雷达数据创建网格

这个脚本相当简单。我们循环遍历 LIDAR 数据中的( x 、 y )点位置，并将它们投影到一米大小的网格中。由于激光雷达数据的精度，我们将在一个单元中有多个点。我们对这些点进行平均，以创建一个共同的高程值。我们必须处理的另一个问题是数据丢失。无论何时对数据进行重新采样，都会丢失信息。

在这种情况下，我们将在光栅的中间出现NODATA孔。为了解决这个问题，我们用周围单元格的平均值填充这些洞，这是一种插值形式。我们只需要两个模块，都可以在 PyPI 上使用，如下面的代码所示:

from laspy.file import File
import numpy as np

# Source LAS file
source = "lidar.las"

# Output ASCII DEM file
target = "lidar.asc"

# Grid cell size (data units)
cell = 1.0

# No data value for output DEM
NODATA = 0

# Open LIDAR LAS file
las = File(source, mode="r")

# xyz min and max
min = las.header.min
max = las.header.max

# Get the x axis distance in meters
xdist = max[0] - min[0]

# Get the y axis distance in meters
ydist = max[1] - min[1]

# Number of columns for our grid
cols = int(xdist) / cell

# Number of rows for our grid
rows = int(ydist) / cell
cols += 1
rows += 1

# Track how many elevation
# values we aggregate
count = np.zeros((rows, cols)).astype(np.float32)

# Aggregate elevation values
zsum = np.zeros((rows, cols)).astype(np.float32)

# Y resolution is negative
ycell = -1 * cell

# Project x, y values to grid
projx = (las.x - min[0]) / cell
projy = (las.y - min[1]) / ycell

# Cast to integers and clip for use as index
ix = projx.astype(np.int32)
iy = projy.astype(np.int32)

# Loop through x, y, z arrays, add to grid shape,
# and aggregate values for averaging
for x, y, z in np.nditer([ix, iy, las.z]):
 count[y, x] += 1
 zsum[y, x] += z

# Change 0 values to 1 to avoid numpy warnings,
# and NaN values in array
nonzero = np.where(count > 0, count, 1)

# Average our z values
zavg = zsum / nonzero

# Interpolate 0 values in array to avoid any
# holes in the grid
mean = np.ones((rows, cols)) * np.mean(zavg)
left = np.roll(zavg, -1, 1)
lavg = np.where(left > 0, left, mean)
right = np.roll(zavg, 1, 1)
ravg = np.where(right > 0, right, mean)
interpolate = (lavg + ravg) / 2
fill = np.where(zavg > 0, zavg, interpolate)

# Create our ASCII DEM header
header = "ncols {}\n".format(fill.shape[1])
header += "nrows {}\n".format(fill.shape[0])
header += "xllcorner {}\n".format(min[0])
header += "yllcorner {}\n".format(min[1])
header += "cellsize {}\n".format(cell)
header += "NODATA_value {}\n".format(NODATA)

# Open the output file, add the header, save the array
with open(target, "wb") as f:
 f.write(bytes(header, 'UTF-8'))
 # The fmt string ensures we output floats
 # that have at least one number but only
 # two decimal places
 np.savetxt(f, fill, fmt="%1.2f")

我们的脚本的结果是一个 ASCIIGRID，当在 OpenEV 中查看时，它看起来像下图。海拔越高越亮，而海拔越低越暗。即使以这种形式，你也能看到建筑物、树木和汽车:

如果我们指定了热图颜色渐变，颜色会让您更清晰地感受到高程差异:

那么，如果我们通过之前的着色地形脚本运行这个输出 DEM，会发生什么呢？直边建筑和斜山有很大的区别。如果您更改着色地形脚本中的输入和输出名称来处理 LIDAR DEM，我们会得到以下坡度结果:

山区平缓起伏的坡度被简化为图像中主要特征的轮廓。在方面图像中，变化非常明显，并且距离如此之短，以至于输出图像看起来非常混乱，如下图所示:

尽管这些图像与粗糙但有些平滑的山地版本有所不同，但我们仍然获得了非常好的阴影浮雕，在视觉上类似于黑白照片:

现在我们知道如何处理 LIDAR 数据，让我们学习如何使用 Python 可视化它。

利用 PIL 实现激光雷达数据可视化

本章前面的数字高程模型图像是使用 QGIS 和 OpenEV 可视化的。我们还可以通过引入前几章没有用到的 Python 图像库 ( PIL )的一些新功能，用 Python 创建输出图像。

在这个例子中，我们将使用PIL.ImageOps模块，它具有直方图均衡化和自动对比度增强的功能。我们将使用 PIL 的fromarray()方法从numpy导入数据。让我们看看在以下代码的帮助下，我们可以多么接近本章中所示的桌面 GIS 程序的输出:

import numpy as np

try:
 import Image
 import ImageOps
except ImportError:
 from PIL import Image, ImageOps

# Source gridded LIDAR DEM file
source = "lidar.asc"

# Output image file
target = "lidar.bmp"

# Load the ASCII DEM into a numpy array
arr = np.loadtxt(source, skiprows=6)

# Convert array to numpy image
im = Image.fromarray(arr).convert("RGB")

# Enhance the image:
# equalize and increase contrast
im = ImageOps.equalize(im)
im = ImageOps.autocontrast(im)

# Save the image
im.save(target)

如您所见，在下图中，增强的阴影浮雕比以前的版本具有更清晰的浮雕:

现在，让我们为我们的阴影浮雕上色。我们将使用内置的 Python colorsys模块进行颜色空间转换。通常，我们将颜色指定为 RGB 值。然而，为了创建热图方案的颜色渐变，我们将使用 HSV (色调、饱和度和值的缩写)值来生成我们的颜色。

HSV 的优点是可以在色轮上将 H 值调整到 0 到 360 度之间。使用单一色调值允许您使用线性渐变方程，这比尝试处理三个独立 RGB 值的组合要容易得多。下图来自网络杂志 Qt 季刊，展示了 HSV 的颜色模型:

colorsys模块允许您在 HSV 和 RGB 值之间来回切换。该模块返回 RGB 值的百分比，然后必须将其映射到每种颜色的 0-255 比例。

在下面的代码中，我们将把 ASCII 数字高程模型转换成 PIL 图像，构建调色板，将调色板应用于灰度图像，并保存图像:

import numpy as np

try:
 import Image
 import ImageOps
except:
 from PIL import Image, ImageOps
import colorsys

# Source LIDAR DEM file
source = "lidar.asc"

# Output image file
target = "lidar.bmp"

# Load the ASCII DEM into a numpy array
arr = np.loadtxt(source, skiprows=6)

# Convert the numpy array to a PIL image.
# Use black and white mode so we can stack
# three bands for the color image.
im = Image.fromarray(arr).convert('L')

# Enhance the image
im = ImageOps.equalize(im)
im = ImageOps.autocontrast(im)

# Begin building our color ramp
palette = []

# Hue, Saturation, Value
# color space starting with yellow.
h = .67
s = 1
v = 1

# We'll step through colors from:
# blue-green-yellow-orange-red.
# Blue=low elevation, Red=high-elevation
step = h / 256.0

# Build the palette
for i in range(256):
 rp, gp, bp = colorsys.hsv_to_rgb(h, s, v)
 r = int(rp * 255)
 g = int(gp * 255)
 b = int(bp * 255)
 palette.extend([r, g, b])
 h -= step

# Apply the palette to the image
im.putpalette(palette)

# Save the image
im.save(target)

上面的代码生成了下面的图像，较高的高度用较暖的颜色，较低的高度用较冷的颜色:

在这张图片中，我们实际上比默认的 QGIS 版本得到了更多的变化。我们可以潜在地通过平滑算法来改善这张图像，该算法可以在颜色相遇的地方混合颜色，并在视觉上柔化图像。如您所见，随着海拔变化的增加，我们的颜色渐变范围从冷色到暖色。

创建不规则三角网

下面的例子是我们迄今为止最复杂的例子。一个不规则三角网 ( 三角网)是一个点数据集的向量表示，该点数据集位于以三角形连接的点的向量表面中。一种算法确定精确表示地形绝对需要哪些点，这与栅格相反，栅格在给定区域内存储固定数量的像元，并且可以在相邻像元中重复高程值，这样可以更有效地存储为多边形。

在地理信息系统中使用三角网时，三角网也可以比栅格更有效地进行动态重采样，栅格需要更少的计算机内存和处理能力。最常见的 TIN 类型是基于 Delaunay 三角剖分，包括所有没有冗余三角形的点。

德劳奈三角测量非常复杂。我们将使用一个由比尔·西蒙斯编写的纯 Python 库，作为史蒂夫·财富的德劳奈三角测量算法voronoi.py的一部分，来计算我们的激光雷达数据中的三角形。您可以从git.io/vOEuJ下载脚本到您的工作目录或site-packages目录。

这个脚本读取 LAS 文件，生成三角形，循环遍历它们，并写出一个形状文件。对于本例，我们将使用激光雷达数据的剪辑版本来减少要处理的区域。如果我们运行 600，000+点的整个数据集，脚本将运行数小时，并生成超过 50 万个三角形。您可以从以下网址下载剪辑的激光雷达数据集作为压缩文件:git.io/vOE62。

由于以下示例的时间密集性，可能需要几分钟才能完成，因此脚本运行时会打印几条状态消息。我们将把三角形存储为多边形类型，这允许顶点具有z高程值。解压 LAS 文件，运行以下代码生成名为mesh.shp的 shapefile:

首先，我们导入我们的库:

import pickle
import os
import time
import math
import numpy as np
import shapefile
from laspy.file import File
# voronoi.py for Python 3: pip install http://git.io/vOEuJ
import voronoi

接下来，我们定义我们的 LIDAR 文件、目标输出文件和 pickle 文件的位置和名称:

# Source LAS file
source = "clippedLAS.las"

# Output shapefile
target = "mesh"

# Triangles pickle archive
archive = "triangles.p"

现在，我们将创建voronoi模块所需的点类:

class Point:
 """Point class required by the voronoi module"""
 def __init__(self, x, y):
   self.px = x
   self.py = y

def x(self):
 return self.px

def y(self):
 return self.py

接下来，我们将创建一个三角形数组来跟踪为网格创建的三角形:

# The triangle array holds tuples
# 3 point indices used to retrieve the points.
# Load it from a pickle
# file or use the voronoi module
# to create the triangles.
triangles = None

接下来，我们需要打开我们的 LIDAR 文件并提取点:

 # Open LIDAR LAS file
 las = File(source, mode="r")
else:

# Open LIDAR LAS file
 las = File(source, mode="r")
 points = []
 print("Assembling points...")

# Pull points from LAS file
 for x, y in np.nditer((las.x, las.y)):
 points.append(Point(x, y))
 print("Composing triangles...")

现在，我们可以对这些点执行德劳奈计算来构建三角形:

# Delaunay Triangulation
 triangles = voronoi.computeDelaunayTriangulation(points)

如果我们再次运行这个精确的脚本，我们会将三角形转储到 pickle 归档中以节省时间:

 # Save the triangles to save time if we write more than
 # one shapefile.
 f = open(archive, "wb")
 pickle.dump(triangles, f, protocol=2)
 f.close()

接下来，我们可以创建一个 shapefile Writer对象，通过设置必要的字段来开始创建我们的输出 shapefile:

print("Creating shapefile...")
 # PolygonZ shapefile (x, y, z, m)
 w = shapefile.Writer(target, shapefile.POLYGONZ)
 w.field("X1", "C", "40")
 w.field("X2", "C", "40")
 w.field("X3", "C", "40")
 w.field("Y1", "C", "40")
 w.field("Y2", "C", "40")
 w.field("Y3", "C", "40")
 w.field("Z1", "C", "40")
 w.field("Z2", "C", "40")
 w.field("Z3", "C", "40")
 tris = len(triangles)

然后，我们循环通过三角形并创建网格:

# Loop through shapes and
 # track progress every 10 percent
 last_percent = 0
 for i in range(tris):
     t = triangles[i]
     percent = int((i/(tris*1.0))*100.0)
     if percent % 10.0 == 0 and percent > last_percent:
         last_percent = percent
         print("{} % done - Shape {}/{} at {}".format(percent, 
         i, tris, time.asctime()))
 part = []
 x1 = las.x[t[0]]
 y1 = las.y[t[0]]
 z1 = las.z[t[0]]
 x2 = las.x[t[1]]
 y2 = las.y[t[1]]
 z2 = las.z[t[1]]
 x3 = las.x[t[2]]
 y3 = las.y[t[2]]
 z3 = las.z[t[2]]

接下来，我们可以消除任何非常长的线段，这些线段是库的误判:

 # Check segments for large triangles
 # along the convex hull which is a common
 # artifact in Delaunay triangulation
 max = 3
 if math.sqrt((x2-x1)**2+(y2-y1)**2) > max:
 continue
 if math.sqrt((x3-x2)**2+(y3-y2)**2) > max:
 continue
 if math.sqrt((x3-x1)**2+(y3-y1)**2) > max:
 continue
 part.append([x1, y1, z1, 0])
 part.append([x2, y2, z2, 0])
 part.append([x3, y3, z3, 0])
 w.poly(parts=[part])
 w.record(x1, x2, x3, y1, y2, y3, z1, z2, z3)
 print("Saving shapefile...")

最后，我们可以保存输出 shapefile:

w.close()
print("Done.")

下图显示了彩色激光雷达数据上 TIN 的放大版本:

网格从点云提供了一个高效、连续的表面，这比点云本身更容易处理。

摘要

高程数据通常可以为分析和衍生产品提供完整的数据集，而无需任何其他数据。在本章中，您学习了如何仅使用 NumPy 读写 ASCII 网格。您还学习了如何创建着色浮雕、坡度网格和坡向网格。我们使用一个鲜为人知的特性创建了高程等高线，这个特性叫做 GDAL 库的等高线，Python 提供了这个特性。

接下来，我们将激光雷达数据转换成易于操作的 ASCII 网格。我们用不同的方法对 PIL 的激光雷达数据进行了可视化实验。最后，我们通过将激光雷达点云转换为多边形的三维形状文件来创建三维表面或三角网。这些是地形分析工具，用于交通规划、建筑规划、水文排水建模、地质勘探等。

在下一章中，我们将结合前三章的构建块来执行一些高级建模，并实际创建一些信息产品。

进一步阅读

您可以通过以下链接找到一些关于 Python 和高程数据的附加教程:https://www . earth data science . org/tutories/Python/elevation/。

八、高级地理空间 Python 建模

在这一章中，我们将在已经学习的数据处理概念的基础上，创建一些全面的信息产品。以前介绍的数据处理方法很少能自己提供问题的答案。您可以结合这些数据处理方法，从多个已处理的数据集构建地理空间模型。地理空间模型是现实世界某个方面的简化表示，有助于我们回答关于项目或问题的一个或多个问题。在本章中，我们将介绍一些在农业、应急管理、物流和其他行业中常用的重要地理空间算法。

我们将创造的产品如下:

作物健康地图
洪水淹没模型
彩色的山影
地形路线图
街道路线图
带有地理位置照片链接的形状文件

虽然这些产品是特定于任务的，但是用于创建它们的算法被广泛应用于地理空间分析。我们将在本章中讨论以下主题:

创建归一化差异植被指数
创建洪水淹没模型
创建彩色山体阴影
执行最低成本路径分析
将路线转换为形状文件
沿着街道走
地理定位照片
计算卫星图像云量

本章中的例子比前几章更长，涉及面更广。因此，有更多的代码注释来使程序更容易理解。在这些例子中，我们还将使用更多的函数。在前面的章节中，为了清晰起见，大部分都避免了函数，但是这些例子足够复杂，以至于某些函数使得代码更容易阅读。这些示例是作为地理空间分析师，您将在工作中使用的实际流程。

技术要求

对于本章，需要满足以下要求:

版本 : Python 3.6 或更高
RAM :最小 6 GB (Windows)，8gb(macOS)；推荐 8 GB
存储:最低 7,200 RPM SATA，可用空间 20 GB，推荐 SSD，可用空间 40 GB。
处理器:最低英特尔酷睿 i3 2.5 GHz，推荐英特尔酷睿 i5。

创建归一化差异植被指数

我们的第一个例子是归一化差异植被指数 ( NVDI )。NDVIs 用于显示感兴趣区域内植物的相对健康状况。一种 NDVI 算法使用卫星或航空图像，通过突出植物中的叶绿素密度来显示相对健康状况。NDVIs 只使用红色和近红外波段。NDVI 公式如下:

NDVI = (Infrared – Red) / (Infrared + Red)

这项分析的目标是，首先生成一个包含红外和红色波段的多光谱图像，最后生成一个使用七个类别的伪彩色图像，这些类别将健康的植物染成深绿色，不太健康的植物染成浅绿色，裸露的土壤染成棕色。

因为健康指数是相对的，所以定位感兴趣的区域很重要。您可以对整个地球执行相对索引，但是广阔的区域，例如低植被极端的撒哈拉沙漠和茂密的森林区域，例如亚马逊丛林，会扭曲中等范围内植被的结果。然而，话虽如此，气候科学家还是会定期创建全球 NDVIs 来研究全球趋势。不过，更常见的应用是用于管理区域，如森林或农田，如本例所示。

我们将首先分析密西西比三角洲的一个农场。为此，我们将从相当大面积的多光谱图像开始，并使用 shapefile 来隔离单个场。下面截图中的图像是我们的广阔区域，感兴趣的领域用黄色突出显示:

您可以从git.io/v3fS9下载该图像和农场区域的形状文件作为 ZIP 文件。

对于这个例子，我们将使用 GDAL，OGR，gdal_array / numpy，以及 Python 图像库 ( PIL )来剪辑和处理数据。在本章的其他示例中，我们将只使用简单的 ASCII 网格和 NumPy。因为我们将使用 ASCII 高程网格，所以不需要 GDAL。在所有示例中，脚本使用以下约定:

导入库。
定义函数。
定义全局变量，如文件名。
执行分析。
保存输出。

我们对作物健康示例的方法分为两个脚本。第一个脚本创建索引图像，这是一个灰度图像。第二个脚本对索引进行分类，并输出彩色图像。在第一个脚本中，我们将执行以下步骤来创建索引映像:

读取红外波段。
读取字段边界形状文件。
将形状文件光栅化为图像。
将 shapefile 图像转换为 NumPy 数组。
使用 NumPy 数组将红色带裁剪到字段中。
对红外波段也这样做。
使用波段数组在 NumPy 中执行 NDVI 算法。
使用gdal_array将生成的索引算法保存到一个 GeoTIFF 文件中。

我们将分部分讨论这个脚本，以便于理解。代码注释还会告诉你每一步都发生了什么。

建立框架

设置框架将帮助我们导入我们需要的模块，并设置我们将用于前面指令的步骤 1 到 5 的功能。imageToArray()功能将 PIL 图像转换为 NumPy 数组，并依赖于gdal_array和 PIL 模块。world2Pixel()功能将地理空间坐标转换为目标图像的像素坐标。该功能使用gdal模块提供的地理参考信息。copy_geo()函数将源图像中的地理参考信息复制到目标数组中，但会考虑到在裁剪图像时产生的偏移。这些功能相当通用，可以在本例以外的各种不同遥感过程中发挥作用:

首先，我们导入我们的库:

import gdal
from osgeo import gdal
from osgeo import gdal_array
from osgeo import ogr
try:
 import Image
 import ImageDraw
except ImportError:
 from PIL import Image, ImageDraw

然后，我们需要一个函数将图像转换成numpy数组:

def imageToArray(i):
    """
    Converts a Python Imaging Library
    array to a gdal_array image.
    """
    a = gdal_array.numpy.fromstring(i.tobytes(), 'b')
    a.shape = i.im.size[1], i.im.size[0]
    return a

现在，我们将设置一个函数，将坐标转换为图像像素:

def world2Pixel(geoMatrix, x, y):
 """
 Uses a gdal geomatrix (gdal.GetGeoTransform())
 to calculate the pixel location of a
 geospatial coordinate
 """
 ulX = geoMatrix[0]
 ulY = geoMatrix[3]
 xDist = geoMatrix[1]
 yDist = geoMatrix[5]
 rtnX = geoMatrix[2]
 rtnY = geoMatrix[4]
 pixel = int((x - ulX) / xDist)
 line = int((ulY - y) / abs(yDist))
 return (pixel, line)

最后，我们将创建一个从图像复制地理元数据的函数:

def copy_geo(array, prototype=None, xoffset=0, yoffset=0):
 """Copy geotransfrom from prototype dataset to array but account
 for x, y offset of clipped array."""
 ds = gdal_array.OpenArray(array)
 prototype = gdal.Open(prototype)
 gdal_array.CopyDatasetInfo(prototype, ds,
 xoff=xoffset, yoff=yoffset)
 return ds

下一步是加载数据，我们将在下一节中检查这些数据。

正在加载数据

在这一节中，我们使用gdal_array加载农场场的源图像，它直接进入一个 NumPy 数组。我们还定义了输出图像的名称，它将是ndvi.tif。这一部分的一个有趣的部分是，我们使用gdal模块第二次加载源图像，而不是gdal_array。

第二个调用是捕获通过gdal而不是gdal_array获得的图像的地理参考数据。幸运的是，gdal仅按需加载栅格数据，因此这种方法避免了将完整数据集加载到内存中两次。一旦我们有了多维 NumPy 数组的数据，我们就可以分离出红色和红外波段，因为它们都将用于 NDVI 方程:

# Multispectral image used
# to create the NDVI. Must
# have red and infrared
# bands
source = "farm.tif"

# Output geotiff file name
target = "ndvi.tif"

# Load the source data as a gdal_array array
srcArray = gdal_array.LoadFile(source)

# Also load as a gdal image to
# get geotransform info
srcImage = gdal.Open(source)
geoTrans = srcImage.GetGeoTransform()

# Red and infrared (or near infrared) bands
r = srcArray[1]
ir = srcArray[2]

现在我们已经加载了数据，我们可以将 shapefile 转换为光栅。

光栅化形状文件

本节开始剪辑过程。然而，第一步是光栅化 shapefile，它勾勒出我们将要分析的特定区域的边界。该区域在更大的field.tif卫星图像内。换句话说，我们将其从向量数据转换为栅格数据。但是我们也想在转换多边形时填充它，以便它可以用作图像遮罩。遮罩中的像素将与红色和红外数组中的像素相关联。

遮罩之外的任何像素都将变成NODATA像素，因此它们不会作为 NDVI 的一部分进行处理。为了进行这种关联，我们需要实心多边形成为一个 NumPy 数组，就像栅格波段一样。这种方法将确保我们的 NDVI 计算仅限于农田。

将 shapefile 多边形转换为 NumPy 数组形式的填充多边形的最简单方法是将其绘制为 PIL 图像中的多边形，填充该多边形，然后使用允许转换的 PIL 和 NumPy 中的现有方法将其转换为 NumPy 数组。

在这个例子中，我们使用ogr模块来读取 shapefile，因为我们已经有了可用的 GDAL。但是，我们也可以使用 PyShp 同样轻松地读取 shapefile。如果我们的农田图像可以作为 ASCII 网格使用，我们就可以完全避免使用gdal、gdal_array和ogr模块:

首先，我们打开 shapefile 并选择唯一的一层:

# Clip a field out of the bands using a
# field boundary shapefile

# Create an OGR layer from a Field boundary shapefile
field = ogr.Open("field.shp")
# Must define a "layer" to keep OGR happy
lyr = field.GetLayer("field")

只有一个多边形，所以我们要抓住这个特征:

# Only one polygon in this shapefile
poly = lyr.GetNextFeature()

现在我们将图层范围转换为图像像素坐标:

# Convert the layer extent to image pixel coordinates
minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

然后，我们计算新图像的像素大小:

# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

接下来，我们以正确的大小创建一个新的空白图像:

# Create a blank image of the correct size
# that will serve as our mask
clipped = gdal_array.numpy.zeros((3, pxHeight, pxWidth),
 gdal_array.numpy.uint8)

现在，我们准备使用边界框来裁剪红色和红外波段:

# Clip red and infrared to new bounds.
rClip = r[ulY:lrY, ulX:lrX]
irClip = ir[ulY:lrY, ulX:lrX]

接下来，我们为图像创建地理配准信息:

# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

然后，我们可以准备将点映射到像素，以创建我们的遮罩图像:

# Map points to pixels for drawing
# the field boundary on a blank
# 8-bit, black and white, mask image.
points = []
pixels = []
# Grab the polygon geometry
geom = poly.GetGeometryRef()
pts = geom.GetGeometryRef(0)

我们遍历所有点要素并存储它们的 x 和 y 值:

# Loop through geometry and turn
# the points into an easy-to-manage
# Python list
for p in range(pts.GetPointCount()):
    points.append((pts.GetX(p), pts.GetY(p)))

现在，我们将这些点转换为像素位置:

# Loop through the points and map to pixels.
# Append the pixels to a pixel list
for p in points:
    pixels.append(world2Pixel(geoTrans, p[0], p[1]))

接下来，我们创建一个新图像作为我们的蒙版图像:

# Create the raster polygon image as a black and white 'L' mode
# and filled as white. White=1
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)

现在我们可以光栅化我们的多边形:

# Create a PIL drawing object
rasterize = ImageDraw.Draw(rasterPoly)

# Dump the pixels to the image
# as a polygon. Black=0
rasterize.polygon(pixels, 0)

最后，我们可以将我们的遮罩转换成一个numpy数组:

# Hand the image back to gdal/gdal_array
# so we can use it as an array mask
mask = imageToArray(rasterPoly)

现在，我们已经将 shapefile 转换为遮罩图像，我们可以剪辑波段了。

修剪乐队

现在我们有了图像遮罩，我们可以将红色和红外波段裁剪到遮罩的边界。对于此过程，我们使用 NumPy 的choose()方法，该方法将遮罩像元与栅格波段像元相关联，并返回该值，或返回0。结果是一个新数组被裁剪到遮罩上，但具有栅格波段中的相关值:

# Clip the red band using the mask
rClip = gdal_array.numpy.choose(mask,
 (rClip, 0)).astype(gdal_array.numpy.uint8)

# Clip the infrared band using the mask
irClip = gdal_array.numpy.choose(mask,
 (irClip, 0)).astype(gdal_array.numpy.uint8)

我们现在只有我们想要的数据，所以我们可以应用我们的 NDVI 相对植被健康公式。

使用 NDVI 公式

我们创建 NDVI 的最后过程是执行方程式红外-红色/红外+红色。我们执行的第一步是消除任何非数字，也称为 NaN ，NumPy 中可能在除法运算中出现的值。在保存输出之前，我们将把任何 NaN 值转换成0。我们将输出保存为ndvi.tif，这将是下一个脚本的输入，以便对 NDVI 进行分类和着色，如下所示:

首先，我们将忽略来自numpy的任何警告，因为我们将在边缘附近得到一些错误:

# We don't care about numpy warnings
# due to NaN values from clipping
gdal_array.numpy.seterr(all="ignore")

现在我们可以执行我们的 NDVI 公式:

# NDVI equation: (infrared - red) / (infrared + red)
# *1.0 converts values to floats,
# +1.0 prevents ZeroDivisionErrors
ndvi = 1.0 * ((irClip - rClip) / (irClip + rClip + 1.0))

如果有任何 NaN 值，我们将其转换为零:

# Convert any NaN values to 0 from the final product
ndvi = gdal_array.numpy.nan_to_num(ndvi)

最后，我们保存我们完成的 NDVI 图像:

# Save the ndvi as a GeoTIFF and copy/adjust 
# the georeferencing info
gtiff = gdal.GetDriverByName( 'GTiff' )
gtiff.CreateCopy(target, copy_geo(ndvi, prototype=source, xoffset=ulX, yoffset=ulY))
gtiff = None

下图是该示例的输出。您需要在地理空间查看器(如 QGIS 或 OpenEV)中查看它。大多数图像编辑器都不会打开该图像。灰色阴影越浅，那块地里的植物就越健康:

现在我们知道如何使用 NDVI 公式，让我们看看如何对它进行分类。

对 NDVI 进行分类

我们现在有一个有效的索引，但是不容易理解，因为它是一个灰度图像。如果我们用直观的方式给图像上色，那么即使是一个孩子也能识别出更健康的植物。在下一节附加功能中，我们读入这个灰度索引，并使用七个类别将其从棕色分类为深绿色。分类和图像处理例程，例如直方图和拉伸功能，几乎与我们在第 6 章、 Python 和遥感中创建直方图部分使用的相同，但这次我们以更具体的方式应用它们。

这个例子的输出将是另一个 GeoTIFF 文件，但这次它将是一个彩色的 RGB 图像。

附加功能

我们不需要之前 NDVI 脚本中的任何函数，但是我们需要添加一个函数来创建和拉伸直方图。这两个函数都适用于 NumPy 数组。我们也将在这个脚本中缩短对gdal_array到gd的引用，因为这是一个很长的名称，我们在整个脚本中都需要它。

让我们看看下面的步骤:

首先，我们导入我们需要的库:

import gdal_array as gd
import operator
from functools import reduce

接下来，我们需要创建一个histogram函数，我们需要它来进行直方图拉伸:

def histogram(a, bins=list(range(256))):
 """
 Histogram function for multi-dimensional array.
 a = array
 bins = range of numbers to match
 """
 # Flatten, sort, then split our arrays for the histogram.
 fa = a.flat
 n = gd.numpy.searchsorted(gd.numpy.sort(fa), bins)
 n = gd.numpy.concatenate([n, [len(fa)]])
 hist = n[1:]-n[:-1]
 return hist

现在，我们创建直方图stretch函数:

def stretch(a):
 """
 Performs a histogram stretch on a gdal_array array image.
 """
 hist = histogram(a)
 lut = []
 for b in range(0, len(hist), 256):
 # step size – create equal interval bins.
 step = reduce(operator.add, hist[b:b+256]) / 255
 # create equalization lookup table
 n = 0
 for i in range(256):
 lut.append(n / step)
 n = n + hist[i+b]
 gd.numpy.take(lut, a, out=a)
 return a

既然我们有了实用功能，我们就可以处理 NDVI 了。

装载 NDVI

接下来，我们将把 NDVI 脚本的输出加载回 NumPy 数组。我们还将输出图像的名称定义为ndvi_color.tif，并创建一个零填充的多维数组作为彩色 NDVI 图像的红色、绿色和蓝色带的占位符。以下代码将把 NDVI TIFF 图像加载到一个numpy数组中:

# NDVI output from ndvi script
source = "ndvi.tif"

# Target file name for classified
# image image
target = "ndvi_color.tif"

# Load the image into an array
ndvi = gd.LoadFile(source).astype(gd.numpy.uint8)

现在我们的图像被加载为一个数组，我们可以拉伸它。

准备 NDVI

我们需要对 NDVI 进行直方图拉伸，以确保图像覆盖了赋予最终产品意义的类别范围:

# Peform a histogram stretch so we are able to
# use all of the classes
ndvi = stretch(ndvi)

# Create a blank 3-band image the same size as the ndvi
rgb = gd.numpy.zeros((3, len(ndvi), len(ndvi[0])), gd.numpy.uint8)

现在我们已经拉伸了图像，我们可以开始分类过程了。

创建类

在这一部分中，我们为我们的 NDVI 类设置了范围，这些范围从 0 到 255 不等。我们将使用七节课。您可以通过在类列表中添加或删除值来更改类的数量。接下来，我们创建一个查找表，或者 LUT ，以便为每个类分配颜色。颜色的数量必须与类别的数量相匹配。

颜色被定义为 RGB 值。start变量定义了第一个类的开始。在这种情况下，0是一个 nodata 值，这是我们在前面的脚本中指定的，所以我们在1开始上课。然后，我们遍历类，提取范围，并使用颜色分配将 RGB 值添加到占位符数组中。最后，我们将彩色图像保存为 GeoTIFF 文件:

# Class list with ndvi upper range values.
# Note the lower and upper values are listed on the ends
classes = [58, 73, 110, 147, 184, 220, 255]

# Color look-up table (lut)
# The lut must match the number of classes
# Specified as R, G, B tuples from dark brown to dark green
lut = [[120, 69, 25], [255, 178, 74], [255, 237, 166], [173, 232, 94],
 [135, 181, 64], [3, 156, 0], [1, 100, 0]]

# Starting value of the first class
start = 1

现在我们可以对图像进行分类:

# For each class value range, grab values within range,
# then filter values through the mask.
for i in range(len(classes)):
 mask = gd.numpy.logical_and(start <= ndvi,
 ndvi <= classes[i])
 for j in range(len(lut[i])):
     rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j]))
     start = classes[i]+1

最后，我们可以保存我们的分类 GeoTIFF 文件:

# Save a geotiff image of the colorized ndvi.
output=gd.SaveArray(rgb, target, format="GTiff", prototype=source)
output = None

这是我们输出的图像:

这是我们这个例子的最终产品。农民可以利用这些数据来确定如何以有针对性、更有效和更环保的方式有效灌溉和喷洒化肥和农药等化学品。事实上，这些类甚至可以转化为向量形状文件，然后加载到田间喷雾器上的全球定位系统驱动的计算机中。然后，当喷雾器在田地周围行驶时，或者在某些情况下，甚至在带有喷雾器附件的飞机上飞过田地时，这将自动在正确的位置应用正确数量的化学品。

还要注意，即使我们将数据剪切到字段中，图像仍然是一个正方形。黑色区域是已转换为黑色的 nodata 值。在显示软件中，您可以使 nodata 颜色透明，而不影响图像的其余部分。

虽然我们创建了一个非常特殊的产品类型，一个分类的 NDVI，这个脚本的框架可以被改变，以实现许多遥感分析算法。NDVIs 有不同的类型，但通过相对较小的更改，您可以将此脚本变成一个工具，用于查找海洋中的有害藻类水华，或者森林中间的烟雾，表明发生了森林火灾。

This book attempts to limit the use of GDAL as much as possible in order to focus on what can be accomplished with pure Python and tools that can easily be installed from PyPI. However, it is helpful to remember that there is a wealth of information on using GDAL and its associated utilities to carry out similar tasks. For another tutorial on clipping a raster with GDAL via its command-line utilities, see joeyklee.github.io/broc-cli-ge….

现在我们已经研究了土地，让我们研究水，以便创建洪水淹没模型。

创建洪水淹没模型

在下一个例子中，我们将开始进入水文学的世界。洪水是最常见和最具破坏性的自然灾害之一，影响到全球几乎每一个人口。地理空间模型是评估洪水影响并在洪水发生前减轻其影响的有力工具。我们经常在新闻上听到一条河流正在达到洪水阶段，但是如果我们不能理解其影响，这些信息就没有意义。

水文洪水模型的开发成本很高，而且可能非常复杂。这些模型对工程师建造防洪系统至关重要。然而，第一反应者和潜在的洪水受害者只对即将到来的洪水的影响感兴趣。

我们可以使用一个非常简单易懂的工具洪水淹没模型开始了解一个地区的洪水影响。这个模型从一个点开始，在一个特定的洪水阶段，用一个洪水盆地所能容纳的最大水量淹没一个区域。通常，这种分析是最坏的情况。成百上千的其他因素被用来计算有多少水将从河流顶部的洪水阶段进入流域。但是我们仍然可以从这个简单的一阶模型中学到很多。

As mentioned in the Elevation data section in Chapter 1, Learning about Geospatial Analysis with Python, the Shuttle Radar Topography Mission (SRTM) dataset provides a nearly-global DEM that you can use for these types of models. More on SRTM data can be found here: www2.jpl.nasa.gov/srtm/.

你可以从 git.io/v3fSg 下载 EPSG 的 ASCII 网格数据:4326 和一个包含该点的形状文件作为.zip文件。shapefile 仅供参考，在这个模型中没有任何作用。下图是一个数字高程模型 ( DEM )源点显示为德克萨斯州休斯顿附近的一颗黄色恒星。在现实世界的分析中，这个点可能是一个测流计，在这里你可以得到河流水位的数据:

我们在这个例子中介绍的算法叫做洪水填充算法。该算法在计算机科学领域众所周知，在经典计算机游戏扫雷中使用，当用户点击方块时清除棋盘上的空方块。也是众所周知的油漆桶工具在 Adobe Photoshop 等图形程序中使用的方法，用于用不同的颜色填充同一颜色相邻像素的区域。

有很多方法可以实现这个算法。最古老和最常见的方法之一是递归地遍历图像的每个像素。递归的问题在于，你最终会不止一次地处理像素，并产生不必要的工作量。递归泛洪填充的资源使用很容易使程序崩溃，即使是中等大小的图像。

该脚本使用基于队列的四路泛洪填充，可以多次访问一个单元格，但确保我们只处理一次单元格。通过使用 Python 的内置集合类型，队列只包含唯一的、未处理的单元，该集合类型只保存唯一的值。我们使用两组:填充，包含我们需要填充的单元格，填充，包含处理后的单元格。

本示例执行以下步骤:

从 ASCII 数字高程模型中提取标题信息。
打开数字高程模型作为numpy数组。
将我们的起点定义为数组中的行和列。
声明洪水高程值。
将地形过滤到所需的高程值及以下。
处理过滤后的数组。
创建一个 1，0，0 数组(即二进制数组)，泛洪像素为 1。
将洪水淹没数组保存为 ASCII 网格。

This example can take a minute or two to run on a slower machine; we'll use the print statements throughout the script as a simple way to track progress. Once again we'll break this script up with explanations, for clarity.

现在我们有了数据，我们可以开始我们的洪水填充功能。

注水功能

我们在这个例子中使用了 ASCII 网格，这意味着这个模型的引擎完全在 NumPy 中。我们首先定义floodFill()函数，这是这个模型的核心和灵魂。这篇维基百科关于洪水填充算法的文章提供了不同方法的极好概述:en.wikipedia.org/wiki/Flood_…。

泛洪填充算法从给定的单元格开始，并开始检查相邻单元格的相似性。相似性因素可能是颜色，或者在我们的例子中是海拔。如果相邻像元与当前像元的高程相同或更低，则该像元将被标记为检查其相邻像元，直到检查完整个网格。NumPy 并不是为了以这种方式在数组上爬行而设计的，但是它在处理多维数组时仍然是有效的。我们穿过每个牢房，检查它的北面、南面、东面和西面的邻居。那些可以被淹没的单元中的任何一个都被添加到填充集合中，并且它们的邻居被添加到填充集合中以被算法检查。

如前所述，如果您试图将相同的值添加到集合中两次，它只会忽略重复的条目并维护一个唯一的列表。通过在数组中使用集合，我们可以有效地只检查一次单元格，因为填充集合包含唯一的单元格。下面的代码实现了我们的floodFill功能:

首先，我们导入我们的库:

import numpy as np
from linecache import getline

接下来，我们创建我们的floodFill函数:

def floodFill(c, r, mask):
 """
 Crawls a mask array containing
 only 1 and 0 values from the
 starting point (c=column,
 r=row - a.k.a. x, y) and returns
 an array with all 1 values
 connected to the starting cell.
 This algorithm performs a 4-way
 check non-recursively.
 """

接下来，我们创建集合来跟踪已经覆盖的单元格:

 # cells already filled
 filled = set()
 # cells to fill
 fill = set()
 fill.add((c, r))
 width = mask.shape[1]-1
 height = mask.shape[0]-1

然后我们创建我们的淹没数组:

 # Our output inundation array
 flood = np.zeros_like(mask, dtype=np.int8)

现在我们可以循环通过细胞并淹没它们，或者不这样做:

 # Loop through and modify the cells which
 # need to be checked.
 while fill:
   # Grab a cell
   x, y = fill.pop()

如果土地高于洪水，跳过它:

   if y == height or x == width or x < 0 or y < 0:
    # Don't fill
    continue

如果地面高程等于或小于洪水，填写:

   if mask[y][x] == 1:
    # Do fill
    flood[y][x] = 1
   filled.add((x, y))

现在，我们检查周围的相邻细胞，看看它们是否需要填充，当细胞用完时，我们返回淹没的矩阵:

   # Check neighbors for 1 values
   west = (x-1, y)
   east = (x+1, y)
   north = (x, y-1)
   south = (x, y+1)
   if west not in filled:
     fill.add(west)
   if east not in filled:
     fill.add(east)
   if north not in filled:
     fill.add(north)
   if south not in filled:
     fill.add(south)
 return flood

现在我们已经设置了floodFill功能，我们可以创建一个洪水。

洪水淹没预测

在脚本的剩余部分，我们从一个 ASCII 网格加载我们的地形数据，定义我们的输出网格文件名，并对地形数据执行算法。洪水填充算法的种子是任意点，如低海拔区域内的sx和sy。在现实世界的应用中，这些点可能是一个已知的位置，比如一个测流计或者大坝的裂缝。在最后一步，我们保存输出网格。

需要执行以下步骤:

首先，我们设置我们的source和target数据名称:

source = "terrain.asc"
target = "flood.asc"

接下来，我们打开源代码:

print("Opening image...")
img = np.loadtxt(source, skiprows=6)
print("Image opened")

我们将为70米以下的所有东西创建一个掩码数组:

# Mask elevations lower than 70 meters.
wet = np.where(img < 70, 1, 0)
print("Image masked")

现在，我们将从头部解析地理空间信息:

# Parse the header using a loop and
# the built-in linecache module
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(h.split(" ")[-1].strip()) for h in hdr]
cols, rows, lx, ly, cell, nd = values
xres = cell
yres = cell * -1

现在，我们将建立一个位于河床上的起点:

# Starting point for the
# flood inundation in pixel coordinates
sx = 2582
sy = 2057

现在，我们触发我们的floodFill功能:

print("Beginning flood fill")
fld = floodFill(sx, sy, wet)
print("Finished flood fill")

header = ""
for i in range(6):
 header += hdr[i]

最后，我们可以保存洪水淹没模型输出:

print("Saving grid")
# Open the output file, add the hdr, save the array
with open(target, "wb") as f:
 f.write(bytes(header, 'UTF-8'))
 np.savetxt(f, fld, fmt="%1i")
print("Done!")

下图显示了分类版数字高程模型的洪水淹没输出，其中较低的高程值为棕色，中等范围的值为绿色，较高的值为灰色和白色:

洪水栅格包括所有小于 70 米的区域，颜色为蓝色。此图像是用 QGIS 创建的，但它可以在 ArcGIS 中显示为 EPSG:4326。您也可以使用 GDAL 将洪水栅格保存为 8 位 TIFF 文件或 JPEG 文件，就像 NDVI 示例一样，以便在标准图形程序中查看它。

下面截图中的图像几乎完全相同，除了被过滤的蒙版，它显示为黄色。这是通过为名为wet的数组而不是fld生成一个文件来完成的，以显示不连续的区域，这些区域不包括在洪水中。这些区域没有连接到源点，因此在洪水期间不太可能到达:

通过更改高程值，可以创建额外的洪水淹没栅格。我们从海拔 70 米开始。如果我们将该值增加到 90，我们可以扩大洪水。以下截图显示了 70 米和 90 米处的洪水事件:

90 米的洪水是浅蓝色的多边形。你可以采取更大或更小的步骤，并以不同的层显示不同的影响。

这个模型是一个优秀的和有用的可视化。但是，您可以通过在洪水遮罩上使用 GDAL 的polygonize()方法来进一步进行分析，就像我们在第 6 章、 Python 和遥感中的从图像中提取特征部分中对岛屿所做的那样。这个操作会给你一个向量淹没多边形。然后，您可以使用我们在第 5 章、 Python 和地理信息系统中执行选择部分讨论的原则，使用多边形选择建筑物以确定人口影响。您也可以将洪水多边形与第 5 章、 Python 和地理信息系统的点密度计算部分中的点密度示例结合起来，以评估洪水对人口的潜在影响。可能性是无穷的。

创建彩色山体阴影

在本例中，我们将结合以前的技术，将来自第 7 章、 Python 和的地形山体阴影与我们在 LIDAR 上使用的颜色分类相结合。对于这个例子，我们将需要命名为dem.asc和relief.asc的 ASCII 网格 DEMs，我们在上一章中使用了它们。

我们将创建一个彩色的数字高程模型和一个山体阴影，然后使用 PIL 将它们混合在一起，以增强高程可视化。代码注释将指导您完成示例，因为您已经熟悉其中的许多步骤:

首先，我们导入我们需要的库:

import gdal_array as gd
try:
 import Image
except ImportError:
 from PIL import Image

对于下一部分，您将需要以下两个文件:https://github . com/GeospatialPython/Learn/raw/master/relief . zip和https://github . com/GeospatialPython/Learn/raw/master/DEM . zip。

然后，我们将为输入和输出设置变量:

relief = "relief.asc"
dem = "dem.asc"
target = "hillshade.tif"

接下来，我们将加载我们的relief图像:

# Load the relief as the background image
bg = gd.numpy.loadtxt(relief, skiprows=6)

然后，我们将加载数字高程模型图像，这样我们将获得高程数据:

# Load the DEM into a numpy array as the foreground image
fg = gd.numpy.loadtxt(dem, skiprows=6)[:-2, :-2]

现在，我们将为我们的彩色化创建一个新图像，其中高程断点在 LUT 中形成类和相应的颜色:

# Create a blank 3-band image to colorize the DEM
rgb = gd.numpy.zeros((3, len(fg), len(fg[0])), gd.numpy.uint8)

# Class list with DEM upper elevation range values.
classes = [356, 649, 942, 1235, 1528,
 1821, 2114, 2300, 2700]

# Color look-up table (lut)
# The lut must match the number of classes.
# Specified as R, G, B tuples
lut = [[63, 159, 152], [96, 235, 155], [100, 246, 174],
 [248, 251, 155], [246, 190, 39], [242, 155, 39],
 [165, 84, 26], [236, 119, 83], [203, 203, 203]]

# Starting elevation value of the first class
start = 1

我们现在可以执行我们的颜色分类:

# Process all classes.
for i in range(len(classes)):
 mask = gd.numpy.logical_and(start <= fg,
 fg <= classes[i])
 for j in range(len(lut[i])):
 rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j]))
 start = classes[i]+1

然后，我们可以将我们的着色浮雕数组转换为图像，以及我们的彩色数字高程模型:

# Convert the shaded relief to a PIL image
im1 = Image.fromarray(bg).convert('RGB')

# Convert the colorized DEM to a PIL image.
# We must transpose it from the Numpy row, col order
# to the PIL col, row order (width, height).
im2 = Image.fromarray(rgb.transpose(1, 2, 0)).convert('RGB')

现在，我们将混合两个图像以获得最终效果，并将其保存到图像文件中:

# Blend the two images with a 40% alpha
hillshade = Image.blend(im1, im2, .4)

# Save the hillshade
hillshade.save(target)

下图显示了输出，这为地理信息系统地图提供了很好的背景:

既然我们可以模拟地形，让我们学习如何在上面导航。

执行最低成本路径分析

计算行车方向是世界上最常用的地理空间功能。通常，这些算法计算点 A 和 B 之间的最短路径，或者它们可以考虑道路的速度限制，甚至当前的交通状况，以便通过行驶时间来选择路线。

但是如果你的工作是修建一条新路呢？或者，如果你负责决定在偏远地区的哪里铺设输电线路或输水线路呢？在基于地形的环境中，最短的路径可能会穿过一座困难的山，或者穿过一个湖。在这种情况下，我们需要考虑障碍，并尽可能避免它们。然而，如果避开一个小障碍让我们走得太远，实施这条路线的成本可能比仅仅翻越一座山还要高。

这种类型的高级分析称为最小成本路径分析。我们在一个区域内寻找一条路线，这条路线是距离与遵循这条路线的成本之间的最佳折衷。我们用于此过程的算法称为 A 星或 A* 算法。最古老的路由方法被称为迪克斯特拉算法，它计算网络中的最短路径，例如道路网络。A*方法也可以做到这一点，但它也更适合遍历网格状的数字高程模型。

You can find out more about these algorithms on the following web pages:

迪克斯特拉算法:en.wikipedia.org/wiki/Dijkst…。
A*算法:en.wikipedia.org/wiki/A-star…。

这个例子是本章中最复杂的。为了更好地理解它，我们有一个简单版本的程序，它是基于文本的，在一个 5×5 的网格上运行，随机生成值。实际上，您可以看到这个程序是如何遵循算法的，然后在具有数千个值的高程网格上进行尝试。

该程序执行以下步骤:

创建一个简单的网格，其中随机生成的伪高程值介于 1 和 16 之间。
在网格的左下角定义一个开始位置。
将端点定义为网格的右上角。
创建一个成本网格，其中包含每个单元的高程，以及单元到终点的距离。
从一开始就检查每个相邻的小区，选择成本最低的一个。
使用选定的单元格重复评估，直到我们到达最后。
返回所选单元格集作为最小成本路径。
设置测试网格。

您只需从命令行运行该程序并查看其输出。该脚本的第一部分将我们的人工地形网格设置为随机生成的 NumPy 数组，名义高程值在 1 到 16 之间。我们还创建了一个距离网格，用于计算每个像元到目标像元的距离。该值是每个单元的成本。

让我们看看以下步骤:

首先，我们将导入numpy并设置网格的大小:

import numpy as np

# Width and height
# of grids
w = 5
h = 5

接下来，我们设置开始位置单元格和结束位置:

# Start location:
# Lower left of grid
start = (h-1, 0)

# End location:
# Top right of grid
dx = w-1
dy = 0

现在，我们可以根据宽度和高度创建一个零网格:

# Blank grid
blank = np.zeros((w, h))

接下来，我们将设置距离网格，以创建阻抗值:

# Distance grid
dist = np.zeros(blank.shape, dtype=np.int8)

# Calculate distance for all cells
for y, x in np.ndindex(blank.shape):
 dist[y][x] = abs((dx-x)+(dy-y))

现在，我们将打印出成本网格中每个单元格的成本值:

# "Terrain" is a random value between 1-16.
# Add to the distance grid to calculate
# The cost of moving to a cell
cost = np.random.randint(1, 16, (w, h)) + dist

print("COST GRID (Value + Distance)\n{}\n".format(cost))

现在我们有一个模拟地形网格可以使用，我们可以测试一个路由算法。

简单的 A*算法

这里实现的 A*搜索算法以类似于上一个示例中的洪水填充算法的方式爬行网格。再一次，我们使用集合来避免使用递归，并避免单元检查的重复。但是这一次，我们没有检查海拔，而是检查通过有问题的小区的路由的距离成本。如果这一举动增加了到达终点的成本，那么我们会选择成本更低的方案。

需要执行以下步骤，如下所示:

首先，我们将通过创建跟踪路径进度的集合来启动我们的 A*函数:

# Our A* search algorithm
def astar(start, end, h, g):
    closed_set = set()
    open_set = set()
    path = set()

接下来，我们将起始单元格添加到打开的单元格列表中，以便处理并开始循环遍历该集合:

    open_set.add(start)
    while open_set:
        cur = open_set.pop()
        if cur == end:
            return path
        closed_set.add(cur)
        path.add(cur)
        options = []
        y1 = cur[0]
        x1 = cur[1]

我们检查周围的单元格作为前进的选项:

        if y1 > 0:
            options.append((y1-1, x1))
        if y1 < h.shape[0]-1:
            options.append((y1+1, x1))
        if x1 > 0:
            options.append((y1, x1-1))
        if x1 < h.shape[1]-1:
            options.append((y1, x1+1))
        if end in options:
            return path
        best = options[0]
        closed_set.add(options[0])

然后，我们检查每个选项，找出最佳选项，并将其附加到路径中，直到到达终点:

        for i in range(1, len(options)):
            option = options[i]
            if option in closed_set:
                continue
            elif h[option] <= h[best]:
                best = option
                closed_set.add(option)
            elif g[option] < g[best]:
                best = option
                closed_set.add(option)
            else:
                closed_set.add(option)
        print(best, ", ", h[best], ", ", g[best])
        open_set.add(best)
    return []

现在我们已经建立了算法，我们可以通过创建一个路径来测试它。

生成测试路径

在本节中，我们将在测试网格上生成一个路径。我们将使用起点、终点、成本网格和距离网格来调用我们的 A*函数:

# Find the path
path = astar(start, (dy, dx), cost, dist)
print()

现在，我们将把我们的路径放在它自己的网格上并打印出来:

# Create and populate the path grid
path_grid = np.zeros(cost.shape, dtype=np.uint8)
for y, x in path:
 path_grid[y][x] = 1
path_grid[dy][dx] = 1

print("PATH GRID: 1=path")
print(path_grid)

接下来，我们将查看这个测试的输出。

查看测试输出

运行该程序时，您将生成一个类似于以下内容的随机编号网格:

COST GRID (Value + Distance)
[[13 10 5 15 9]
 [15 13 16 5 16]
 [17 8 9 9 17]
 [ 4 1 11 6 12]
 [ 2 7 7 11 8]]

(Y,X), HEURISTIC, DISTANCE
(3, 0) , 4 , 1
(3, 1) , 1 , 0
(2, 1) , 8 , 1
(2, 2) , 9 , 0
(2, 3) , 9 , 1
(1, 3) , 5 , 0
(0, 3) , 15 , 1

PATH GRID: 1=path
[[0 0 0 1 1]
 [0 0 0 1 0]
 [0 1 1 1 0]
 [1 1 0 0 0]
 [1 0 0 0 0]]

网格足够小，以便您可以轻松地手动跟踪算法的步骤。该实现使用曼哈顿距离，这意味着该距离不使用对角线，仅使用左、右、上、下测量值。为了保持简单，搜索也不会沿对角线移动。

现实世界的例子

现在我们已经对 A算法有了基本的了解，让我们来看一个更复杂的例子。对于地形示例，我们将使用位于加拿大不列颠哥伦比亚省温哥华附近的相同 DEM，我们在创建阴影地形*部分的第 7 章、 Python 和高程数据中使用了该 DEM。该网格的空间参考是 EPSG:26910 NAD 83/UTM 10N 区。您可以从git.io/v3fpL下载形状文件的数字高程模型、地形起伏以及起点和终点作为压缩包。

我们将实际使用阴影浮雕进行可视化。我们在本练习中的目标是以尽可能低的成本从起点到终点:

单看地形，有两条路径走低海拔路线，方向没有太大变化。下面的截图说明了这两条路线:

所以，我们期望当我们使用 A*算法时，它会很接近。请记住，该算法只查看紧邻区域，因此它不能像我们一样查看整个图像，也不能根据前面已知的障碍物在路线的早期进行调整。

我们将从我们的简单示例中扩展这个实现，并使用欧几里德距离，或乌鸦飞翔时的测量，我们还将允许搜索向八个方向而不是四个方向看。我们将优先考虑地形作为主要决策点。为了确保我们朝着目标前进，并且不会偏离轨道太远，我们还会把从终点到起点的距离作为较低的优先级。除了这些区别之外，步骤与简单的示例相同。输出将是一个栅格，路径值设置为 1，其他值设置为零。

*既然理解了问题，那就来解决吧！

正在加载网格

在这一节和接下来的几节中，我们将创建可以在地形上创建路线的脚本。剧本一开始很简单。我们将网格从 ASCII 网格加载到 NumPy 数组中。我们命名输出路径网格，然后定义开始单元格和结束单元格:

首先，我们导入我们的库:

import numpy as np
import math
from linecache import getline
import pickle

接下来，我们将定义输入和输出数据源:

# Our terrain data
source = "dem.asc"

# Output file name for the path raster
target = "path.asc"

然后，我们可以跳过标题加载网格:

print("Opening %s..." % source)
cost = np.loadtxt(source, skiprows=6)
print("Opened %s." % source)

接下来，我们将解析地理空间和网格大小信息的标题:

# Parse the header
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(ln.split(" ")[-1].strip()) for ln in hdr]
cols, rows, lx, ly, cell, nd = values

最后，我们将定义我们的开始和结束位置:

# Starting column, row
sx = 1006
sy = 954

# Ending column, row
dx = 303
dy = 109

现在我们的网格已经加载，我们可以设置我们需要的功能。

定义助手函数

我们需要三个功能来在地形上走路线。一个是 A*算法，另外两个辅助算法选择下一步。我们将简要讨论这些助手函数。首先，我们有一个简单的欧几里德距离函数e_dist，它返回两点之间的直线距离作为地图单位。接下来，我们有一个名为weighted_score的重要函数，它根据邻居和当前小区之间的海拔变化以及到目的地的距离，返回相邻小区的分数。

这个功能比单独的距离或高度更好，因为它减少了两个单元之间存在联系的机会，从而更容易避免回溯。这个评分公式大致基于一个名为尼森评分的概念，该概念通常用于这些类型的算法中，并在本章前面提到的维基百科文章中引用。这个函数的伟大之处在于，它可以用任何你想要的值给相邻的单元格打分。您还可以使用实时提要来查看相邻单元格中的当前天气，并避开有雨或雪的单元格。

下面的代码将创建我们穿越地形所需的距离函数和权重函数:

首先，我们将创建一个欧几里德距离函数，它将给出点之间的距离:

def e_dist(p1, p2):
 """
 Takes two points and returns
 the Euclidian distance
 """
 x1, y1 = p1
 x2, y2 = p2
 distance = math.sqrt((x1-x2)**2+(y1-y2)**2)
 return int(distance)

现在，我们将创建权重函数，以便为每个节点的移动适宜性评分:

def weighted_score(cur, node, h, start, end):
 """
 Provides a weighted score by comparing the
 current node with a neighboring node. Loosely
 based on the Nisson Score concept: f=g+h
 In this case, the "h" value, or "heuristic",
 is the elevation value of each node.
 """

我们从0的score开始，检查节点距离终点和起点的距离:

 score = 0
 # current node elevation
 cur_h = h[cur]
 # current node distance from end
 cur_g = e_dist(cur, end)
 # current node distance from
 cur_d = e_dist(cur, start)

接下来，我们检查相邻节点，并决定移动到哪里:

 # neighbor node elevation
 node_h = h[node]
 # neighbor node distance from end
 node_g = e_dist(node, end)
 # neighbor node distance from start
 node_d = e_dist(node, start)
 # Compare values with the highest
 # weight given to terrain followed
 # by progress towards the goal.
 if node_h < cur_h:
 score += cur_h-node_h
 if node_g < cur_g:
 score += 10
 if node_d > cur_d:
 score += 10
 return score

现在我们的助手函数已经完成，我们可以构建 A*函数了。

真实世界的 A*算法

这个算法比我们前面例子中的简单版本更复杂。我们使用集合来避免冗余。它还实现了我们更先进的评分算法，并在进行额外计算之前进行检查，以确保我们没有到达路径的末端。与我们的上一个示例不同，这个更高级的版本还检查八个方向的单元格，因此路径可以对角移动。这个函数的末尾有一个print语句被注释掉了。您可以取消对它的注释，以便观察搜索在网格中的爬行。下面的代码将实现我们将在本节剩余部分使用的 A*算法:

首先，我们通过接受一个起点、一个终点和一个分数来打开函数:

def astar(start, end, h):
 """
 A-Star (or A*) search algorithm.
 Moves through nodes in a network
 (or grid), scores each node's
 neighbors, and goes to the node
 with the best score until it finds
 the end. A* is an evolved Dijkstra
 algorithm.
 """

现在，我们设置了跟踪进度的设置:

 # Closed set of nodes to avoid
 closed_set = set()
 # Open set of nodes to evaluate
 open_set = set()
 # Output set of path nodes
 path = set()

接下来，我们从起点开始处理:

 # Add the starting point to
 # to begin processing
 open_set.add(start)
 while open_set:
 # Grab the next node
 cur = open_set.pop()

如果我们到达终点，我们返回完整的路径:

 # Return if we're at the end
 if cur == end:
 return path

否则，我们将继续通过网格工作并消除各种可能性:

 # Close off this node to future
 # processing
 closed_set.add(cur)
 # The current node is always
 # a path node by definition
 path.add(cur)

为了保持移动，我们会在移动时抓取所有需要处理的邻居:

 # List to hold neighboring
 # nodes for processing
 options = []
 # Grab all of the neighbors
 y1 = cur[0]
 x1 = cur[1]
 if y1 > 0:
 options.append((y1-1, x1))
 if y1 < h.shape[0]-1:
 options.append((y1+1, x1))
 if x1 > 0:
 options.append((y1, x1-1))
 if x1 < h.shape[1]-1:
 options.append((y1, x1+1))
 if x1 > 0 and y1 > 0:
 options.append((y1-1, x1-1))
 if y1 < h.shape[0]-1 and x1 < h.shape[1]-1:
 options.append((y1+1, x1+1))
 if y1 < h.shape[0]-1 and x1 > 0:
 options.append((y1+1, x1-1))
 if y1 > 0 and x1 < h.shape[1]-1:
 options.append((y1-1, x1+1))

我们检查每个邻居是否是目的地:

 # If the end is a neighbor, return
 if end in options:
 return path

我们将第一个选项作为best选项，并处理其他选项，边走边升级:

 # Store the best known node
 best = options[0]
 # Begin scoring neighbors
 best_score = weighted_score(cur, best, h, start, end)
 # process the other 7 neighbors
 for i in range(1, len(options)):
 option = options[i]
 # Make sure the node is new
 if option in closed_set:
 continue
 else:
 # Score the option and compare 
 # it to the best known
 option_score = weighted_score(cur, option, 
 h, start, end)
 if option_score > best_score:
 best = option
 best_score = option_score
 else:
 # If the node isn't better seal it off
 closed_set.add(option)
 # Uncomment this print statement to watch
 # the path develop in real time:
 # print(best, e_dist(best, end))
 # Add the best node to the open set
 open_set.add(best)
return []

既然我们有了路由算法，我们就可以生成真实世界的路径。

生成真实世界的路径

最后，我们将现实世界的路径创建为零网格中的一串 1。然后，可以将该栅格引入到 QGIS 等应用程序中，并在地形网格上进行可视化。在下面的代码中，我们将使用我们的算法和助手函数来生成路径，如下所示:

首先，我们将起点和终点以及地形网格发送给路由功能:

print("Searching for path...")
p = astar((sy, sx), (dy, dx), cost)
print("Path found.")
print("Creating path grid...")
path = np.zeros(cost.shape)
print("Plotting path...")
for y, x in p:
 path[y][x] = 1
path[dy][dx] = 1
print("Path plotted.")

一旦我们有了路径，我们就可以把它保存为 ASCII 网格:

print("Saving %s..." % target)
header = ""
for i in range(6):
 header += hdr[i]

# Open the output file, add the hdr, save the array
with open(target, "wb") as f:
 f.write(bytes(header, 'UTF-8'))
 np.savetxt(f, path, fmt="%4i")

现在，我们希望保存路径数据，因为从起点到终点，点的顺序是正确的。当我们将它们放入栅格中时，我们会失去这个顺序，因为它都是一个栅格。我们将使用内置的 Python pickle模块将列表对象保存到磁盘。我们将在下一节中使用这些数据来创建路线的向量形状文件。因此，我们将路径数据保存为一个腌制的 Python 对象，以后可以重用，而无需运行整个程序:

print("Saving path data...")
with open("path.p", "wb") as pathFile:
 pickle.dump(p, pathFile)
print("Done!")

以下是我们搜索的输出路径:

如您所见，A*搜索非常接近我们手动选择的路线之一。在一些情况下，算法选择处理一些地形，而不是试图绕过它。有时，轻微的地形被认为比绕过它的距离花费更少。您可以在路线右上角的放大部分看到这种选择的示例。红线是我们的程序通过地形生成的路线:

我们只使用了两个值:地形和距离。但是你也可以添加数百个因素，如土壤类型、水体和现有道路。所有这些项目都可以作为阻抗或直接墙。您只需修改示例中的评分函数，以考虑任何附加因素。请记住，您添加的因素越多，就越难追踪 A实现在选择路线时在想什么*。

这种分析的一个明显的未来方向是将这条路线创建为一条线的向量版本。该过程将包括将每个单元映射到一个点，然后在将其保存为 shapefile 或 GeoJSON 文件之前，使用最近邻分析对这些点进行正确排序。

将路线转换为形状文件

最低成本路径路线的栅格版本对于可视化很有用，但它对于分析来说并不太好，因为它嵌入在栅格中，因此很难与其他数据集相关联，就像我们在本书中多次做的那样。我们的下一个目标是使用我们在创建路径时保存的路径数据来创建一个 shapefile，因为保存的数据是以正确的顺序保存的。下面的代码将把我们的栅格路径转换成一个更容易在地理信息系统中用于分析的形状文件:

首先，我们将导入我们需要的模块，这些模块并不多。我们将使用pickle模块来恢复路径data对象。然后，我们将使用linecache模块从路径栅格中读取地理空间标题信息，以便将路径行和列映射到地球坐标。最后，我们将使用shapefile模块导出 shapefile:

import pickle
from linecache import getline
import shapefile

接下来，我们将创建一个函数，将行和列转换为 x 和 y 坐标。该函数接受来自路径栅格文件的元数据头信息，以及列号和行号:

def pix2coord(gt,x,y):
 geotransform = gt
 ox = gt[2]
 oy = gt[3]
 pw = gt[4]
 ph = gt[4]
 cx = ox + pw * x + (pw/2)
 cy = oy + pw * y + (ph/2)
 return cx, cy

现在，我们将从腌制对象中恢复path对象:

with open("path.p", "rb") as pathFile:
 path = pickle.load(pathFile)

然后，我们将解析路径栅格文件中的元数据信息:

hdr = [getline("path.asc", i) for i in range(1, 7)]
gt = [float(ln.split(" ")[-1].strip()) for ln in hdr]

接下来，我们需要一个列表对象来保存转换后的坐标:

coords = []

现在，我们将每个栅格位置从最低成本路径对象转换为地理空间坐标，并将其存储在我们创建的列表中:

for y,x in path:
 coords.append(pix2coord(gt,x,y))

最后，只需几行，我们就可以写出一个线条形状文件:

with shapefile.Writer("path", shapeType=shapefile.POLYLINE) as w:
 w.field("NAME")
 w.record("LeastCostPath")
 w.line([coords])

干得好！您已经创建了一个程序，该程序可以根据一组规则自动通过障碍物，并将其导出到一个文件中，您可以在地理信息系统中显示和分析该文件！我们只使用了三个规则，但是您可以通过添加其他数据集，例如天气或水体，或者您可以想象的任何其他东西，来添加对程序如何选择路径的附加限制。

既然我们理解了在任意表面上开辟一条路径，我们就来看看通过网络的路由。

计算卫星图像云量

卫星图像给我们提供了一个强有力的地球鸟瞰图。它们有多种用途，我们在第 6 章、Python 和遥感中看到过。然而，它们有一个缺点——云。当卫星绕地球运行并收集图像时，它不可避免地会对云成像。除了阻碍我们对地球的观察，云数据还会通过在无用的云数据上浪费 CPU 周期对遥感算法产生不利影响，或者通过引入不需要的数据值来扭曲结果。

解决方案是创建云遮罩。云遮罩是将云数据隔离在单独栅格中的栅格。然后，您可以在处理图像时使用该栅格作为参考，以避免云数据，或者您甚至可以使用它从原始图像中移除云。

在本节中，我们将使用rasterio模块和rio-l8qa插件为 Landsat 图像创建一个云遮罩。云遮罩将被创建为单独的图像，仅包含云:

首先，我们需要从bit.ly/landsat8dat…下载一些样本 Landsat 8 卫星图像数据作为 ZIP 文件。
点击右上角的下载图标，将数据下载为 ZIP 文件，解压到名为l8的目录下。
接下来，通过运行pip，确保您拥有我们需要的栅格库:

pip install rasterio
pip install rio-l8qa

现在，我们将通过首先导入所需的库来创建云遮罩:

import glob
import os
import rasterio
from l8qa.qa import write_cloud_mask

接下来，我们需要提供我们卫星图像目录的参考:

# Directory containing landsat data
landsat_dir = "l8"

现在，我们需要找到卫星数据的质量保证元数据，它为我们提供了生成云遮罩所需的信息:

src_qa = glob.glob(os.path.join(landsat_dir, '*QA*'))[0]

最后，我们使用质量保证文件来创建云遮罩 TIFF 文件:

with rasterio.open(src_qa) as qa_raster:
 profile = qa_raster.profile
 profile.update(nodata=0)
 write_cloud_mask(qa_raster.read(1), profile, 'cloudmask.tif')

下面的图像只是来自 Landsat 8 数据集的波段 7(短波红外)图像:

下一个图像是仅包含云和阴影位置的云遮罩图像:

最后，这是图像上的面具，显示云是黑色的:

这个例子刷表面，你可以用图像蒙版。另一个rasterio模块rio-cloudmask允许您从头开始计算云遮罩，而无需使用质量保证数据。但是它需要一些额外的预处理步骤。您可以在这里了解更多信息:github.com/mapbox/rio-…

沿着街道走

沿着街道走的路线使用一个由线组成的连接网络，这被称为图形。图中的线可以有阻抗值，这阻碍了路由算法将它们包括在路由中。阻抗值的例子通常包括交通量、速度限制甚至距离。布线图的一个关键要求是所有的线(称为边)都必须连接起来。为映射而创建的道路数据集通常具有节点不相交的线。

在这个例子中，我们将根据距离计算通过图表的最短路径。我们将使用起点和终点，它们不是图中的节点，这意味着我们必须首先找到离我们的起点和终点最近的图节点。

为了计算最短路线，我们将使用一个强大的纯 Python 图形库，称为网络。NetworkX 是一个通用的网络图形库，可以创建、操作和分析复杂的网络，包括地理空间网络。如果pip没有在你的系统上安装 NetworkX，那么你可以在networkx.readthedocs.org/en/stable/找到不同操作系统的 NetworkX 下载安装说明。

您可以从git.io/vcXFQ下载位于美国墨西哥湾沿岸的公路网以及起点和终点的 ZIP 文件。然后，您可以按照以下步骤操作:

首先，我们需要导入将要使用的库。除了网络之外，我们还将使用 PyShp 库来读写形状文件:

import networkx as nx
import math
from itertools import tee
import shapefile
import os

接下来，我们将当前目录定义为我们将要创建的 route shapefile 的输出目录:

savedir = "."

现在，我们需要一个能够计算点与点之间距离的函数，以便填充图形的阻抗值，并找到距离路线起点和终点最近的节点:

def haversine(n0, n1):
 x1, y1 = n0
 x2, y2 = n1
 x_dist = math.radians(x1 - x2)
 y_dist = math.radians(y1 - y2)
 y1_rad = math.radians(y1)
 y2_rad = math.radians(y2)
 a = math.sin(y_dist/2)**2 + math.sin(x_dist/2)**2 \
 * math.cos(y1_rad) * math.cos(y2_rad)
 c = 2 * math.asin(math.sqrt(a))
 distance = c * 6371
 return distance

然后，我们将创建另一个函数，该函数从列表中返回成对的点，以给出我们将用来构建图形边的线段:

def pairwise(iterable):
 """Return an iterable in tuples of two
 s -> (s0,s1), (s1,s2), (s2, s3), ..."""
 a, b = tee(iterable)
 next(b, None)
 return zip(a, b)

现在，我们将定义我们的道路网络形状文件。该道路网是美国地质调查局的美国州际公路文件 shapefile】美国地质勘探局)的子集，该文件已经过编辑，以确保所有道路都已连接:

shp = "road_network.shp"

接下来，我们将使用 NetworkX 创建一个图形，并将 shapefile 段添加为图形边:

G = nx.DiGraph()
r = shapefile.Reader(shp)
for s in r.shapes():
 for p1, p2 in pairwise(s.points):
 G.add_edge(tuple(p1), tuple(p2))

然后，我们可以提取连接的组件作为子图。然而，在这种情况下，我们已经确保了整个图是连通的:

sg = list(nx.connected_component_subgraphs(G.to_undirected()))[0]

接下来，我们可以阅读我们想要导航的start和end点:

r = shapefile.Reader("start_end")
start = r.shape(0).points[0]
end = r.shape(1).points[0]

现在，我们循环通过图表，并使用我们的haversine公式为每条边分配距离值:

for n0, n1 in sg.edges_iter():
 dist = haversine(n0, n1)
 sg.edge[n0][n1]["dist"] = dist

接下来，我们必须在图中找到离我们的起点和终点最近的节点，以便通过循环遍历所有节点来开始和结束我们的路线，并测量到我们的终点的距离，直到我们找到最短的距离:

nn_start = None
nn_end = None
start_delta = float("inf")
end_delta = float("inf")
for n in sg.nodes():
 s_dist = haversine(start, n)
 e_dist = haversine(end, n)
 if s_dist < start_delta:
 nn_start = n
 start_delta = s_dist
 if e_dist < end_delta:
 nn_end = n 
 end_delta = e_dist

现在，我们准备计算通过我们的道路网的最短距离:

path = nx.shortest_path(sg, source=nn_start, target=nn_end, weight="dist")

最后，我们将结果添加到 shapefile 并保存我们的路线:

w = shapefile.Writer(shapefile.POLYLINE)
w.field("NAME", "C", 40)
w.line(parts=[[list(p) for p in path]])
w.record("route")
w.save(os.path.join(savedir, "route"))

下面的截图显示了浅灰色的道路网，起点和终点，以及黑色的路线。您可以看到，路线穿过道路网络，以便在最短的距离内到达距离终点最近的道路:

现在我们知道了如何创建各种类型的路线，我们可以看看如何定位您在路线上旅行时可能拍摄的照片。

地理定位照片

使用支持全球定位系统的相机拍摄的照片，包括智能手机，以一种称为 EXIF 标签的格式将位置信息存储在文件的标题中。这些标签主要基于 TIFF 图像标准使用的相同标题标签。在本例中，我们将使用这些标签创建一个形状文件，其中包含照片的点位置和照片的文件路径作为属性。

我们将在这个例子中使用 PIL，因为它能够提取 EXIF 数据。大多数用智能手机拍摄的照片都是带有地理标签的图像；但是，您可以从git.io/vczR0下载本例中使用的器械包:

首先，我们将导入我们需要的库，包括用于图像元数据的 PIL 和用于形状文件的 PyShp:

import glob
import os
try:
 import Image
 import ImageDraw
except ImportError:
 from PIL import Image
 from PIL.ExifTags import TAGS
import shapefile

现在，我们需要三个函数。第一个函数提取 EXIF 数据。第二个函数将度、分、秒 ( DMS )坐标转换为十进制度(EXIF 数据将 GPS 数据存储为 DMS 坐标)。第三个函数提取全球定位系统数据并执行坐标转换:

def exif(img):
 # extract exif data.
 exif_data = {}
 try: 
 i = Image.open(img)
 tags = i._getexif()
 for tag, value in tags.items():
 decoded = TAGS.get(tag, tag)
 exif_data[decoded] = value
 except:
 pass
 return exif_data

def dms2dd(d, m, s, i):
 # convert degrees, min, sec to decimal degrees
 sec = float((m * 60) + s)
 dec = float(sec / 3600)
 deg = float(d + dec)
 if i.upper() == 'W':
 deg = deg * -1
 elif i.upper() == 'S':
 deg = deg * -1
 return float(deg)

def gps(exif):
 # get gps data from exif
 lat = None
 lon = None
 if exif['GPSInfo']: 
 # Lat
 coords = exif['GPSInfo']
 i = coords[1]
 d = coords[2][0][0]
 m = coords[2][1][0]
 s = coords[2][2][0]
 lat = dms2dd(d, m, s, i)
 # Lon
 i = coords[3]
 d = coords[4][0][0]
 m = coords[4][1][0]
 s = coords[4][2][0]
 lon = dms2dd(d, m, s, i)
 return lat, lon

接下来，我们将循环浏览照片，提取坐标，并将坐标和文件名存储在字典中:

photos = {}
photo_dir = "./photos"
files = glob.glob(os.path.join(photo_dir, "*.jpg"))
for f in files:
 e = exif(f)
 lat, lon = gps(e)
 photos[f] = [lon, lat]

现在，我们将照片信息保存为形状文件:

with shapefile.Writer("photos", shapefile.POINT) as w:
    w.field("NAME", "C", 80)
    for f, coords in photos.items():
        w.point(*coords)
        w.record(f)

shapefile 中照片的文件名现在是照片拍摄点位置的属性。当您单击照片路径或点时，包括 QGIS 和 ArcGIS 在内的 GIS 程序具有将这些属性转换为链接的工具。以下来自 QGIS 的截图显示，其中一张照片在使用运行要素操作工具点击关联点后打开:

要查看结果，请使用以下说明:

从qgis.org下载 QGIS，按照安装说明进行。
打开 QGIS，将photos.shp文件拖到空白地图上。
在左侧的图层面板中，右键单击名为照片的图层，然后选择属性。
在操作选项卡上，单击绿色加号打开新的操作对话框。
在类型下拉菜单中，选择打开。
在描述字段中，输入打开图像。
单击右下角的“插入”按钮。
单击“确定”按钮，然后关闭属性对话框。
单击运行特征操作工具右侧的黑色小箭头，这是一个带有绿色中心和白色箭头的齿轮图标。
在弹出的菜单中，选择“打开图像”。
现在，点击地图上的一个点来查看地理标记的图像弹出窗口。

现在，让我们通过使用卫星图像，从地球上拍摄的图像转移到地球本身拍摄的图像。

摘要

在本章中，我们学习了如何创建三个真实世界的产品，它们每天都在政府、科学和工业中使用。除了通常使用黑盒包进行分析的地方——花费数千美元——我们能够使用非常少且免费的跨平台 Python 工具。除了本章中的例子之外，您现在还有一些更可重用的函数、算法和处理框架用于其他高级分析，这将允许您解决在交通、农业和天气等领域遇到的新问题。

在下一章中，我们将进入一个相对较新的地理空间分析领域:实时和接近实时的数据。*