Python 地理空间数据入门

在本教程中，您将了解两个用于处理地理空间数据的常用软件包：geopandas 和 Shapely。随后，您将把这两个软件包应用于使用 Python 读取地理空间数据，并绘制飓风佛罗伦萨在 8 月 30 日至 9 月 18 日期间的轨迹。

什么是地理空间数据？

Spatial data、Geospatial data、GIS data 或 geodata，都是指根据地理坐标系来标识建筑、街道、城镇、城市、国家等物理对象地理位置的数值数据。通过空间数据，您不仅可以了解位置，还能获知任何对象的长度、大小、面积或形状。您可以获取的一类空间数据示例包括：对象的坐标，如纬度、经度和海拔。可以使用地理信息系统（GIS）或其他专业软件来访问、可视化、操作和分析地理空间数据。

为什么要用地理空间数据？

地理空间数据在我们的日常生活中有着广泛应用。一个例子是您用于从一个位置导航到另一个位置的地图应用。另一个您每天都可能看到的例子是在天气频道上。

没错，地理空间数据用于表示某物相对于其周围事物的位置关系：城市地图上的您的住址、世界地图上的飓风位置，等等。因此，在本教程中，您将关注这场具有破坏性的飓风佛罗伦萨并跟踪其位置。

软件包要求

不要跳过这里。这是重要部分。如果您不确定是否满足要求，请再检查一遍。首先，您需要安装以下所有软件包。

• Pandas：提供数据结构和数据分析工具
• Numpy：Python 科学计算的基础包
• SciPy：（发音为“Sigh Pie”）面向数学、科学与工程的 Python 开源软件生态系统
• RTree：libspatialindex 的 ctypes Python 封装，提供多种高级空间索引功能
• GDAL：栅格与矢量地理空间数据格式的转换库
• Fiona：用于读取和写入空间数据文件
• Shapely：几何对象、谓词与操作
• GeoPandas：扩展 pandas 的数据类型以支持对几何类型的空间操作
• PySAL：用 Python 编写的空间分析函数库，旨在支持高级应用开发
• Matplotlib：Python 2D 绘图库
• Missingno：Python 缺失数据可视化模块

按上述顺序安装所有软件包以确保一切正常运行。有些包是其他包的前置依赖，例如：安装 GeoPandas 需要 Shapely，而安装 Shapely 之前应先安装 RTree、GDAL 和 Fiona。最简单的安装方式为：pip install PACKAGE_NAME

如果您使用 Windows 且无法通过上述方式安装，请查看此网站下载相关包，然后执行 pip install PATH_TO_PACKAGE。

关于数据

您可能知道，可怕的飓风佛罗伦萨刚刚横扫了美国东海岸的部分地区，造成估计170 亿美元的损失。本教程将帮助您了解它从何而来、何时何地增强，并借助 Python 更深入地理解这一自然灾害并对其进行分析。有许多网站提供关于这场飓风的信息。例如，这个网站提供了 1902 年至 2018 年间美国多场风暴和飓风的数据，因此有大量可供您后续使用的数据。在本教程中，您只会使用飓风佛罗伦萨数据

您还将使用来自互联网的美国地图地理空间数据。Eric Celeste 的博文包含美国县与州各种边界数据文件。本教程使用的是 US States、5m、GeoJSON 文件。

 # Load all importance packages
import geopandas
import numpy as np
import pandas as pd
from shapely.geometry import Point

import missingno as msn

import seaborn as sns
import matplotlib.pyplot as plt

% matplotlib inline技术支持

首先，让我们看看第一个地理空间数据框：美国各州地理数据

# Getting to know GEOJSON file:
country = geopandas.read_file("data/gz_2010_us_040_00_5m.json")
country.head()
技术支持

检查刚加载的数据框类型，您会看到它是 Geo Data Frame，具备 Pandas DataFrame 的所有常规特性。

type(country)
技术支持

geopandas.geodataframe.GeoDataFrame
技术支持

检查包含坐标的列的数据类型：它是 GeoSeries。

type(country.geometry)
技术支持

geopandas.geoseries.GeoSeries
技术支持

GeoSeries 中的每个值都是一个 Shapely 对象。它可以是：

