Python 地理空间数据分析（二）
原文：
zh.annas-archive.org/md5/1e5dd3cb52782f4aa8fc65ed40d9831a译者：飞龙
第六章：ArcGIS Python API
由 Esri 开发的客户端软件和 GIS 组成的套件被称为 ArcGIS。它不仅是一个 API，还是一个 Python 包，允许用户查询存储在 ArcGIS Online 或 ArcGIS Enterprise 中的信息。它不是开源的，但 Esri 在行业中的领导地位已经产生了大量免费内容和教程，供你访问和探索。我将分享一些学习资源、可访问的工具和信息，你可以通过 ArcGIS Python API 访问它们。你需要使用 Jupyter Notebook 和 ArcGIS Online 来跟进。
设置
ArcGIS Pro 是一个桌面应用程序，还有一个基于浏览器的平台叫做 ArcGIS Online；在本章中，你将使用 ArcGIS Online 和 ArcGIS Python API 进行工作。
ArcGIS Python API 中可用的模块
为了扩展标准 Python 库，ArcGIS Python API 允许通过点符号访问其他模块。你可以在 文档 中探索它们，随着我们在编码中使用它们，我会提供更多信息。
Arcgis.gis
将你连接到 ArcGIS Online
Arcgis.features
提供用于处理地理元素组的功能
Arcgis.geometry
允许输入和输出不同的几何体（例如点或多边形）
Arcgis.geocoding
为地图可视化分配位置坐标
Arcgis.geoenrichment
为区域添加属性
在教授地理空间技能时，我使用开源或低成本选项，以便让尽可能多的人获得这些知识。由于大多数 Esri 用户都是专业或企业订阅者，很少有人了解个人使用的成本负担。无论你选择什么，请在使用之前仔细阅读文档，以避免不必要的费用。
ArcGIS API 是作为 arcgis 包分发的。在下面的示例中，你将通过 Conda 安装它，这是建议的协议，就像你在本书的几个章节中已经做过的那样。
安装 ArcGIS Pro
虽然 ArcGIS Pro 只适用于 Windows，但 API 允许非 Windows 用户访问 ArcGIS 功能，而无需打开 Windows 环境。我假设如果你已经有完整的 ArcGIS Pro 许可证，你可以按照文档安装 Python Package Manager or Python Command Prompt。要在桌面环境中使用 Pro 许可证并与 API 一起工作，你需要在运行命令的同一台计算机上安装 ArcGIS。
提示
如果你和我一样拥有 ArcGIS Pro 账户，但更喜欢在 MacOS 上工作，你可以离线使用 ArcGIS Pro。为此，请登录你的 ArcGIS Pro 账户，在“设置 >> 许可证”下，选中允许你离线工作的复选框。
即使没有许可证，也有访问公开资源的选项，因此让我们探索一下。
设置你的环境
接下来，打开您的终端。我建议为下载 ArcGIS 创建一个环境。（您可能记得，环境允许您安装兼容版本的软件包。）在这里，我称我的环境为 esriENV。
提示
虽然可以在环境中更新单个软件包，但如果在环境中工作时生成持久错误，我会将整个环境删除并从头开始创建。这是一个新的 Python 地理空间软件包，其中包含您可能希望考虑添加到任何创建的空间中的常用软件包：
mamba install -c conda-forge geospatial
您可能对这里列出的 Mamba 安装感兴趣。自本书编写以来，mambaforge 已成为 Conda-forge 社区项目。在操作上，Mamba 承诺速度和与 Conda 软件包的兼容性。在您的实验中，您将找到自己喜欢的工作流程。我已将 Mamba 安装到我的基础环境中，并在自己的项目中使用它，但本文是在长期可靠性的基础上编写的 Conda 环境。
让我们回顾一下如何创建您的环境并安装本笔记本所需的软件包。在终端中输入以下代码，但是在 *myenv* 的地方插入 您 的环境名称：
conda create --name myenv
要创建一个环境并包含一个特定版本的 Python 来解决软件包依赖关系，请输入以下代码：
#conda create --name myenv python=x.x
conda create --name esriENV python=3.9
现在，您必须激活您的新环境。目前我们只有一个频道或地方可以托管安装的软件包，即 conda-forge，因此我们希望所有需要安装软件包的依赖都来自 conda-forge 频道，除非它们只存在于默认频道上。
或者，您可能希望在列表中的任何频道上拥有最新版本的软件包。如果是这种情况，请使用conda config --set channel_priority strict：
conda activate esriENV
#conda config --set channel_priority strict
安装软件包
您已准备好安装您的软件包：
conda install -c esri arcgis
接下来，您将需要安装 Jupyter Notebook；IPython 内核，以在 Conda 环境中执行 Jupyter 的 Python 后端；以及 Jupyter Notebook 扩展，以简化工具（如自动完成）的编码任务：
conda install -c conda-forge notebook
conda install -c conda-forge nb_conda_kernels
conda install -c conda-forge jupyter_contrib_nbextensions
conda install -c conda-forge bokeh
在您的终端中，使用您给它起的名称激活 ArcGIS 环境：
Conda activate esriENV
下载完成后，您可以再次从终端打开 Jupyter Notebook，以便逐步记录您的工作并评估输出结果：
jupyter notebook
笔记本将与 API 一起安装，并在新窗口中打开。
连接到 ArcGIS Python API
安装完成后，现在是时候登录了。首先，我将向您展示如何匿名登录；随后的部分将解释如何使用开发者或个人账户登录。
作为匿名用户连接到 ArcGIS Online
有许多方法可以免费或几乎不收费地访问 ArcGIS 的部分内容；即使没有 ArcGIS 账户，您也可以使用资源、创建地图并共享它们，但功能有限。您可以通过运行代码或查阅文档在 Jupyter Notebook 中探索这些选项。公共账户的另一个好处是它允许匿名访问。
您将首先在您的 Jupyter Notebook 中导入 Python API。在这个示例中，您将使用一个公共账户，因此只显示免费资源。没问题——这里有很多东西可以让您探索。要从gis模块导入 GIS，运行：
from arcgis.gis import GIS
gis = GIS()
使用凭据连接到 ArcGIS 用户账户
如果您有用户账户，可以通过两种方式之一登录：使用内置登录或使用 API 密钥。登录的优势在于您可以将创建的任何地图保存到您的账户中。
无论哪种方式，您首先需要使用终端导入 GIS：
from arcgis.gis import GIS
要使用内置登录，运行：
gis = GIS(username="someuser", password="secret1234")
要使用 API 密钥登录，请使用以下代码（我已经缩短了令牌）：
gis = GIS(api_key="Request YOUR KEY AND INSERT IT HERE",
referer="https")
当您连接凭据并运行单元格时，输出将是您组织的 ArcGIS Online 用户名。我的账户用户名是“datalchemy”；如果您创建一个账户，您需要用您自己的账户凭据替换这些内容：
gis = GIS(username="datalchemy", password="xxx")
gis
如果您拥有账户和登录 URL（列在您组织的设置中），可以通过密码保护方式登录，再次用您自己的凭据替换：
import arcgis
import getpass
username = input ("Username:")
password = getpass.getpass ("Password:")
gis = arcgis.gis.GIS("https://datalchemy.maps.arcgis.com", username, password)
运行单元格后，将弹出一个框，提示您输入用户名和密码。
现在，您可以开始探索一些图像层次。
探索图像层：城市热岛地图
在我们开始之前，先说一下命名约定。图 6-1 是一幅芝加哥地图，我称之为map1。
这个指定是任意的：您可以根据自己的喜好为地图命名，但重要的是在涉及多个地图的代码编写时避免混淆。通常会指定一个带有整数后缀的地图变量，或者简单地使用m。选择一种方法并保持一致。它应与map()函数有所区别。您可以通过在 Jupyter Notebook 中运行以下代码生成map1：
map1 =gis.map("Chicago,Illinois")
map1
这将使您的地图居中显示在伊利诺伊州的芝加哥。您可以通过缩放到特定位置或绘制多边形来设置地图覆盖的区域，称为其范围。
图 6-1. 芝加哥地图
注意
要访问 ArcGIS Online 的帮助功能，请输入：
gis?
如果您键入一个点并按 Tab 键，将显示一个下拉菜单，其中包含您可以点击以获取更多信息的附加属性：
gis = GIS.
现在您可以尝试将一些属性层叠到您的地图上。图像层显示来自图像服务的数据。这些资源允许您应用规则来显示图像。您可以通过将项目类型标识为Imagery Layer来搜索它们。
您可以使用gis.content.search()中的搜索词来描述您想在地图中包含的要素。创建一个分配给gis.content.search()的变量可以帮助您查找可用的要素图层或其他item_type数据。使用以下代码片段，您可以探索可用的 Landsat 9 卫星视图，以用作影像层。不同的影像层具有不同的属性。在这里，我已将搜索限制为图 6-2 中显示的两个选项，这些选项是带有指定术语的公开可用要素图层：
from IPython.display import display
items = gis.content.search("Landsat 9 Views", item_type="Imagery Layer",
max_items=2)
for item in items:
display(item)
图 6-2. 使用gis.content.search()搜索 Landsat 9 视图
代码块中的item_type包括 Web 场景、要素图层、地理空间数据、底图和具有地形的交互式 3D 环境。您可以在文档中了解更多关于项目和项目类型的信息。
Web 场景允许您可视化和分析包含所有配置设置的地理空间内容，例如底图、样式和范围。要探索这些内容，您可以用web scene替换item_type。
我们将使用 Web 场景来探索芝加哥的热岛。根据美国环保局（EPA），热岛是：
城市化区域的温度比外围区域高。建筑物、道路和其他基础设施等结构比森林和水体等自然景观更多地吸收并重新发射太阳的热量。在这些结构高度集中且绿化有限的城市区域中，相对于外围区域，成为更高温度的“岛屿”。
城市热岛严重性测量城市基础设施如何吸收并重新发射这些热量。
要定位图像服务器的 URL 并开始在地图上探索热岛，请单击搜索结果中的标题（“城市热岛”），或通过在下一个单元格中调用img_svc_url来列出：
img_svc_url = 'https://server6.tplgis.org/arcgis6/rest/services/
heat_severity_2019/ImageServer'
这输出：
[<Item title:"Urban Heat Island Severity for U.S. cities - 2019" type:
Imagery Layer owner:TPL_GIS_Support>,
<Item title:"Multispectral Landsat" type:Imagery Layer owner:esri>]
您也可以使用此代码尝试栅格函数。更多关于这些内容的信息在下一节中，但现在，请尝试这个：
from arcgis.raster import ImageryLayer
landsat_urbanheat = ImageryLayer(img_svc_url)
landsat_urbanheat.properties.name
这输出：
‘Heat_severity_2019’
您已经访问了城市热岛严重性影像层。现在将其添加到您的芝加哥地图中：
map = gis.map('Chicago', zoomlevel=13)
map
map.add_layer(landsat_urbanheat)
landsat_urbanheat图层现在已添加到您的地图中，这应与图 6-3 类似。
图 6-3. 添加到地图范围的地图图层
现在，您将从多光谱 Landsat 图层的 URL 中添加影像层，并将其分配给变量landsat_ms：
img_svc_url ='*https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/
ImageServe*'
from arcgis.raster import ImageryLayer
landsat_ms = ImageryLayer(img_svc_url)
landsat_ms.properties.name
最后一行代码仅确认了图层的名称。
这将生成以下输出：
'Landsat/MS'
图像层的描述也可通过一行代码获取：
landsat_ms.properties['description']
输出是详细的：
Multispectral Landsat image service covering the landmass of the World. This
service includes scenes from Landsat 8 and Global Land Survey (GLS) data from
epochs 1990, 2000, 2005 and 2010 at 30 meter resolution as well as GLS 1975
at 60 meter resolution. GLS datasets are created by the United States
Geological Survey (USGS) and the National Aeronautics and Space Administration
(NASA) using Landsat images. This service can be used for mapping and change
detection of urban growth, change of natural resources and comparing Landsat
8 imagery with GLS data. Using on-the-fly processing, the raw DN values are
transformed to scaled (0 - 10000) apparent reflectance values and then
different service based renderings for band combinations and indices are applied.
The band names are in line with Landsat 8 bands; GLS data band names are mapped
along the same lines.
要探索更多的栅格函数对象，请访问常见数据类型文档。在学习栅格函数和图像图层的过程中，我还建议您随时查阅ArcGIS Python API 参考指南。
栅格函数
栅格函数处理来自图像的像素以进行一项或多项光栅处理，无需下载或创建中间文件进行分析。访问栅格函数受到限制，需要额外的凭据，但您仍然有机会探索 Landsat 多光谱图像。
这次，我们将探索洛杉矶地图，逐一了解其中的一些亮点：
from arcgis.gis import GIS
from arcgis.geocoding import geocode
from arcgis.raster.functions import *
from arcgis import geometry
import ipywidgets as widgets
import pandas as pd
您将在 gis.content.search() 中再次运行搜索。这次，您要查找组织外的图像（如果您已创建组织），并设置 outside_org=True。这是在访问公开可用数据集时的设置：
landsat_item = gis.content.search("Landsat Multispectral tags:'Landsat on AWS',
'landsat 9', 'Multispectral', 'Multitemporal', 'imagery', 'temporal', 'MS'",
'Imagery Layer', outside_org=**True**)[0]
landsat = landsat_item.layers[0]
df = **None**
您的搜索结果应包括多光谱 Landsat，如图 6-4 所示。要查看 landsat_item，请在单元格中调用它：
landsat_item
图 6-4. 显示 Esri 的多光谱 Landsat 成像图层的搜索结果
欲了解更多信息，请单击链接或运行以下代码：
from IPython.display import HTML
HTML(landsat_item.description)
本图展示了多光谱光带，这些光带由比人眼更敏感的仪器检测到。正如您在第一章 中学到的，多光谱光带通过电磁波谱中特定波长范围内的图像数据来突出显示不同的地表覆盖特征。您可以检查每个波段的名称、最小和最大波长以及其他属性：
pd.DataFrame(landsat.key_properties()['BandProperties'])
这些波段指示了与陆地覆盖、植被及其他类型地表有关的详细信息。表 6-1 是输出的摘录，以表格形式呈现以便阅读。
表 6-1. 多光谱和热红外传感器波段
| 波段 | 波长（微米） | 分辨率（米） |
|---|---|---|
| 第一波段—近岸气溶胶 | 0.43–0.45 | 30 |
| 第二波段—蓝色 | 0.45–0.51 | 30 |
| 第三波段—绿色 | 0.53–0.59 | 30 |
| 第四波段—红色 | 0.64–0.67 | 30 |
| 第五波段—近红外（NIR） | 0.85–0.88 | 30 |
| 第六波段—短波红外（SWIR）1 | 1.57–1.65 | 30 |
| 第七波段—短波红外（SWIR）2 | 2.11–2.29 | 30 |
| 第八波段—全色 | 0.50–0.68 | 15 |
| 第九波段—卷云 | 1.36–1.38 | 30 |
| 第十波段—热红外（TIRS）1 | 10.6–11.19 | 100 |
| 第十一波段—热红外（TIRS）2 | 11.50–12.51 | 100 |
使用您的 Jupyter Notebook 创建地图：
landsat = landsat_item.layers[0]
landsat
m = gis.map('los angeles')
m
现在添加要素图层：
m.add_layer(landsat)
您现在可以查看您的地图（图 6-5）。您还可以开始使用栅格函数。
图 6-5. 使用多光谱 Landsat 数据创建的洛杉矶地图
在这里，您将使用一个栅格函数来突出显示地图中的彩色可视化效果。例如，color infrared显示鲜艳的红色条带来指示健康的植被，而natural color显示我们通常看到的地形：绿色的植被，蓝色的水域，棕色的土壤。
许多选项可通过*动态范围调整（DRA）*来进行调整，以增强细节并在匹配您的计算机或其他设备的范围时提高可见性。您可以使用以下代码生成可用栅格函数的列表，如随后显示的输出所示：
for rasterfunc in landsat.properties.rasterFunctionInfos:
print(rasterfunc.name)
这将输出一个列表：
Agriculture with DRA
Bathymetric with DRA
Color Infrared with DRA
Geology with DRA
Natural Color with DRA
Short-wave Infrared with DRA
Agriculture
Bathymetric
Color Infrared
Geology
Natural Color
Short-wave Infrared
NDVI Colorized
Normalized Difference Moisture Index Colorized
NDVI Raw
NBR Raw
None
熟悉其中几个函数应该为您在独立探索中深入了解其他选项铺平了道路。
让我们尝试查看color_infrared。您将使用apply函数：
From arcgis.raster.functions import apply
Color_infrared = apply (landsat, 'Color Infrared with DRA')
要可视化地图，请调用map函数：
m = gis.map('los angeles')
m.add_layer(color_infrared)
m
健康的植被现在以鲜艳的红色显示（图 6-6）。
图 6-6. 用color_infrared查看的洛杉矶
用于查看植被的最佳栅格函数通常是归一化差异植被指数（NDVI）。图 6-7 展示了NDVI_colorized函数，显示绿色的植被：
ndvi_colorized = apply(landsat, 'NDVI Colorized')
ndvi_colorized
图 6-7. 用于 NDVI 的洛杉矶
您已经学会了如何搜索图像和要素图层，以及如何对特定位置应用栅格函数来突出显示特定功能。接下来，我们将探讨属性，并查看arcgis.geometry模块。
探索图像属性
在这一部分中，您将使用 ArcGIS 栅格模块get_samples来探索您的洛杉矶地图的几何形状。该模块允许您查看所选几何形状的样本点位置、空间分辨率和像素值。它是ArcGIS REST API的一个示例，该 API 提供有关 Web 应用程序架构的信息。与主要关注访问点的一般 ArcGIS API 文档不同，ArcGIS REST API 文档提供了地理空间功能的全部范围。它还描述了get_samples操作可以调用的参数：
get_samples(geometry, geometry_type=None, sample_distance=None, sample_count=None,
mosaic_rule=None, pixel_size=None, return_first_value_only=None,
interpolation=None, out_fields=None, slice_id=None)
在这里，您指定地图的几何形状和样本计数，并根据您定义的范围返回有关其图像属性的信息―在本例中是洛杉矶的范围：
import arcgis
g = arcgis.geometry.Geometry(area['extent'])
为什么要使用几何形状？当您创建 3D 模型时，像高度和太阳方位角这样的数据点对于计算阴影坡度非常有用。在使用hillshade函数时，基本上是将 2D 表面渲染为逼真的 3D 地形。
接下来，指定几何形状和样本计数：
samples = landsat.get_samples(g, sample_count=50,
out_fields='AcquisitionDate,OBJECTID,GroupName,
Category,SunAzimuth,SunElevation,CloudCover')
您可以指定要查看的样本列表中的项目，或查看所有项目：
samples[10]
输出包括位置、对象 ID、云层覆盖计算和像素值等信息：
{'location': {'x': -13150297.20625444,
'y': 4059732.8562727477,
'spatialReference': {'wkid': 102100, 'latestWkid': 3857}},
'locationId': 10,
'value': '1158 991 843 769 2850 1675 1030 44 21824 22410 22691',
'rasterId': 3508198,
'resolution': 30,
'attributes': {'AcquisitionDate': 1647023302000,
'OBJECTID': 3508198,
'GroupName': 'LC09_L1TP_041036_20220311_20220311_02_T1_MTL',
'Category': 1,
'SunAzimuth': 144.81874084,
'SunElevation': 45.83943939,
'CloudCover': 0.0024},
'values': [1158.0,
991.0,
843.0,
769.0,
2850.0,
1675.0,
1030.0,
44.0,
21824.0,
22410.0,
22691.0]}
在 Python 中，我们使用datetime对象进行时间序列数据或在不同时间点收集的数据。再次，在 Python 中，您正在对datetime类进行采样，并可以通过运行以下代码渲染获取日期。（我将在本章后面提供更多详细信息，但这是我们在示例数据上使用该功能的方式。）
由于 Python 基于零索引访问数据，使用[0]请求列表中的第一个索引或示例中的第一个值：
import datetime
value = samples[0]['attributes']['AcquisitionDate']
datetime.datetime.fromtimestamp(value /1000).strftime("Acquisition Date: %d %b,
%Y")
这将输出：
'Acquisition Date: 11 Mar, 2022'
通过在 Los Angeles 的特定位置采样值，您可以估算地图上特定点的光谱配置。
m = gis.map('los angeles')
m
m.add_layer(landsat)
如果选择特定像素，光谱配置将绘制在该位置反射的所有波段，如图 6-8 所示。
下一个代码允许您在 Jupyter Notebook 中的画布上选择一个点以识别光谱配置。请注意，如果您从 ArcGIS Python API 指南运行此示例并请求 Landsat 9，则会生成错误，因为添加了额外的波段。
图 6-8. 点击生成光谱配置在图 6-9
get_samples()方法将收集样本数据中包含的像素值。像素值是*数字。*为了正确计算它们，您需要先将其转换为浮点数，然后转换为整数。数字记录像素的电磁强度。
您已将Bokeh 包安装到 Conda 环境中以允许交互式绘图和制图。Landset 规格基于 Landsat 8 数据，并仅包括了八个波段。这导致了 Bokeh 的错误。为解决此问题，请使用包含新波段的运行：
from bokeh.models import Range1d
from bokeh.plotting import figure, show, output_notebook
from IPython.display import clear_output
output_notebook()
def get_samples(mw, g):
clear_output()
m.draw(g)
samples = landsat.get_samples(g, pixel_size=30)
values = samples[0]['value']
vals = [float(int(s)/100000) for s in values.split(' ')]
x = ['1','2', '3', '4', '5', '6', '7', '8','9','10','11']
y = vals
p = figure(title="Spectral Profile", x_axis_label='Spectral Bands',
y_axis_label='Data Values', width=600, height=300)
p.line(x, y, legend_label="Selected Point", line_color="red", line_width=2)
p.circle(x, y, line_color="red", fill_color="white", size=8)
p.y_range=Range1d(0, 1.0)
show(p)
print('Click anywhere on the map to plot the spectral profile for that location.')
m.on_click(get_samples)
您应该收到 Bokeh 2.4.2（或您的版本号）已成功加载的消息。现在返回地图并选择您的点。点击生成图 6-8 中显示的地图（您的地图将根据您选择的点有不同的值）。
您生成的地图（图 6-8）现在是交互式的。在新地图上任意点击以绘制该位置的光谱配置。现在您可以选择不同的点并查看生成的值。图 6-9 展示了在我选择的位置上用于探测卷云的第 9 波段的光谱配置。
图 6-9. 使用 Bokeh 可视化的光谱配置
根据点击 Landsat 9 图像图层时选择的独特位置，您的光谱配置文件将不同。请参阅表 6-1 以识别包含的波段。这些代表了 ArcGIS Python API 中get_samples方法的几个应用。
改进图像
光栅函数允许您执行诸如提取特定波段以检查土地利用、植被或火灾；连续数据（如温度）；扫描图像；以及卫星图像等操作。拉伸函数允许您调整地图上的亮度和对比度。以下代码从可见光谱中选择 3、2 和 1 波段（红、绿和蓝），并导入stretch raster.function：
from arcgis.raster.functions import stretch, extract_band
naturalcolor = stretch(extract_band(landsat, [3,2,1]),
stretch_type='percentclip', min_percent=0.1, max_percent=0.1,
gamma=[1, 1, 1], dra=True)
naturalcolor
此代码设置的百分比剪切最小值(percentclip)将排除拉伸应用于栅格的最低 10%值。如果大多数像素在特定范围内，这是非常有用的。设置stretch_type会修剪掉异常值，重新分布值的直方图，从而生成图 6-10 中的图像。
图 6-10. 自然色波段展示了肉眼在没有增强的情况下所看到的图像
图 6-10 中的图像是肉眼在没有任何增强的情况下所看到的。根据您要查找的特征，可能有其他波段可以探索。我们将在接下来的部分查看其他波段。
比较同一位置的多个时间点
ArcGIS API 允许您使用称为时间滑块的地图小部件比较同一位置在不同时间点的图像，该小部件使您可以使用可配置的属性（如开始和结束时间及其间隔）对地图进行动画化。对于此示例，我选择了一个城市和一个卫星底图。
地图缩放小部件接受 0 到 23 之间的值。级别 9 提供了全球概览；级别 10 介于大型都市区域和城市之间；级别 20 是个体建筑物的级别。我从未缩放到 23，但这样做会呈现出最详细的视图。然而，根据可视化参数的规模，缩放的限制程度有所不同（转换率可在文档中查看）。您可以更改缩放级别并处理数据，但分辨率将保持不变。例如，所有 Landsat 图像的分辨率为 15 米。这意味着卫星图像捕捉到的地面细节为 15 米或更大。（作为比较，洛杉矶的好莱坞标志和标准半挂拖车均长约 15 米。）
尝试将缩放级别调整到 10：
map = gis.map("los angeles")
map.basemap = "satellite"
map.zoom = 10
map
item title的名称将作为输出打印以供参考：
landsat_item = gis.content.search("Landsat Multispectral tags:'Landsat on AWS',
'landsat 9', 'Multispectral', 'Multitemporal', 'imagery', 'temporal', 'MS'",
'Imagery Layer', outside_org=**True**)[0]
print(landsat_item)
<Item title:"Multispectral Landsat" type:Imagery Layer owner:esri>
您已搜索到一个 Landsat 项，并且通过其标题在输出中确认了该项。不同的子层将提供不同的详细信息。
将图层添加到地图中：
map.add_layer(landsat_item)
请求确认map.time_slider是否可用：
map.time_slider
它应该返回输出True。
您可以将日期调整为特定间隔。以下代码测量所选开始和结束时间之间的 10 天间隔：
from datetime import datetime
map.time_slider = True
map.set_time_extent(start_time=datetime(2021, 12, 12), end_time=datetime(2022, 4,
12), interval=10, unit='days')
此代码将在 Figure 6-11 中输出地图，其中包括时间滑块小部件。
图 6-11. 时间滑块小部件
在代码单元中选择map.draw(polygon)选项，然后使用鼠标绘制一个多边形。如果您有感兴趣的特定位置，请输入其坐标：
from arcgis.geometry import Point
pt = Point({"x" : 34.092809, "y" : -118.328661,
"spatialReference" : {"wkid" : 3309}})
空间参考或*WKID（Well-Known ID）*定义了您选择的位置的公差和分辨率（参见空间参考列表）。
map.draw选项包括选择点（pt）、polyline或polygon：
map.draw(pt)
map.draw('polyline')
map.draw('polygon')
当您选择多边形并运行代码单元时，只需返回到地图，您就可以跟踪一个多边形。时间滑块小部件允许在感兴趣的区域跨时间序列进行比较。在处理图层时，通常需要筛选出出现在画布上的内容。
过滤图层
您可能希望过滤您的图层：例如，也许您只想看到具有少于 10%云覆盖率的特定集合中的 Landsat 图层。为此，熟悉 Landsat 的数据字典有助于获取关于数据库或系统中数据元素、属性和名称的完整信息。
在以下代码片段中，WRS_Row提供了有关拍摄图像的卫星轨道路径的信息：
selected = landsat.filter_by(where="(Category = 1) AND (CloudCover <=0.10) AND
(WRS_Row = 36)",
geometry=arcgis.geometry.filters.intersects(area['extent']))
值36表示北半球。（您可以在 Landsat 数据字典、集合背景或其他 Landsat 信息中找到这些详细信息。）
您可以访问的信息显示在您之前生成的 HTML 项目描述中。现在您可以在数据帧（df）中查看刚刚创建的表格：
fs = selected.query(out_fields="AcquisitionDate, GroupName, Best, CloudCover,
WRS_Row, Month, Name",
return_geometry=True,
return_distinct_values=False,
order_by_fields="AcquisitionDate")
因为您正在比较日期，您将希望看到最旧的获取日期以及最近的获取日期。要清楚地查看输出，请在代码单元中运行它。重要的是看列标题以及在哪里找到标识信息（见 Figure 6-12）。您可以通过调用df.head()来查看前五行：
df = fs.sdf
df.head()
图 6-12. 采集日期输出
现在运行df.tail()并选择最近的获取日期。您可以在 Jupyter Notebook 中运行代码：
df = fs.sdf
df.tail()
以下命令给出了数据的形状：
df.shape
输出告诉我们，数据有 9 列和 193 行：(193,9)。
如果您只想要获取日期，可以运行此数据，或从数据帧中的任何其他列中获取：
df['Time'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df['Time'].head(10)
这个输出：
0 1977-05-30
1 1979-06-08
2 1989-06-28
3 1990-05-07
4 2000-04-24
5 2000-05-01
6 2005-05-15
7 2005-05-24
8 2009-04-16
9 2009-05-11
Name: Time, dtype: datetime64[ns]
波段具有特定波长，用于统计计算。在图像图层中要避免像素重叠，因为这可能会扭曲计算结果。您可以使用默认方法，其中像素值是从最后一个数据集计算的（见 图 6-13）：
m3 = gis.map('los angeles', 7)
display(m3)
m3.add_layer(selected.last())
图 6-13. 使用 last() 方法重叠像素
或者，您可以请求 first() 方法，并从第一个光栅数据集计算像素值（输出显示在 图 6-14 中）：
m3 = gis.map('los angeles', 7)
display(m3)
m3.add_layer(selected.first())
图 6-14. 使用 first() 方法重叠像素
通过 OBJECTID 查询数据，辅以 df(head) 和 df(tail) 数据，您可以比较不同时期捕获的图像：
old = landsat.filter_by('OBJECTID=3106827')
new = landsat.filter_by('OBJECTID=3558253')
from arcgis.raster.functions import *
stretch 函数 通过扩展像素值来提高可见性。有不同类型的拉伸，但在这里，您只需通过计算标准偏差 (stddev) 删除任何极端值：
diff = stretch(composite_band([ndvi(old, '5 4'),
ndvi(new, '5 4'),
ndvi(old, '5 4')]),
stretch_type='stddev', num_stddev=2.5, min=0,
max=255, dra=True, astype='u8')
diff
这将输出 图 6-15。绿色条带表示植被密度增加，洋红色显示所选 OBJECTID 时间范围内的减少。
图 6-15. ArcGIS 函数显示随时间图像差异
或许有您想要捕获的特定阈值。尝试仅测量阈值变化高于 10%的区域。使用以下代码：
threshold_val = 0.1
masked = colormap(remap(ndvi_diff,
input_ranges=[threshold_val, 1],
output_values=[1],
no_data_ranges=[-1, threshold_val], astype='u8'),
colormap=[[1, 124, 252, 0]], astype='u8')
Image(masked.export_image(bbox=area['extent'], size=[1200,450], f='image'))
输出（图 6-16）将这些区域呈现为绿色。
图 6-16. 植被指数（NDVI）变化的掩码阈值
让我们再次渲染地图：
m = gis.map('los angeles')
m
最终图像显示在 图 6-17 中。您现在可以通过向地图添加图层来查看掩码阈值：
m.add_layer(diff)
m.add_layer(masked)
图 6-17. 在请求的阈值处显示的组合图像，展示了掩码阈值
摘要
现在，您已经根据 ArcGIS API for Python 定制并探索了一些公开可用资源，应该对栅格函数和图像图层有所了解。继续学习时，请尝试独立使用这些函数和图层，看看您可以访问哪些其他查询。
第七章：GeoPandas 和空间统计
地图很美丽。它们讲述的故事可能如此迷人，以至于很容易无意中忽略其中的地理空间统计数据。但地理空间地图不仅仅是静态图像。其中嵌入了信息，例如与 GIS 图层中特定要素相关联的属性或在光栅图像中观察到的像素密度。
Python 有多种用于处理地理空间数据的包。如果您熟悉 Python，您可能已经了解 pandas，这是一个专为 Python 构建的数据分析工具。Pandas 允许我们将各种数据类型读入数据框中：一组包含行（表示记录）和列（表示属性）的表格。GeoPandas 是 pandas 的扩展，允许您使用其称为 GeoDataFrame 的方式操作几何和地理空间数据框架：每行是空间要素，如点、线或多边形。
本章将向您展示如何使用 GeoPandas 和 GeoDataFrame 分析数据并创建地图，以及一些其他重要的包，如 matplotlib 用于可视化数据和人口普查数据 API。您还将学习如何访问地理空间文件，并通过创建和比较人口地图深入了解人口统计数据。
安装 GeoPandas
要安装 GeoPandas，您将在终端中使用 conda 和 conda-forge。与前几章一样，我还将向您展示如何为本章中所有需要处理的数据创建环境。
安装 GeoPandas 可能有点棘手。可以把它想象成一个情绪化的青少年。它通常想要先行（导入）并喜欢最新的趋势（注意依赖项的版本）。它还喜欢一个响亮的昵称：通常情况下，GeoPandas 被称为 gpd，但有时代码会写出 geopandas。请注意并调整您的变量。
要设置您的环境，请从这里开始：
conda create --name geop_env python = 3.X
您可以在此处添加您的 Python 版本。
MacBook-Pro-8:~ USERNAME$ conda activate geop_env
现在您可以开始将文件安装到 geop_env 中：
conda install -c conda-forge geopandas
现在您可以将任何想要访问的包添加到 geop_env 中。如果您希望获取一组用于地理空间分析的包，我建议下载地理空间包，但现在如果您愿意，也可以逐个添加包：
conda install -c conda-forge geospatial
GeoPandas 简化了在 Python 中处理地理空间数据的工作。让我们探索扩展操作，使您能够在几何数据上执行空间操作。
处理 GeoJSON 文件
GeoPandas 让您能够处理 GeoJSON 文件。您可能已经听说过 JSON，这是一种常见的数据交换格式。GeoJSON 是一种基于 JSON 的数据交换格式，专门设计用于表示地理数据结构及其非空间属性。使用 GeoJSON 文件工作使得识别地图中的坐标和位置变得相当简单。
图 7-1 显示了来自geojson.io的地图。Geojson.io 是一个开源工具，简化了创建、查看和分享地图。你可以在互联网上执行各种功能，而无需将任何数据下载到你的计算机。
图 7-1. 缩放到巴黎第一区并选择边界以创建 GeoJSON 文件
点击geojson.io链接并放大到一个国家。我选择了法国，因为，嗯，它是法国。查找特定位置的坐标的简单方法是选择你感兴趣的区域并拖动多边形周围的区域。选择多边形的工具可见于地图仪表板的右边缘。尝试在地图上找到法国巴黎。找到杜乐丽花园并在其周围绘制一个多边形。保存多边形后，你可以通过点击边界查看格式化选项和样式属性。
顶部菜单中的 Meta 标签（显示在图 7-1 中）允许你生成一个边界框（bbox），其中包含你想要映射的坐标和点集。它定义了感兴趣对象的空间位置和其坐标。你还可以加载以 well-known text（WKT）标记语言表示的字符串，用于矢量几何对象。
左下角的按钮，在图 7-2 中可见，允许你在 Mapbox（矢量瓦片）、卫星（光栅数据）和 OSM（OpenStreetMap，在第五章中介绍过）之间切换。
图 7-2. 地图查看器选项
如果你想要保存这个文件，你可以以几种不同的格式保存，包括*.csv*、.txt和 shapefile。然后你可以将其导入到 QGIS、ArcGIS 或你选择的平台进行进一步的探索和分析。电子表格数据存储为*.csv或.txt*，通常包含位置信息，如经度和纬度的坐标数据、地址或邮政编码。Shapefile 是 Esri 的矢量文件，捕捉位置、形状和有关位置的相关文件。
创建一个 GeoDataFrame
GeoDataFrame 类似于带有额外几何数据列的 pandas 数据帧。许多类型的矢量数据描述了具有固定位置的离散数据，而连续数据通常由栅格数据表示（尽管两者都可以是任意类型）。
GeoPandas 提供了一些数据集示例供探索。你将访问缩写为nybb的 NYC Boroughs 示例：
geopandas.datasets.available = ['naturalearth_cities', 'naturalearth_lowres',
'nybb']
这些文件代表城市的位置、不同国家的轮廓以及纽约市行政区的边界。要访问它们，你需要将数据集名称作为参数包含进去：
geopandas.datasets.get_path("nybb")
您可以从终端打开 Jupyter Notebook：
(geop_env) MacBook-Pro-8:~ USERNAME$ jupyter notebook
现在，笔记本已打开，您可以创建一个数据帧。首先，将您需要的软件包导入到您的环境中：
%matplotlib inline
import requests
import pandas as pd
import geopandas as gpd
from scipy.spatial.distance import cdist
接下来，您需要检索数据集并调用活动几何列（在本例中为“geometry”）来创建一个名为world.geometry.name的数据帧。当您调用head()时，将返回数据集中的前五行和列：
boros_world = gpd.read_file(gpd.datasets.get_path('nybb'))
print(f"{type(boros_world)}, {boros_world.geometry.name}")
print(boros_world.head())
print(boros_world.geometry.geom_type.value_counts())
您应该得到以下输出：
<class 'geopandas.geodataframe.GeoDataFrame'>, geometry
BoroCode BoroName Shape_Leng Shape_Area \
0 5 Staten Island 330470.010332 1.623820e+09
1 4 Queens 896344.047763 3.045213e+09
2 3 Brooklyn 741080.523166 1.937479e+09
3 1 Manhattan 359299.096471 6.364715e+08
4 2 Bronx 464392.991824 1.186925e+09
geometry
0 MULTIPOLYGON (((970217.022 145643.332, 970227....
1 MULTIPOLYGON (((1029606.077 156073.814, 102957...
2 MULTIPOLYGON (((1021176.479 151374.797, 102100...
3 MULTIPOLYGON (((981219.056 188655.316, 980940....
4 MULTIPOLYGON (((1012821.806 229228.265, 101278...
MultiPolygon 5
dtype: int64
注意第一列中的索引分配。通过引用其索引位置，您将能够隔离单个区域。BoroCode和BoroName是您可以指定的附加标识符，但是您可能会注意到在随后的代码片段中通过索引调用的简易性。
您现在正在使用一个 GeoDataFrame。您可以在 Jupyter Notebook 中运行以下代码生成图 7-3 的绘图：
boros = gpd.read_file(gpd.datasets.get_path('nybb'))
ax = boros.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
图 7-3。纽约市行政区域的位置
使用figsize参数可以更改图像的大小。在示例中，alpha 值小于 1 可调整透明度。边缘颜色也可以自定义。你可以在matplotlib 文档中找到详细信息。
在大多数情况下，我更喜欢没有可见轴的地图。要删除框架，只需添加：
ax.set_axis_off()
您将通过交互地图进一步探索数据，因此不需要区域坐标。这将生成图 7-4 中的图像。
图 7-4。没有轴框架的纽约市行政区域的位置
通过调用行号可以请求单个区域。斯泰顿岛在您几页前生成的 GeoDataFrame 中的第 0 行。您可以使用可调用函数DataFrame.loc根据其索引号访问位置。我认为loc是“列标签”的缩写，因为它需要输入标签，如“索引”或“几何”。缩写iloc用于整数位置索引/列的整数。替换您创建的数据帧名称：
boros.loc[0,'geometry']
这将输出显示在图 7-5 中的斯泰顿岛的图像。
图 7-5。使用.loc函数调用单个区域
您可以在图 7-5 中看到斯泰顿岛区域的边界，但是现在您知道如何访问 GeoDataFrame 后，获取更多信息会更有帮助。您还可以找到区域的面积：
gdf = gdf.set_index("BoroName")
gdf["area"] = gdf.area
gdf["area"]
这将输出：
BoroName
Staten Island 1.623822e+09
Queens 3.045214e+09
Brooklyn 1.937478e+09
Manhattan 6.364712e+08
Bronx 1.186926e+09
Name: area, dtype: float64
要查看一个行政区域的边界，请运行：
gdf['boundary'] = gdf.boundary
gdf['boundary']
这将输出：
BoroName
Staten Island MULTILINESTRING ((970217.022 145643.332, 97022...
Queens MULTILINESTRING ((1029606.077 156073.814, 1029...
Brooklyn MULTILINESTRING ((1021176.479 151374.797, 1021...
Manhattan MULTILINESTRING ((981219.056 188655.316, 98094...
Bronx MULTILINESTRING ((1012821.806 229228.265, 1012...
Name: boundary, dtype: geometry
在交互式地图 (Figure 7-6) 中，您可以通过悬停在市区上来检索信息。参考距离是 0，因为史泰登岛的索引号是 0. 点击任何其他市区将告诉您与史泰登岛的距离。gdf.explore 函数生成一个交互式地图：
gdf.explore("area", legend=True)
图 7-6. 纽约市行政区交互地图
将像几何形状这样的地理空间能力添加到 pandas 数据框中仅仅是个开始。人口普查数据是开放数据的宝库，了解如何使用 GeoPandas 将帮助您扩展访问和查询这一重要数据集的能力。让我们再看看更多人口普查数据。
使用美国人口普查数据：洛杉矶人口密度地图
在这个练习中，您将计算加利福尼亚州洛杉矶县各个小区的人口密度。
警告
当使用 Python 人口普查包时要小心。如果开发者失去了维护它们的兴趣或者被其他任务压倒，您的地图和图形将不再工作。在写这一章节时，我就是这样学到的。好消息是，它提醒我始终可靠的备用基础知识。
如果您需要一个专业服务或研究的人口普查数据问题的高质量地图，您可以在像 QGIS 这样的开源平台上或使用像 ArcGIS 这样的专有软件构建一个。这个例子将带您了解基础知识。一旦您理解了人口普查 API 的结构并能够在一个不花哨的 Python 包之外编写代码，您就可以轻松访问这个最强大的公开可用数据集之一。
第一步是在 census.gov 上找到数据。根据您的兴趣领域，2020 年人口普查的数据已经可以用于分析。我依赖于 Census Reporter 来识别和下载数据集。大多数仓库都有 ACS 数据（通常称为“长表人口普查”），每年从全国一部分收集，并以年度或五年时段的形式捆绑。API 默认获取最近五年的 ACS 时段。这是一个易于导航的第一站，全面了解捕获数据类型的方法。
Census Reporter 为许多大都市区提供即用的社区数据，这些示例中您将看到这些数据。
小贴士
对于更高级的用户，IPUMS（前身为集成公共使用微数据系列，现以其首字母缩写命名）提供了来自人口普查和调查数据的全球协调数据集，并附有支持文件。您可以使用它将 ACS 和人口普查数据与其他来源的合并数据集集成。
通过人口普查 API 和 FTP 访问小区和人口数据
通过 FTP 站点或直接访问人口普查 API 是能够在人口普查 GUI 外工作的基础。这将帮助您解决更大的问题，并且数据已准备好在笔记本或 QGIS 中进行分析。
FTP 站点提供了快速简便访问大文件的途径，这些文件可能难以导航。如果您在 MacOS 环境下工作，可能需要调整一些设置，但我更喜欢通过 FTP 访问人口普查数据的可靠性和无杂乱。
选择“系统偏好设置 >> 共享”并选中“远程登录”。选择允许访问的用户。
接下来，导航到您的 Finder 窗口并选择 Go。在下拉菜单中，您将看到一个连接到服务器的选项（图 7-7）。
图 7-7. 连接到美国人口普查局 FTP 站点
FTP 站点 可通过链接访问。您将能够浏览可用的文件。您还可以通过服务器站点获得访问权限。总结文件是大多数感兴趣的文件所在之处。
您将要探索美国人口普查的十年人口数据。这里包含的地理数据来自 2019 年 ACS 一年调查。人口普查表来自 2010 年人口普查。^(1) Google Colab 在这项练习中表现良好，因为文件结构是直接可用的，并且会话结束后，文件会被移除。
您需要下载以下在 Colab 中尚未包含的包：
!pip install geopandas #install separately to avoid errors
import requests # accessing from the CensusAPI **import** sys # checking python version **import** pandas **as** pd # working with dataframes **import** geopandas **as** gpd # creating maps **import** folium **as** fm # creating HTML map complete with OSM basemap
当处理人口普查数据时，我学到了一条重要的经验教训：始终了解您的包版本，以防存在依赖关系（或冲突！）。这是一个有用的实践方法，可以添加到您的笔记本中。使用以下代码显示版本：
# Display package versions
print("Python Version", sys.version)
print("requests version:", requests.__version__)
print("pandas version:", pd.__version__)
print("geopandas version:", gpd.__version__)
print("folium version:", fm.__version__)
现在您可以从 FTP 服务器请求人口普查区信息。您正在寻找人口普查区的形状文件。这些是 .zip 文件；请注意，您需要的文件可能捆绑在一个较大的文件中。下载加利福尼亚州的拓扑集成地理编码和参照（TIGER）文件：FIPS 06.（联邦信息处理标准出版物，或 FIPS 文件指定州和县的等效物）。TIGER 形状文件包含地理信息但不包括人口统计信息，因此我们将很快添加该信息。链接标识了美国人口普查局 FTP 站点，并指定了加利福尼亚州 2019 年按区划级别的地理文件：
'ftp://ftp2.census.gov/geo/tiger/TIGER2019/TRACT/tl_2019_06_tract.zip'
一旦您在 FTP 中确定文件，要检索它，您需要 !wget 调用。wget 实用程序从 Web 服务器检索 FTP 文件：
!wget ftp://ftp2.census.gov/geo/tiger/TIGER2019/TRACT/tl_2019_06_tract.zip
对于 !wget 调用，.zip 文件路径上的引号是不需要的，并且文件将解压缩到您的 Colab 笔记本中。您可以在笔记本中直接观察解压缩进度，直到完成到 100%。您的输出应如下所示：
--2022-08-22 15:52:00-- ftp://ftp2.census.gov/geo/tiger/TIGER2019/TRACT/
tl_2019_06_tract.zip
=> 'tl_2019_06_tract.zip'
Resolving ftp2.census.gov (ftp2.census.gov)... 148.129.75.35, 2610:20:2010:a09:
1000:0:9481:4b23
Connecting to ftp2.census.gov (ftp2.census.gov)|148.129.75.35|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done. ==> PWD ... done.
==> TYPE I ... done. ==> CWD (1) /geo/tiger/TIGER2019/TRACT ... done.
==> SIZE tl_2019_06_tract.zip ... 29388806
==> PASV ... done. ==> RETR tl_2019_06_tract.zip ... done.
Length: 29388806 (28M) (unauthoritative)
tl_2019_06_tract.zi 100%[===================>] 28.03M 262KB/s in 69s
2022-08-22 15:53:10 (416 KB/s) - 'tl_2019_06_tract.zip' saved [29388806]
图 7-8 显示文件在解压缩到文件夹层次结构后的情况。导航到县级数据，使用与查找加利福尼亚州小区相同的步骤。
图 7-8。Google Colab 文件结构
现在文件解压缩后，我们可以从 Web 服务器检索 FTP 文件：
'ftp://ftp2.census.gov/geo/tiger/TIGER2019/COUNTY/tl_2019_us_county.zip'
!wget ftp://ftp2.census.gov/geo/tiger/TIGER2019/COUNTY/tl_2019_us_county.zip
现在您可以使用 GeoPandas 绘制加利福尼亚小区（ca_tracts），从您的 Colab 窗口读取文件。每个文件名旁边有三个点。右键单击点并选择复制路径以复制文件链接：
ca_tracts = gpd.read_file('/content/tl_2019_06_tract.shp')
ca_tracts.plot()
输出显示在图 7-9 中。
图 7-9。加利福尼亚人口普查小区
让我们花点时间弄清楚表 ID 的含义，以帮助您找到正确的数据。
在浏览器中访问人口普查 API 的数据
接下来，您需要打开美国人口普查局网站。在搜索窗口中输入**DP03**，这是人口普查数据的表 ID（图 7-10）。该 ID 可能看起来很随机，但实际上它提供了重要信息。从结果中可以立即看出，DP 指示一个数据概况表。表 ID 描述此表包含“广泛的社会、经济、住房和人口统计信息”。ID 中的后续数字 03 还涉及选定的经济特征。
图 7-10。人口普查 API：数据
选择 Enter 以查看可用表格列表。图 7-11 显示来自此搜索可用的不同年份和类型的 ACS 调查。这些被称为数据概况，稍后我们将再次提及它们。
图 7-11。来自美国人口普查局的数据概况表
census.gov 上的可用 API 页面（图 7-12）是一个重要资源，包含广泛主题的公开可用数据集。从这里选择美国社区调查（ACS）。
图 7-12。可用 API 及其他美国人口普查局数据集
在美国社区调查（ACS）下，选择 1 年估计，并滚动到 2019 年。滚动表格列表直至看到数据概况。复制第一个链接并粘贴到浏览器中。您将编辑此链接并生成一个*.csv*文件。
提示
当你在“可用 API”页面时，你可以请求一个用于运行美国人口普查局服务请求的密钥。对于这个练习，你不需要一个，但我建议请求一个并将其保存在一个可靠的位置。浏览列出的其他数据集，并熟悉公开可用数据的范围。另一个有用的美国人口普查局资源是开发者论坛，这是一个接受数据集和咨询请求的活跃社区。
使用数据概况
在数据概况部分，你有将要编辑的示例调用以及关于数据概况表的其他详细信息：
-
示例调用：
api.census.gov/data/2019/acs/acs1/profile?get=group(DP02)&for=us:1&key=YOUR_KEY_GOES_HERE
点击链接以收集更多信息。你的第一个编辑将包括替换你感兴趣的表格 DP03，并从 URL 中删除占位符 &key=YOUR_KEY_GOES_HERE：
api.census.gov/data/2019/acs/acs1/profile?get=group(DP03)&for=us:1
然后将 URL 粘贴到浏览器中并按 Enter 键。现在你将想要定制 URL 的地理部分，包括所有都市区域。从前一页选择示例和支持的地理位置，并滚动到地理层次结构图 7-13 中所示的都市区域。
图 7-13. 地理层次结构
注意到 *（通配符选项），复制 URL，确保再次删除 YOUR_KEY_GOES_HERE 部分。另一种选项将为你提供一个单一都市区域的数据：俄亥俄州阿克伦都会区。
一旦你熟悉了 URL 和数据表，只需进行一些快速编辑即可获得所需的人口或经济数据的可用文件。按如下修改 URL：
现在将其粘贴到浏览器中并按 Enter 键。该文件以 JSON 格式保存在我的 Mac 上。与其将其保存为 .csv 文件，这会迫使你在 Microsoft Excel 中进行编辑，不如运行以下代码行将其从 .json 转换为 .csv：
df = pd.read_json(‘link to json file’)
df.to_csv()
df
这将输出一个格式化的 .csv 文件。
创建地图
要完成我们创建人口密度项目，你将返回到 2010 年人口普查的 URL。
让我们从查看 2010 年人口普查的参考 URL 开始：
现在您已经熟悉了人口普查 API 的调用，您可以在 Python 中分解链接，计算密度并自定义表格。在加利福尼亚州人口普查区段的 URL 调用中包括表P001001；sf1指的是十年一次的人口普查：
# The built out request for the URL https://api.census.gov/data/2010/dec/
sf1?get=LSAD_NAME,P001001&for=tract:*&in=state:06 HOST = "https://api.census.gov/data" year = "2010" dataset = "dec/sf1" base_url = "/".join([HOST, year, dataset])
predicates = {}
get_vars = ["LSAD_NAME", "P001001"]
predicates["get"] = ",".join(get_vars)
predicates["for"] = "tract:*" predicates["in"] = "state:06" r = requests.get(base_url, params=predicates)
现在您已经发出了请求，让我们查看数据的列名称：
print(r.json()[0])
['LSAD_NAME', 'P001001', 'state', 'county', 'tract']
记住index函数。您只需读取数据的第一行（r）来查看标题。您还可以创建新的名称：
# Create user friendly column names tract_name = ["tract_name", "tract_pop", "state_fips", "county_fips",
"tract_fips"]
# Reading the json into pandas df tractdf = pd.DataFrame(columns=tract_name, data=r.json()[1:])
# Changing data types to integer tractdf["tract_pop"] = tractdf["tract_pop"].astype(int)
tractdf.head()
在笔记本中跟随操作。以下代码将选择所有洛杉矶县的属性，删除任何空值，并创建一个新的geoid列。运行代码，并验证确实已添加了geoid列。
首先，我们需要从洛杉矶县中选择属性：
onecounty_tractdf = tractdf.loc[tractdf['county_fips'] ==
'037'].copy()Onecounty_tractdf
然后，我们创建包含新geoid列的新数据框架：
onecounty_tractdf['tract_fips'] = onecounty_tractdf['tract_fips'].str.ljust(6,'0')
onecounty_tractdf['geoid'] = onecounty_tractdf['state_fips'] +
onecounty_tractdf['county_fips'] + onecounty_tractdf['tract_fips']
onecounty_tractdf
数据框中有 2,346 条记录。您还可以通过一行代码快速计算记录数量：
onecounty_tractdf.count()
要将区段级别和县级文件连接起来，您需要选择一个连接属性。要回顾可用的选择，请使用以下代码检查列标题：
ca_tracts.info
该连接将合并来自两个数据集的列标题GEOID和geoid：
attr_joined = pd.merge(ca_tracts, onecounty_tractdf, left_on='GEOID',
right_on='geoid')
# Check that all 2345 Census Tracts joined
attr_joined.count()
双重检查坐标参考系统，确保地理文件和人口普查数据匹配。ALAND 是面积，以平方米计。现在您可以将分组块添加到您的地图上：
map = fm.Map(location=[center_y, center_x], zoom_start=10)
# Add Study Area Block Groups to Map
fm.Choropleth(
geo_data = ca_prj,
data=ca_prj,
columns=['tract_pop','ALAND'],
key_on= 'feature.properties.tract_pop',
fill_color='YlGnBu',
name = 'Population Density',
legend_name='Population Density'
).add_to(map)
map
人口密度衡量在给定区域内居住的人数，通常以平方公里计算。您可以在代码单元格中直接将研究区域转换为平方公里：
# Create a new column for Census Tract area in square Kilometers ca_prj['AreaLandKM2'] = (ca_prj['ALAND'] * .000001)
ca_prj[['geoid','TRACTCE','ALAND','AWATER','AreaLandKM2']].head()
ca_prj['ppl_perKM2']=(ca_prj['tract_pop']/ca_prj['AreaLandKM2'])
ca_prj[['geoid','TRACTCE','tract_pop','AreaLandKM2','ppl_perKM2']].head(16)
接下来，选择您希望地图中心的位置：
center_x = (ca_prj.bounds.minx.mean() + ca_prj.bounds.maxx.mean())/2
center_y = (ca_prj.bounds.miny.mean() + ca_prj.bounds.maxy.mean())/2
print(f'The center of the data file is located at {center_x} {center_y}')
现在您可以通过运行以下代码查看地图和图例（自动生成），显示洛杉矶县内的人口密度：
map = fm.Map(location=[center_y, center_x], zoom_start=10)
# Add Study Area Block Groups to Map
fm.Choropleth(
geo_data = ca_prj,
data=ca_prj,
columns=['TRACTCE','ppl_perKM2'],
key_on= 'feature.properties.TRACTCE',
fill_color='YlGnBu',
name = 'Population Density',
legend_name='Population Density'
).add_to(map)
map
输出结果，您完成的地图显示在图 7-14 中。
图 7-14. 洛杉矶县人口普查区段密度
从人口普查 API FTP 站点导入 TIGER 文件可以保证单一的位置精度源。在这里，您查看了人口普查区段级别，但通过简单编辑 URL 调用，您可以更改到县级、国会选区或市政区等级。虽然这些文件类型不包括人口统计或经济信息，但您学会了如何发现数据资源并在笔记本中查看它们。将 JSON 文件保存为*.csv*格式还意味着您可以将其上传到如 QGIS 等 GUI 进行进一步探索。
总结
GeoPandas 是一个功能强大的开源 Python 项目，允许用户将地理数据与 pandas 对象集成，并进行空间操作。它在地理空间计算中默默地发挥作用。在本章中，您已经学习了文件的读取和写入、数据的选择、地图和图表的制作、交互式地图制作以及几个有用的包，以便将其包含在您未来的地理空间分析中。
^(1) 截至本文撰写时，2020 年人口普查的最终数据尚未公布。一旦数据公布，您只需在代码单元格中更新链接，即可获得最新信息。
第八章：数据清洗
处理数据时的一个通用问题是了解数据的完整性。数据工程依赖于清理、处理和可视化数据的能力。现在你已经熟悉了将数据与基于笔记本的代码编辑器集成的基本功能，无论是在本地使用 Jupyter Notebook 还是在云端使用 Google Colab，现在是学习如何清理你的数据的时候了。数据经常是不完整的（缺失），格式不一致，或者其他方面不准确——通常称为杂乱的数据。数据清洗是解决这些问题并准备数据进行分析的过程。
在这一章中，我们将探索一些公开可用的数据集，使用一些可以加载到 Colab 笔记本中的包来找到并清理杂乱的数据。你将使用来自纽约市开放数据门户NYC Open Data，更新于 2021 年 7 月 7 日的 NYPD_Complaint_Data_Historic 数据集。我筛选了 2020 年的数据，以便更容易查看和操作。你可以根据你的数据问题筛选数据并导出为 CSV 文件。本章将向你展示如何管理、删除、更新和 consoldate 数据，并用一些有用的 Python 包来处理数据。数据分析的准确性取决于数据集或数据库的质量，这一章将提供工具来评估和解决常见的不一致性问题。
检查缺失数据
如果你曾经参加过数据竞赛，比如在Kaggle上提供的那些竞赛，你可能会注意到数据集被设计成帮助你集中注意力解决特定任务，比如构建可视化。通常情况下，数据已经被清理过了。现实生活要复杂一些，缺失数据是一个持续存在的问题。让我们看看你如何让你的数据世界变得更加整洁。
使用电子表格和表格数据格式，评估数据的形状是一项直截了当的任务。稍微滚动一下就可以轻松地发现缺少数据的行和/或列。当你审查数据集以查找缺失数据的模式时，实际上是在评估它们的空值，或者说空数据（缺失值）的存在。地理空间分析，像一般的分析一样，通常涉及多个表格，因此学会如何识别这些表格数据之间的数据模式非常重要。
上传至 Colab
我在这个例子中使用了Google Colab，但你可以使用任何你喜欢的环境。将NYPD Complaint Data Historic CSV 文件上传到你的笔记本上，并点击文件名旁边的点。这将为你提供包含在笔记本中的路径（参见图 8-1）。
图 8-1 上传文件至 Colab
Missingno是一个 Python 库，用于检测数据集中的缺失值，并可视化它们在数据框架中的分布，这样可以更容易地发现模式。^(1) Missingno 在底层运行时依赖 NumPy、pandas、SciPy、matplotlib 和 seaborn，这使得底层代码片段更加熟悉。安装 missingno：
pip install missingno
根据数据集大小的不同，您可以调整样本大小。这里我设置了样本为 1000。导入 Python 库 pandas 将帮助您进行一些数据清洗和准备工作：
import pandas as pd
NYPD = pd.read_csv("/content/NYPD_Complaint_Data_Historic.csv")
import missingno as msno
%matplotlib inline
msno.matrix(NYPD.sample(1000))
一旦您上传数据集，您将需要查看数据字典，了解列标题表示的内容。
我们数据框架中列出的列显示在表 8-1 中。CMPLNT_NUM 似乎是每行的关键标识符。这一点很重要，因为每行都是数据集中一个事件的记录。
表 8-1. 纽约市警察局投诉数据集数据字典
| Field name | 描述 |
|---|---|
| CMPLNT_NUM | 每个投诉的随机生成的持久标识符 |
| CMPLNT_FR_DT | 报告事件发生的确切日期（或者如果存在 CMPLNT_TO_DT，则为事件开始日期） |
| CMPLNT_FR_TM | 报告事件发生的确切时间（或者如果存在 CMPLNT_TO_TM，则为事件开始时间） |
| CMPLNT_TO_DT | 如果确切时间不明确，则报告事件结束日期 |
| CMPLNT_TO_TM | 如果确切时间不明确，则报告事件结束时间 |
| ADDR_PCT_CD | 发生事件的警区 |
| RPT_DT | 向警察报告事件的日期 |
| KY_CD | 三位数的犯罪分类代码 |
| OFNS_DESC | 与关键代码对应的犯罪描述 |
| PD_CD | 三位数的内部分类代码（比关键代码更详细） |
| PD_DESC | 与 PD 代码对应的内部分类描述（比犯罪描述更详细） |
| CRM_ATPT_CPTD_CD | 犯罪行为是否成功完成、尝试但失败、或过早中断的指示符 |
| LAW_CAT_CD | 犯罪级别：重罪、轻罪或违规行为 |
| BORO_NM | 事件发生的市区名称 |
| LOC_OF_OCCUR_DESC | 发生地点具体描述：内部、相对于、前部或后部 |
| PREM_TYP_DESC | 场所的具体描述：杂货店、住宅、街道等 |
| JURIS_DESC | 管辖区代码的描述 |
| JURISDICTION_CODE | 事件责任管辖区：可以是内部机构，如警察（0）、运输（1）和住房（2），也可以是外部（3），如矫正局、港务局等。 |
| PARKS_NM | 如果适用，纽约市公园、游乐场或绿地发生事件的名称（不包括州立公园） |
| HADEVELOPT | 如果适用，纽约市住房管理局发生事件的住房开发项目名称 |
| HOUSING_PSA | 发展级别代码 |
| X_COORD_CD | 纽约州平面坐标系统的x坐标，长岛区，NAD 83，单位为英尺（FIPS 3104） |
| Y_COORD_CD | 纽约州平面坐标系统的y坐标，长岛区，NAD 83，单位为英尺（FIPS 3104） |
| SUSP_AGE_GROUP | 嫌疑人年龄组 |
| SUSP_RACE | 嫌疑人种族描述 |
| SUSP_SEX | 嫌疑人性别描述 |
| TRANSIT_DISTRICT | 事件发生的运输区 |
| Latitude | 全球坐标系统的中部纬度坐标，WGS 1984，十进制度数（EPSG 4326） |
| Longitude | 全球坐标系统的中部经度坐标，WGS 1984，十进制度数（EPSG 4326） |
| Lat_Lon | 地理空间位置点（纬度和经度的组合） |
| PATROL_BORO | 事件发生的巡逻区名称 |
| STATION_NAME | 运输站名称 |
| VIC_AGE_GROUP | 受害人年龄组 |
| VIC_RACE | 受害人种族描述 |
| VIC_SEX | 受害人性别描述 |
数据集中的变量需要确切地按照数据框中呈现的方式命名。如果它们全部大写或用下划线分隔，您需要确保复制它们（或者可以类似地重命名它们，就像您在第七章中重命名人口普查数据列一样）。
接下来，还有其他需要在数据集中考虑的问题。
空值和非空值
您首先需要检查数据集中的缺失值。通过使用像NYPD.info()和NYPD.describe()这样的 pandas 功能快速查看数据（参见图表 8-2 和 8-4）。
使用.info()可以显示列名、非空值数和数据类型（Dtype）：
NYPD.info()
在这个数据集中，范围索引有 1115 条记录；因此，任何具有少于 1115 个值的列都存在缺失数据。列的非空计数可以帮助您查看缺失的数据，并决定哪些列未捕捉到足够的数据以提供洞察。
该数据集的文档指出，空值可能归因于部门表格的变更和数据收集的不一致性。此外，如果 2020 年（我选择用来过滤数据的年份）的信息不可用或在收集时不明确，则将其分类为未知/不可用/未报告。重要的是要阅读支持文档，以了解数据的限制，这可能会限制您能够有用地提出的问题。
数据类型
数据清洗中的首要任务之一应该是检查数据集中每一列的数据类型。最常见的数据类型包括布尔值、整数、浮点数、对象和日期时间，您应该对这些类型非常熟悉。
这些数据类型中最棘手的是对象。Python 对象可以包含各种数据类型。它们通常是字符串，但您应该知道它们还可以包括整数、浮点数、列表和字典。在 pandas 库中，不同的数据类型被检测并分类为 NaN（默认标记为“不是数字”）或 NaT 表示缺少日期时间数据。NYPD.info的输出包括数据类型（Dtype），或者您可以使用dtypes DataFrame 属性，通过在其后添加点号来使用NYPD.dtypes。
元数据
您还可以返回有关数据形状的信息。这与数据类型信息一起被称为元数据，或者称为数据的数据。要查看 NYPD 数据的形状，请运行：
NYPD.shape
(7375993, 35)
简化起见，在此示例中，我在源数据中筛选了仅包含 2020 年数据，以便我们处理较少的记录：1,115 条记录（图 8-2）。
图 8-2. NYPD.info()的数据框摘要
每次投诉发生的确切时间，或者 CMPLNT_FR_TM，在图 8-2 中列出为对象或字符串。当 datetime 或其他数据类型被分配为字符串或对象时，Pandas 通常会猜测列中的数据类型。这通常会由于low_memory=True等默认设置而变得复杂。大文件会被“分块”，并且该文件通常会生成包含多个数据类型的列。我建议您不要简单地将参数切换为False。手动更改dtype更为高效，如下面的代码单元格所示：
NYPD = pd.read_csv('/content/drive/MyDrive/NYPD_Complaint_Data_Historic.csv',
parse_dates=['CMPLNT_FR_TM'])
NYPD.dtypes.head()
CMPLNT_NUM int64
CMPLNT_FR_DT object
CMPLNT_FR_TM datetime64[ns]
NYPD = NYPD.drop(columns =['CMPLNT_TO_DT','CMPLNT_FR_TM','PARKS_NM','HADEVELOPT', 'HOUSING_PSA'])
NYPD.head()将更新您保留的列，并确认您的“删除”列不再显示在数据帧中。
如果不与您感兴趣的探索相关，您可以删除具有不完整数据的列。但我强烈建议您不要删除数据。通常，数据缺失的原因和选择包含的数据一样引人注目或重要。
让我们通过识别几个样本列来尝试这个功能从数据帧中删除。
摘要统计
您可以使用.describe()查看摘要统计，或者汇总您的样本数据的统计信息。
查看行计数是识别数据帧中不同特征缺失数据的另一个机会。在df.describe()的输出中，排除了 NaN 值，因此您可以看到数据的实际分布以及每列中值的计数。
中心趋势测量
正态分布是指平均值、众数和中位数都是相同的值。概率分布或中心趋势测量描述了分布中的中心值，当分布不正常时，如图 8-3 所示。
Figure 8-3. 中心趋势测量的回顾
还有一些指标可以表明数据分布的偏斜程度和异常值的存在。 图 8-4 是2021 年能源和水数据披露法律 84 号（2020 年日历年数据）的一个快照，这是一份详细考察纽约市私人建筑物能源和水使用情况的大型数据集。表中的中位数或第二四分位数（Q2），即表中的 50%指标，表示中位数值，即一半数值在其上方，一半数值在其下方。
Figure 8-4. df.describe() 函数提供了数值数据列的描述性统计信息。
快速查看数据形状的一种方法是观察 Q2 或中位数的位置。如果查看 图 8-4 中的第二列，您可以确定中位数的值更接近 25%还是 75%。如果更接近 25%或 Q1，则数据右偏，意味着有 25%的值低于 Q1 值。将 Q1 视为最小值和中位数之间的值。Q3 的值位于中位数和最大值之间，指示 75%的报告值将低于 Q3。这将导致数据具有左偏。
您可以通过 df.describe() 函数快速查看数据集中的数值值。将以下代码段输入到您的笔记本中以查看：
df.describe()
绘制直方图有助于更好地可视化数据的形状。如 图 8-5 所示，存在右偏尾或正偏态。这表明均值大于中位数，中位数大于众数。
Figure 8-5. 数据分布的可视化
让我们回到 NYPD 数据。我们可以可视化数据并查看其分布：
size, scale = 1000, 10
Complaints = pd.Series(np.random.gamma(scale, size=size) ** 1.5)
Complaints.plot.hist(grid=True, bins=20, rwidth=0.9,
color='#607c8e')
plt.title('Distribution of Complaint Types')
plt.xlabel('CMPLNT_NUM')
plt.ylabel('OFNS_DESC')
plt.grid(axis='y', alpha=0.75)size, scale = 1000, 10
Complaints = pd.Series(np.random.gamma(scale, size=size) ** 1.5)
Complaints.plot.hist(grid=True, bins=20, rwidth=0.9,
color='#607c8e')
plt.title('Distribution of Complaint Types')
plt.xlabel('PATROL_BORO')
plt.ylabel('OFNS_DESC')
plt.grid(axis='y', alpha=0.75)
可视化还显示了投诉类型数据的右偏分布（见 图 8-5）。
图 8-6 简单统计了数据集中每一列缺失的数值。在分析数据之前了解记录缺失的时机提供了透明度，并且在分享数据可视化时，传达这些信息至关重要。
Figure 8-6. NYPD.isna().sum() 列出了数据集中每一列的缺失值统计。
NYPD.isna() 返回带有布尔值指示缺失值的数据框，如 图 8-7 所示。在计算中，布尔值只有两个可能的值：True 和 False，通常表示为 1 和 0。在这种情况下，True 表示缺失数据。一目了然，您可以看到数据缺失的位置。
图 8-7. 数据框 NYPD.isna()
这些只是探索数据框并了解数据完整性或空值范围的几个选项。
现在您已经找到了缺失数据，接下来该怎么办？
替换缺失值
检查数据框的另一个原因是查看缺失值并确定它们的特征。一旦您知道表中如何记录缺失值，您就可以使用 na_values 参数将它们替换为您选择的变量，例如 NaN、unknown 或其他感兴趣的值：
df_test = pd.read_csv("/content/NYPD_Complaint_Data_Historic.csv",
na_values = 'NaN')
print(df)
探索数据集通常会发现数据表中缺失值报告的不一致性。代码片段将查询替代变量，如 ? 和 999。
使用 Missingno 可视化数据
在本节中，我们将再次使用 Missingno。Missingno 可以将表格数据矩阵转换为布尔掩码，因为 掩码 基于底层操作数据的某些标准返回数据。数据根据指定标准返回 True/False 布尔值。在计算中，布尔值只有两个可能的值。Missingno 将包含数据的单元格标记为 True，空单元格标记为 False。它可以将这些数据可视化为 空值矩阵（图 8-8），这种图表使用空白空间来直观表示缺失数据并显示模式。数据存在的地方将被阴影化。
图 8-8. 纽约警察局投诉数据集的空值矩阵
在 图 8-8 中，每列代表数据中的一个特征，每行代表一个观察结果。矩阵右侧的小火花线显示了每个记录中缺失数据的量。
您在这里看到了任何缺失数据模式吗？我注意到第四和第五栏，即投诉结束日期和时间，具有相同的模式：当这两个值中的一个缺失时，另一个也缺失。您可以使用 msno.bar 将相同信息显示为列，而不是简化的矩阵。
msno.bar(NYPD.sample(1000))
柱的高度等于空值或缺失数据的级别（图 8-9）。较高的柱表示值完整（没有缺失数据）或几乎完整。
图 8-9. 通过列显示的空值
评估 msno.heatmap(NYPD) 可以检查 空值相关性，即变量的存在或缺失以及这如何强烈地影响另一个变量的存在。相关性的度量范围为 –1 到 1，如图 8-10 所示。零相关 意味着一个变量的存在似乎与另一个变量的存在无关；相关分数为 1 表示一个变量在另一个存在时也存在。负值表示 负相关：一个变量的存在意味着另一个变量不会存在。
图 8-10. 空值数据热图（NYPD）
始终存在或从不存在的值不会在热图中可视化。例如，有两列——RPT_DT（报告给警察的日期）和 JURIS_DESC（辖区描述）——报告没有缺失数据，它们不会在热图中可视化。其他似乎不测量任何内容的值可能是由错误数据生成的。对数据的更详细检查可能会发现一些有趣的东西。请注意小于 1 的值：例如，图 8-10 中的 –0.3 值值得更仔细地观察。例如，小于 –1 的值表示相关性几乎是 100% 的负相关。
树状图，或称树状图，是一种分支图，用于直观地表示不同组之间的关系。分支称为 类群，末端称为 叶子。类群的排列告诉我们哪些叶子最相似，基于它们的接近程度。分支点的高度指示它们之间的相似性或差异性：高度越大，差异性越大。树状图中分组的越接近 0，表示一列中的空值出现与其他列的空值出现或缺失有更密切的关系。
scipy.cluster.hierarchy.dendrogram 通过绘制连接类群的 U 形链接来组成聚类。如果列数少于 50，则树状图的方向默认为自上而下；如果列数超过 50，则默认为左右方向。
在图 8-11 中，表示嫌疑人种族和性别的列（SUSP_RACE 和 SUSP_SEX）比表示嫌疑人年龄组和转运区（SUSP_AGE_GROUP 和 TRANSIT_DISTRICT）的列更相似。
图 8-11. NYPD 数据集中的树状图关系
要查看层次聚类，请将以下代码片段写入代码单元格中：
msno.dendrogram(NYPD)
在关注数据集内数据分布时，树形图是非常有用的。层次聚类定量估计每个数据点与数据集中的每个其他点的关系。聚类距离在垂直轴上，不同的数据点水平显示。数据样本的最高粒度在数据样本水平上，树形图是一种快速可视化的方法。
映射模式
QGIS 是查看数据和查找数据模式的另一种方式。图 8-12 中的数据被筛选，仅显示入室盗窃案，以帮助根据位置指导额外洞察。让我们看看如何生成这样的地图。
图 8-12. 在地图上可视化 NYPD 筛选后的犯罪数据
纬度和经度
数据集有纬度和经度列。你可以将它们转换为指出特征的点，然后使用 GeoPandas 创建地图。让我们使用相同的文件导入数据并添加一个 shapefile。首先，导入 pandas：
import pandas as pd
接下来，read_csv 将数据文件读入你的笔记本：
df = pd.read_csv('/content/drive/MyDrive/NYPD_Complaint_Data_Historic.csv')
如果你使用的是 Google Colab，可以直接上传文件。你可以将 df 变量修改为任何你喜欢的内容。df.head() 函数将返回每列的前几行：
df.head()
作为最佳实践，请确保精确捕获数据中出现的拼写和大写，编写代码时始终将经度列在纬度列之前。
Shapefiles
你将使用 Python 库 Shapely 进行计算几何，或处理诸如点、曲线和曲面等几何对象。上传 shapefile 时，请确保将与 shapefile 相关的所有文件（它们将具有 .dbf、.shx 和 .prj 扩展名）上传到同一目录中。图 8-13 中的 shapefile 来自NYC 开放数据区界。
图 8-13. NYC 开放数据区界
首先，导入你的库和包：
!pip install geopandas
import geopandas as gpd
from shapely.geometry import Point, Polygon
import matplotlib.pyplot as plt
street_map = gpd.read_file('/content/tl_2021_36_place.shp')
提示
Python 库在代码中通常被缩写。你已经看到 pandas 被缩写为 pd，而 GeoPandas 被缩写为 gpd。请注意，plt 用于 matplotlib —— 不要误以为它是函数或参数的一部分。
现在 Python 只需要知道如何在定义的空间中应用坐标 —— 所以你需要定义它：
# designate coordinate system
crs = {'init': 'epsg:4326'}
# zip x and y coordinates into single feature
geometry = [Point(xy) for xy in zip(df['Longitude'], df['Latitude'])]
# create GeoPandas dataframe
geo_df = gpd.GeoDataFrame(df,
crs = crs,
geometry = geometry)
GeoPandas 是笛卡尔坐标参考系统，这意味着每个点由一对数字坐标定义，比如我们示例中的纬度和经度。它将地理数据分配给 pandas 对象。GeoPandas 数据框从定义我们的 CRS 和组合我们的经度和纬度列而创建。Shapely 中的 Point 函数使用经度和纬度列创建一个名为几何的新列，如图 8-14 所示。
图 8-14. 添加新列，几何图形
在前述代码单元格中，我们定义了geo_df，现在我们可以查看新创建的几何列：
geo_df.head()
这里我扩大了figsize以便更好地查看，但在较小的尺寸下加载速度会更快。默认单位是英寸，但也可以进行厘米和像素的转换。street_map被分配给轴，因为这是我们分配形状文件的位置：
# create figure and define axes, assign to subplot (matplotlib) fig, ax = plt.subplots(figsize=(30,30))
# add .shp mapfile to axes street_map.plot(ax=ax, alpha=0.4,color='grey')
# add geodataframe to axes
# assign 'OFNS_DESC' variable to represent coordinates on graph
# add legend
# make data points transparent using alpha
# assign size of points using markersize geo_df.plot(column='OFNS_DESC',ax=ax,alpha=0.5, legend=**True**,markersize=10)
# add title to graph plt.title('Reported Offenses', fontsize=15,fontweight='bold')
# set Latitude and Longitude boundaries for map display plt.xlim(-74.02,-73.925)
plt.ylim( 40.7,40.8)
# show map plt.show()
图 8-15 是从 GeoPandas 数据框创建的地图。您可以过滤和编辑报告的事件，以创建更具针对性的地图。（参考Jupyter Notebook获取更多选项。）plt.xlim 和 plt.ylim 命令可让您选择特定的边界来进一步编辑您的投影。
如果您想从数据集中选择一种犯罪类型，请使用df.loc函数来定位所有实例。以下是一个显示入室盗窃的示例：
fig, ax = plt.subplots(figsize=(15,15))
street_map.plot(ax=ax, alpha=0.4,color='grey',legend=True,markersize=20)
geo_df.loc[df['OFNS_DESC'] == 'BURGLARY']
geo_df.plot(ax=ax, alpha=0.5,)
plt.xlim(-74.02,-73.925)
plt.ylim( 40.7,40.8)
plt.show()
或者您可能希望列出几种不同的违法行为：
df.loc[df['OFNS_DESC'].isin([value1, value2, value3, ...])]
或者一组不同的参数：
df.loc[(df['OFNS_DESC'] == value) & (df['TRANSIT_DISTRICT'] == value)]
图 8-15. GeoPandas 在纽约市各区地图上可视化的数据框
尝试提出不同的问题！
摘要
您已学习了一系列在清理数据时非常有用的方法。这一点尤为重要，特别是在处理开放源码数据时。例如，社区级别的数据通常是手动输入的，拼写错误、缺少日期和几何变量的缺失可能会限制这些宝贵资源的实用性。
^(1) Bilogur, Aleksey. 2018. “Missingno: A Missing Data Visualization Suite.” Journal of Open Source Software 3 (22): 547. https://doi.org/10.21105/joss.00547
第九章：探索地理空间数据抽象库（GDAL）
讨论开源平台和库，如果没有介绍Geospatial Data Abstraction Library (GDAL)，将是不完整的。这是一个有效处理栅格和矢量地理空间数据的资源。处理栅格和矢量数据需要一系列工具，而 GDAL 在本书中使用的许多程序的内部工作中起着关键作用，包括 ArcGIS、QGIS 和 Google Earth。当依赖图形用户界面的程序化便利时，很容易忘记将多种数据类型和格式结合起来，以便以统一的数据模型高效地工作。这就是为什么 GDAL 是一个重要工具。它简化了在各种格式和空间参考系统中处理地理空间数据的工作。
在本章中，我们将处理栅格数据，并看看如何使用 Spyder IDE 处理 GDAL，包括如何使用warp函数更改地图投影，处理栅格波段，转换文件并创建二进制掩模。我们还将快速介绍另外三个有用的数据集资源：EarthExplorer、Copernicus 开放访问中心和 Google Earth Engine（GEE）。
首先，我会向你展示如何使用命令行界面，也称为终端，快速读取、转换和重新投影你的地理空间数据。为什么要使用命令行？如果你正在处理多个需要相同功能的文件，你不必手动逐个操作，试图在 QGIS 或 ArcGIS 中回忆必要的步骤和过程。使用命令行，你可以将一行简单的代码保存为脚本，用来处理各种文件和数据集。正如本章稍后将看到的那样，这些脚本也可以在 QGIS 中运行。在命令行工作时，我建议在 QGIS 中查看结果，或者在你对终端工作感到舒适后，可以参考 Spyder IDE 部分。
其次，我会向你展示如何在 IDE 中使用 GDAL，这是一种集成了编程和应用开发的工具。许多初学者 Python 课程都是在 IDE 中教授的，所以如果你有一些 Python 编程背景，可能对至少一个 IDE 比较熟悉。Python 有很多 IDE，但在本章中，我们将使用一个叫做 Spyder 的 IDE。像 Spyder 这样的 IDE 通过实时代码内省（允许你直接查看函数和类是如何执行以及它们包含的信息）、与 matplotlib 一起显示的内联图形以及我最喜欢的变量资源管理器（提供了编辑和可视化代码中对象的内置支持）等功能，简化了与代码的工作和编写。不管你使用哪个操作系统，IDE 看起来都是一样的，这在你的最终用户访问不同的工作流程和资源时非常方便。
配置 GDAL
从历史上看，这个库被称为GDAL/OGR，其中 GDAL 是栅格库的一部分，而 OpenGIS 简单要素参考实现（OGR）是库的矢量部分。现在，这些组合库通常被称为GDAL。我知道有许多数据科学家几乎完全在命令行上工作，所以我想快速演示一下如何实现这一点。对我来说，在命令行工作使我更容易熟悉调用 GDAL 函数的语法。它的类和函数与你从 Python 中了解到的有些不同，尽管在使用 Python API 时，语法会很熟悉。GDAL 保持对 Python 的绑定，但这与你迄今为止使用 Python 库的方式有所不同。即使使用 Python 绑定，你也需要知道如何调用 GDAL 函数并理解它们的功能。^(1) 在使用 IDE 工作的部分中，我们将回到 Python。
安装 Spyder
Spyder IDE 完全由 Python 编写。它以被数据科学家、数据分析师和数据工程师使用而闻名。其功能包括自动补全和实时分析。让我们首先安装它。
我建议为 Spyder 安装创建一个新环境。也可以将 GDAL 安装到现有环境中；你可以输入**conda env list**查看已有的环境。你还需要添加 NumPy 和 matplotlib。
在相同的环境和目录中安装 Spyder 和 GDAL，用你自己的文件路径替换我的：
spyder-env/Users/bonnymcclain/opt/miniconda3/envs/spyder-env
conda install -c conda-forge spyder
安装完 Spyder 后，可以通过在终端中输入**spyder**来启动它。
安装 GDAL
你还将在命令行安装 GDAL。打开终端并在提示符下运行以下命令：
Conda install -c conda-forge gdal
要检查安装是否成功，请运行：
gdalinfo —-version
可能需要设置安装文件夹的路径以从终端访问它。接下来将介绍如何完成这些步骤。我还建议查阅GDAL 文档。
在命令行中使用 GDAL
如今，普通计算机用户主要在 GUI 中工作：这是一个图形层，使得可视化文件和执行任务变得更加容易。但命令行也有其优点，包括更快的处理速度。不需要加载大量文件，如果想在不同数据集上重复分析，可以保存小型 shell 脚本执行重复操作。
激活你首选的环境，conda activate minimal_ds，或者你刚创建的用于处理 GDAL 的 Spyder 环境。我想要在我的 TIFF 文件夹中操作文件，所以我已经输入cd来将其设置为当前目录。现在你可以看到 TIFF 文件被添加为我的工作目录：
(minimal_ds) MacBook-Pro-8:~ bonnymcclain$ cd TIFFs
(minimal_ds) MacBook-Pro-8:TIFFs bonnymcclain$ ls
要在不同目录之间移动，请使用ls/深入了解您的目录结构。例如，如果您想要处理 Downloads 文件夹（假设它不是您当前的目录），请使用以下代码来浏览 Downloads 目录中的文件夹：
ls/Downloads/another folder
(minimal_ds) MacBook-Pro-8:~ bonnymcclain$ cd TIFFs
(minimal_ds) MacBook-Pro-8:TIFFs bonnymcclain$ ls
花点时间在不同目录之间移动，因为理解结构是在终端中使用 GDAL 包进行工作时最重要的一课。再次强调，访问 Spyder IDE 中的包、库和文件时，您需要理解结构。
如果您在正确的目录中，冒号后面和您的用户名之前将会看到目录名。这将输出该目录中的文件列表。
TIFFs bonnymcclain$ ls
我已经包括了一个关于使用GEE创建您自己文件资源的部分，以便您可以跟进。Copernicus 公开访问中心和EarthExplorer也是获取栅格数据的来源，但由于较大文件的大小通常渲染速度较慢。
现在您的安装已经完成，您可以学习一些命令来帮助您入门。这里是GDAL 文档，可以帮助您解决下载和安装中的任何问题。
使用 GDAL 编辑您的数据
gdalinfo命令行参数展示了关于您的栅格数据集的可用信息。您可以使用gdalinfo命令并提供活动目录中文件的名称来查看这些参数。在终端中输入以下代码：
gdalinfo [--help-general] [-json] [-mm] [-stats | -approx_stats] [-hist] [-nogcp]
[-nomd]
[-norat] [-noct] [-nofl] [-checksum] [-proj4]
[-listmdd] [-mdd domain|`all`]* [-wkt_format WKT1|WKT2|...]
[-sd subdataset] [-oo NAME=VALUE]* [-if format]* datasetname
在这个示例中，我从 GEE 数据片段中选择了一个保存的文件。您可以在文档中了解更多关于命令行参数的信息。
输入以下代码，用您目录中*.tif*文件的名称替换（注意不需要括号）：
(minimal_ds) MacBook-Pro-8:TIFFs bonnymcclain$ gdalinfo Sentinel2_RGB20200501.tif
这里是输出的一部分片段：
(minimal_ds) MacBook-Pro-8:TIFFs bonnymcclain$ gdalinfo Sentinel2_RGB20200501.tif
Driver: GTiff/GeoTIFF This is the format of the saved file
Files: Sentinel2_RGB20200501.tif
Size is 5579, 4151
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 29N",
BASEGEOGCRS["WGS 84",
…
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (566320.000, 4133200.000) (8d15' 4.53"W, 37d20'35.24"N)
Lower Left (566320.000, 4091690.000) (8d15'17.78"W, 36d58' 8.31"N)
Upper Right (622110.000, 4133200.000) (7d37'17.50"W, 37d20'14.86"N)
Lower Right (622110.000, 4091690.000) (7d37'41.89"W, 36d57'48.20"N)
Center (594215.000, 4112445.000) (7d56'20.39"W, 37d 9'13.16"N)
您可以从第二行看出保存文件的格式是 GTiff/GeoTIFF。
小贴士
要查看其他格式的列表，请在命令行中输入gdal_translate --formats。您还可以查看文件大小、坐标系、像素大小和坐标。
这里是显示关于颜色波段信息的输出部分：
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
Description = B2
Band 2 Block=256x256 Type=Float32, ColorInterp=Undefined
Description = B3
Band 3 Block=256x256 Type=Float32, ColorInterp=Undefined
Description = B4
注意 Sentinel 卫星数据波段的颜色解释。它们被设置为Gray用于 B2 波段，对于剩余的两个波段则是ColorInterp=Undefined。由于未下载元数据，您需要帮助 GDAL 解释这些波段。
由于这是Sentinel 2 数据，您知道波段是蓝色（B2）、绿色（B3）和红色（B4）。您将直接在原地编辑数据集，使用命令gdal_edit.py，接着是选项（如colorinterp），最后是值（这里是波段的颜色值）。最后一步是提供输入文件，Sentinel2_RGB20200501.tif：
gdal_edit.py -colorinterp_1 blue -colorinterp_2 green -colorinterp_3 red
Sentinel2_RGB20200501.tif
如果您的功能需要输出文件，则还需要添加output.tif。这里只需要输入文件。
当您重新运行 gdalinfo Sentinel2_RGB20200501.tif 时，您会看到颜色波段已更新，现在显示为蓝色、绿色和红色。
变形功能
您还可以通过学习使用 gdal_warp 函数来更改栅格投影。（一旦熟悉 warp 函数，您可以用它来学习其他常见函数；要探索它们，建议从 Python API 的 GDAL 文档开始。）
在以下代码中，参数 -t_srs 指定您要定位的坐标系统。例如，每个地理坐标系统都分配了唯一的 EPSG 代码。*xxxxx* 的位置，请输入您希望更改为的 EPSG 代码。输入是栅格文件，输出是重命名的修改文件，如 output_rd.tif 所示：
gdalwarp -t_srs EPSG: xxxxx Sentinel2_RGB20200501.tif output_rd.tif
gdalinfo Sentinel2_RGB20200501.tif
您的输出将展示新的投影。
捕获输入栅格波段
接下来，您将使用 GDAL 命令行捕获投影波段的图像统计信息：
(minimal_ds) MacBook-Pro-8:TIFFs bonnymcclain$ gdalinfo -stats srtm_41_19_4326.tif
输出摘录：
STATISTICS_MAXIMUM=640
STATISTICS_MEAN=256.70790240722
STATISTICS_MINIMUM=27
STATISTICS_STDDEV=119.8746675927
STATISTICS_VALID_PERCENT=100
为什么要捕获输入栅格波段进行统计分析？主要用于分类：在查看不同位置时识别群集。例如，波段特征有助于确定您是否观察到植被、水体或者住宅区。
我在本章前面提到，GDAL 在 QGIS 中以及多种其他程序中运行。 图 9-1 展示了代码在 QGIS 平台上的运行，上传您的 .tif 文件，并选择栅格信息。
图 9-1. 在 QGIS 中显示 gdalinfo
现在您已经体验了在命令行中使用 GDAL，我鼓励您进一步探索更多可用的任务。这些技能对处理和编辑地理空间数据和应用程序的核心算法非常有用。在处理各种矢量和栅格数据格式时，跨平台和多样化格式转换的能力尤为重要。
接下来，我们将与 Python 集成的 GDAL 进行操作。
使用 Python 处理 GDAL 库
当您启动 Spyder（再次在终端中键入 **spyder**）时，Spyder 控制台将在新的浏览器窗口中打开，如 图 9-2 所示，但是会是空白的。
图 9-2. Spyder 控制台
在 Spyder 中的初始定位
图 9-2 的左侧是 Spyder 的脚本编辑器，您可以在其中创建和打开文件。（您可以自由排列面板和控制台。）右侧的两个控制台是您可以探索生成的变量和定位文件的地方。显示图像的面板是绘图窗格。您可以在底部控制台中编写简单的代码。如需帮助，请选择上窗口的帮助选项卡。
您可以在 Python 控制台或编辑器中运行代码。您还可以使用 Python 控制台控制 Spyder 的调试器。您生成的图像和图表将显示在绘图面板中或嵌入到控制台中。
变量资源管理器显示在图 9-2 右上窗格中，是我在集成开发环境中喜欢工作的一个原因。这些对象是在运行代码时生成的。点击一个变量以深入探索其详细信息。
每个控制台窗口还有一个“汉堡包”图标，或☰，展开成一个菜单（图 9-3），您可以在其中找到有关窗格的其他信息。
图 9-3. 在控制台目录中工作
这些选项是可定制的。我建议访问Spyder 网站获取界面详细信息以及探索所有可用窗格。如果您想在安装之前尝试 Spyder，还有一个Binder 选项，可以在浏览器中运行。
虽然这并非必须，但我喜欢在 Spyder 中启动一个新项目，就像在图 9-4 中演示的那样。这简化了返回到之前文件的过程，而不管当前工作目录如何。项目窗格也将变得可见。这让我想起 QGIS 的便利，因为文件很容易到手。
图 9-4. 在 Spyder 中创建新项目
项目窗格还可以与 Git 集成，用于版本控制。查阅Spyder 文件文档以了解如何在您的存储库中设置此功能。
在 Spyder 中探索您的数据
您距离在集成开发环境中探索您的栅格数据只有几行代码。第一步是导入必要的包：GDAL、NumPY 和 matplotlib：
from osgeo import gdal
import numpy as np #to manipulate the array from raster
import matplotlib.pyplot as plt #to view figure
NumPY 将操作从栅格文件创建的数组，而 matplotlib 将允许您查看文件。您可以探索的栅格程序包含在GDAL 文档中。
实用脚本属于实用类，这是一个跨应用程序可用的相关方法集合。Python 实用脚本位于osgeo_utils.samples子模块中，如下面的代码所示：
ds = gdal.Open("slope.tif")
gt = ds.GetGeoTransform()
proj =ds.GetProjection()
band = ds.GetRasterBand(1)
array = band.ReadAsArray()
如果您不确定如何找到*.tif*文件，请跳转到“探索开源栅格文件”，我将介绍如何查找栅格文件以进行探索。这些是简单的数字高程模型（DEMs），其中每个点或像素都具有高程值。通常，它们表示为 DEM .tif文件。如果您在阅读本书时一直在探索，那么您的下载文件夹可能已经拥有不少文件了。
你将写入编辑器的代码脚本在文本中有解释；更多信息可从GDAL Python API获取。
你可以随意给你的变量取任何名字，但为了简单起见，我使用 ds 来表示数据集。导入你的 .tif 文件并创建变量 ds 后，你会看到变量出现在变量资源管理器中。ds 文件的格式是 gdal.Dataset。
GDAL 中的文件转换
你需要将你的 .tif 文件从图 9-6 中所见的表格格式转换为地理坐标，因此接下来你将定义 gt，地理转换。
图中图 9-5 中的六个系数，从上到下依次对应：
0
左上角像素的x坐标（262846.525725）
1
像素分辨率，从西到东（25.0）
2
行旋转（通常为 0）
3
左上角的像素的y坐标（4464275.0）
4
列旋转（再次，通常为 0）
5
北西像素分辨率和高度（通常为负值表示北向图像），-25.0
图 9-5. 地理转换为地理参考坐标
投影信息在变量资源管理器中可见（图 9-2），但在图 9-6 中显示更大。在通用横轴墨卡托投影（UTM）中，投影是 UTM 30 北区；EPSG 代码是 EPSG:32630。
图 9-6. GDAL 读取的波段的投影
GetRasterBand 将波段提取到数据集中。要确定你有多少个波段，写入控制台：
ds.RasterCount
这会输出：1。
输入波段数到函数中：
band = ds.GetRasterBand(1)
array = band.ReadAsArray()
探索数组变量，你会看到 Python 已经将 GDAL 文件读入 NumPY 数组。数组共享我们 DEM 中的高程，如图 9-7 所示。
图 9-7. DEM .tif 文件的 NumPY 对象数组
在 GDAL 中使用二进制掩码
掩模 对于裁剪边界或值到特定范围非常有用。 二进制掩模（通常称为binmask) 是一个过滤掉高于或低于某个数字的高度值的工具（例如，如果你只想查看低于海平面的位置）。如果你想保存一切大于或等于平均值的内容，你将把它赋值为 1；否则，它将是 0. 输出显示在图 9-8 中。暗值处于较低高度。
代码中的 binmask 变量指的是基于条件（如 mean）返回数组中元素的 NumPy (np) 函数：
binmask = np.where((array >=np.mean(array)),1,0)
plt.figure()
plt.imshow(binmask)
要将这些数据保存为 GeoTIFF 文件，您将需要一个 GDAL 驱动程序来支持您选择的栅格文件格式。您可以查看不同驱动程序的长列表，包括 GeoTIFF，在GDAL 文档中找到。
binmask.shape
这输出为（410，601）。
图 9-8. 二进制掩模输出
回顾如何在 Python 中使用index函数。binmask 的形状为 410 行 601 列。在下面的代码片段中，xsize指的是列数。您通过索引调用列，在本例中为[1]，以获取计数（在本例中为 601）。您可以用ysize和索引[0]来获取行数（在本例中为 410）：
#Now we want to save as GeoTiff so we need a GDAL driver
driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tif", xsize = binmask.shape[1],
ysize = binmask.shape[0], bands =1,
eType=gdal.GDT_Int16)
gdalconst定义了图像中的数据类型。由于您正在使用0和1定义 binmask 的形状，您需要将值设置为整数，这将使用gdal.GDT_Int16完成。现在，您只需要打印数据类型属性。
最后一步是另一个地理转换。由于您没有更改任何内容，设置可以保持不变。一旦关闭文件，它们将可供您使用。如果您忘记了，您将无法使用它们。您也不需要读取数组，因为您还没有将其存储在任何地方。WriteArray(binmask)将提供输出：
#need geotransform to be complete
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()
#close
outband = None
outds = None
总结一下，您已经从 DEM 识别并生成了一幅图像，并将输出转换为 GeoTIFF 文件。栅格图像与关于图像在地球表面像素级位置的任何信息或元数据一起保存。
完整脚本
这是完整的代码；当您想查看图像时，请取消注释打印选项：
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
ds = gdal.Open("slope.tif")
gt = ds.GetGeoTransform()
proj=ds.GetProjection()
band = ds.GetRasterBand(1)
array = band.ReadAsArray()
#plt.figure()
#plt.imshow(array)
binmask = np.where((array >=np.mean(array)),1,0)
plt.figure()
plt.imshow(binmask)
driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tif", xsize = binmask.shape[1],
ysize = binmask.shape[0], bands =1,
eType=gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()
#DON’T FORGET TO CLOSE FILE
outband = None
outds = None
探索开源栅格文件
开源地理空间社区拥有大量的公开数据集，供您继续学习并支持当前和未来的项目。这本书旨在成为一个不断发展的资源，邀请您进行额外的学习和技能发展。因此，在接下来的几节中，我将提供一些基础知识，让您开始探索这些资源，而不是进行完整的练习。目标是让您能够开始自主探索。
USGS EarthExplorer
美国地质调查局（USGS）托管了最大的免费卫星和航空图像数据库之一，名为 EarthExplorer。
您需要注册一个免费帐户在EarthExplorer。如果您有一个包含 shapefile（.shp）的压缩 ZIP 文件，您可以使用左上角的 KML/Shapefile 上传按钮上传它。免费的开源数据资源通常提供下载压缩的 shapefile 的选项。必须上传 shapefile 的所有附带文件。
要使用 EarthExplorer，你可以上传一个你想处理的 shapefile。图 9-9 展示了我上传的名为“Bodega 海洋实验室和保护区”的文件。当我在我想查看的区域绘制一个多边形（显示为红色）时，我可以下载 GeoTiff 文件。你可以更改坐标来重新定义多边形的形状。
图 9-9. 上传 shapefile 到 EarthExplorer
除了上传一个 shapefile 外，还有几种创建图像的方式：
-
缩放到你希望探索的区域并绘制多边形或圆圈
-
搜索地址（图 9-10）
-
双击地图并选择使用地图按钮
-
选择一个日期范围
你还可以在菜单中输入日期范围或云覆盖范围，如图 9-10 所示。在这个练习中，使用你选择的方法跟随我来到 Bodega 海洋实验室和保护区。将位置输入到搜索条件中。
图 9-10. 搜索感兴趣区域的卫星图像
一旦你到达 Bodega 海洋实验室和保护区，选择数字高程数据，然后选择 SRTM Void Filled，如图 9-11 所示。SRTM 是空中雷达地形测量任务的缩写。Void Filled SRTMs 经过额外处理以填补缺失数据。
图 9-11. 地球测量数字高程 SRTM 数据
如果有符合你条件的图像，它们将作为搜索结果加载（图 9-12）。一旦调整了参数，你将看到数据集的缩略图。找到最适合你需求的数据集并下载 GeoTIFF。将其保存到工作目录中的文件夹或使用绝对路径（完整路径）将其上传到 Spyder 控制台中。^(2)
图 9-12. EarthExplorer 中的搜索结果
哥白尼开放获取中心
我要向你展示的下一个数据资源是哥白尼开放获取中心。你可以通过仪表板中的设置来访问它，类似于你访问 EarthExplorer 的方式。虽然哥白尼的界面可能不如 EarthExplorer 直观，但它提供了一些很棒的数据。尝试搜索 Sentinel 卫星数据（图 9-13）。
图 9-13. 哥白尼开放获取中心 Sentinel 数据
谷歌地球引擎
你在第八章中学习了 GEE，所以这里我只简要提一下如何与 GDAL 结合使用它。GEE 数据以全球范围在云端提供，因此无需下载。你甚至可以使用诸如裁剪卫星数据之类的函数。
运行搜索 GEE DEM 文件。找到 Earth Engine 代码片段 ee.Image(“USGS/sDEP/10m”)，如 图 9-14 所示。您可以复制 JavaScript 代码并粘贴到 GEE 控制台中。
图 9-14. GEE DEM 文件目录
在 GEE 控制台中，选择运行。这将生成图中显示的地图（见 图 9-15）。您可以选择一个多边形并创建几何导入。图层面板允许您更改图层的不透明度或切换图层的显示和隐藏。只需将文件保存为 GeoTIFF，您就有了另一个 DEM 文件选项。
图 9-15. Google Earth Engine DEM .tif 文件
GDAL 可能很复杂，但学习如何使用这个资源库来扩展您的地理空间技能是值得的。
总结
您已经了解了如何在终端和集成开发环境中工作，这两者都根据个人喜好具有各自的用途和通常的可访问性。这些技能突显了这两个选项的效用，以及跨平台应用程序中企业（Esri）和开源社区依赖的强大界面。
在处理地理空间数据和广泛的可用工具时，特别是在开源社区中，重要的是阅读用户文档，以提升您的技能和与数据交互的能力。不要犹豫，如果有问题或见解，请向社区寻求帮助。
^(1) 您可能还记得 绑定 是连接两种编程语言的库，使得一个语言编写的库可以在另一种语言中使用。
^(2) 相对文件路径 是相对于您当前工作目录的路径；绝对文件路径 是从根目录提供的路径。