• Point
• Line
• Polygon
• MultiPolygon

每种对象可用于不同类型的物理对象，例如：Point 表示建筑，Line 表示街道，Polygon 表示城市，MultiPolygon 表示包含多个城市的国家。关于每种几何对象的更多信息，请阅读此文：https://shapely.readthedocs.io/en/stable/manual.html#geometric-objects

type(country.geometry[0])
技术支持

shapely.geometry.multipolygon.MultiPolygon
技术支持

与 Pandas DataFrame 类似，GeoDataFrame 也有 plot 属性，它利用数据框中的几何信息来绘制地图：

country.plot()
技术支持

<matplotlib.axes._subplots.AxesSubplot at 0x1cfe68c1358>
技术支持

如您所见，美国地图相对画布较小。原因是信息包含阿拉斯加、夏威夷和波多黎各，它们分布较散。就本教程而言，您可以排除阿拉斯加和夏威夷，因为飓风并未接近这两个州。您还可以添加图形尺寸和颜色以自定义绘图：

# Exclude Alaska and Hawaii for now
country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(figsize=(30,20), color='#3B3C6E');
技术支持

并不难，对吧！现在您已有美国地图，让我们加载飓风数据：

florence = pd.read_csv('data/florence.csv')
florence.head()
技术支持

探索性数据分析

这始终是加载任何数据集时要做的第一件事：

• 检查信息与数据类型
• 是否存在缺失值
• 统计信息

florence.info()
技术支持

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 105 entries, 0 to 104
Data columns (total 11 columns):
AdvisoryNumber    105 non-null object
Date              105 non-null object
Lat               105 non-null float64
Long              105 non-null float64
Wind              105 non-null int64
Pres              105 non-null int64
Movement          105 non-null object
Type              105 non-null object
Name              105 non-null object
Received          105 non-null object
Forecaster        104 non-null object
dtypes: float64(2), int64(2), object(7)
memory usage: 9.1+ KB
技术支持

使用 missingno 包检查缺失值。这个包通过可视化展示缺失数据非常有用。如下所示，只有 “Forecaster” 列缺失一个值，而本教程并不需要它，所以可以暂时忽略。

# Notice you can always adjust the color of the visualization
msn.bar(florence, color='darkolivegreen');
技术支持

查看一些统计信息，其中一些会很有用，比如该飓风的平均风速、最大和最小风速等。

# Statistical information
florence.describe()
技术支持

对于大多数数据，您需要清理并只保留需要处理的部分。这里，您只需要时间、坐标（纬度和经度）、风速、气压和名称。Movement 与 Type 可选，其他可以删除。

# dropping all unused features:
florence = florence.drop(['AdvisoryNumber', 'Forecaster', 'Received'], axis=1)
florence.head()
技术支持

通常，如果单独绘制数据，不必对坐标做额外处理。不过，如果您希望其看起来与地图上的方向一致，检查经纬度就很重要。此处经度为西经，您需要在数字前加 “-” 以正确绘制数据：

# Add "-" in front of the number to correctly plot the data:
florence['Long'] = 0 - florence['Long']
florence.head()
技术支持

接着，您可以将纬度与经度合并以创建飓风坐标，随后将其转为 GeoPoint 以用于可视化。

# Combining Lattitude and Longitude to create hurricane coordinates:
florence['coordinates'] = florence[['Long', 'Lat']].values.tolist()
florence.head()
技术支持

# Change the coordinates to a geoPoint
florence['coordinates'] = florence['coordinates'].apply(Point)
florence.head()
技术支持

检查 florence 数据框及其 coordinates 列的类型。它们分别是 pandas 的 DataFrame 与 Series。

type(florence)
技术支持

pandas.core.frame.DataFrame
技术支持

type(florence['coordinates'])
技术支持

pandas.core.series.Series
技术支持

将数据转换为地理空间数据后，我们再检查一次 florence 数据框和 coordinates 列的类型。现在分别是 GeoDataFrame 和 GeoSeries。

# Convert the count df to geodf
florence = geopandas.GeoDataFrame(florence, geometry='coordinates')
florence.head()
技术支持

type(florence)
技术支持

geopandas.geodataframe.GeoDataFrame
技术支持

type(florence['coordinates'])
技术支持

geopandas.geoseries.GeoSeries
技术支持

请注意，尽管现在是 GeoDataFrame 和 GeoSeries，它们仍然像普通的 DataFrame 和 Series 那样工作。这意味着您仍可对数据框执行过滤、分组，或对列提取最小值、最大值或平均值。

# Filtering from before the hurricane was named.
florence[florence['Name']=='Six']
技术支持

# Groupping by name to see how many names it has in the data set:
florence.groupby('Name').Type.count()
技术支持

Name
FLORENCE     6
Florence    85
SIX          4
Six         10
Name: Type, dtype: int64
技术支持

求飓风佛罗伦萨的平均风速：

print("Mean wind speed of Hurricane Florence is {} mph and it can go up to {} mph maximum".format(round(florence.Wind.mean(),4),
                                                                                         florence.Wind.max()))
技术支持

Mean wind speed of Hurricane Florence is 74.4286 mph and it can go up to 140 mph maximum
技术支持

因此，飓风佛罗伦萨的平均风速为每小时 74.43 英里（每小时 119.78 公里），最大风速为每小时 140 英里（每小时 225.308 公里）。为了感受这种风速有多可怕，由英国皇家海军制定的蒲福风级展示了风对水面和陆地的影响。当风速达到每小时 48 至 55 英里时，就已经可以折断并连根拔起树木，并造成“相当严重的结构性破坏”。

在这种时刻，您绝不想身处其中。

可视化

与 pandas DataFrame 类似，GeoDataFrame 也有 .plot 属性。不过，该属性会利用 GeoDataFrame 中的坐标来绘图。我们来看看：

florence.plot(figsize=(20,10));
技术支持

发生了什么？您只看到一堆没有参照的点。这有问题吗？
没有，一切正常。因为该数据框只包含飓风佛罗伦萨在各个时间点的位置（坐标）信息，所以我们只能在空白底图上绘制这些位置。

因此，下一步是在美国地图上绘制飓风位置，以查看它在哪里登陆以及当时有多强。为此，您将使用（在开头加载的）美国地图坐标作为底图，并在其上叠加飓风佛罗伦萨的位置。

 # Plotting to see the hurricane overlay the US map:
fig, ax = plt.subplots(1, figsize=(30,20))
base = country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(ax=ax, color='#3B3C6E')

# plotting the hurricane position on top with red color to stand out:
florence.plot(ax=base, color='darkred', marker="*", markersize=10);技术支持

看起来很棒！现在，我们通过更多细节来完善它，例如：

• 添加标题
• 根据风速为飓风位置着色，以查看其袭击各城市时的强度
• 移除坐标轴
• 添加图例
• 将结果保存为图像文件以便后续使用

 fig, ax = plt.subplots(1, figsize=(20,20))
base = country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(ax=ax, color='#3B3C6E')
florence.plot(ax=base, column='Wind', marker="<", markersize=10, cmap='cool', label="Wind speed(mph)")
_ = ax.axis('off')
plt.legend()
ax.set_title("Hurricane Florence in US Map", fontsize=25)
plt.savefig('Hurricane_footage.png',bbox_inches='tight');技术支持

可见，飓风在东海岸近海时最强。接近陆地后开始减弱，但当风速在每小时 60 至 77 英里的范围内时，仍会造成可怕的破坏。

结语

做得好！您已经学习了在 Python 中处理地理空间数据的一些必要步骤。您还学会了如何绘制地理空间数据，并通过自定义形状、颜色和叠加来讲述一个故事。如果您想练习技能，网上有大量地理空间数据可供尝试。一个示例就是我上面提供的天气网站。

如需参考资料，请查看 Karim Bahgat（2015）所著的“Python Geospatial Development Essentials”一书以获得更深入的指导。

如果您想与我联系，可以发邮件至 [email protected]，或通过LinkedIn 与我建立联系。祝您学习愉快，注意安全！

如果您希望进一步学习 Python 中的地理空间数据，请参加 DataCamp 的Visualizing Geospatial Data in Python课程。

以下是本教程使用的数据：

• export.csv
• florence.csv
• gz_2010_us_040_00_5m
• hurricane_data