Python地理空间数据分析

NSDT工具推荐： Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 - 可编程3D场景编辑器 - REVIT导出3D模型插件 - 3D模型语义搜索引擎 - Three.js虚拟轴心开发包 - 3D模型在线减面 - STL模型在线切割

地理空间数据具有很大的价值。 本文介绍处理地理空间数据所需的思维方式和工具。 它还包括了有史以来第一次空间数据分析的重现：约翰·斯诺 (John Snow) 对 1854 年布罗德街霍乱爆发的调查。‌‌

1、地理空间数据简介

地理空间数据描述地球表面上的任何物体或特征。 常见的例子包括：

• 品牌应该将下一家门店开在哪里？
• 天气如何影响区域销售？
• 开车的最佳路线是什么？
• 哪个地区将受到飓风的打击最严重？
• 冰盖融化与碳排放有何关系？
• 哪些地区发生火灾的风险最高？

这些问题的答案很有价值，使空间数据技能成为任何数据科学家工具集的重要补充。

2、基础

让我们首先学习地理空间数据的语言。 在本节结束时，你将了解：

• 矢量数据与栅格数据
• 地理参考系统 (CRS)
• 地理配准和地理编码之间的区别。

2.1 矢量数据

矢量数据（vector data）代表世界的几何形状。 当你打开导航地图时，你会看到矢量数据。 道路网络、建筑物、餐馆和 ATM 都是具有相关属性的矢量。

注意：矢量是数学对象。 与栅格不同，你可以放大矢量而不会损失分辨率。

矢量数据主要分为三种类型：

• 点
• 线。 连接点创建一条线。
• 多边形。 连接线与封闭区域生成多边形。

我们可以使用矢量来呈现地球表面的特征和属性。 你最常看到存储在 shapefile (.shp) 中的矢量。

定义属性的特定属性通常伴随着向量。 例如，建筑物的属性（例如，其名称、地址、价格、建造日期）可以伴随多边形。

2.2 光栅数据

光栅数据（raster data）是像素网格。 栅格中的每个像素都有一个值，例如颜色、高度、温度、风速或其他测量值。

Google 地图中的默认视图包含矢量，而卫星视图包含拼接在一起的光栅卫星图像。 卫星图像中的每个像素都有一个与之关联的值/颜色。 高程图中的每个像素代表一个特定的高度。 光栅 = 带像素的图像

这些不是你常见的图像。 它们包含我们眼睛可以看到的 RGB 数据，以及来自可见电磁频谱之外的多光谱甚至高光谱信息。 我们可以获得具有多个通道的图像，而不是仅限于 3 个通道/颜色 (RGB)。

肉眼看不见的东西，只吸收电磁频谱的一小部分，却可以在其他电磁频率中显现出来。

光栅 vs. 矢量：

2.3 坐标参考系统 (CRS)

为了确定地球表面的确切位置，我们使用地理坐标系。

例如，尝试在 Google 地图上搜索 37.971441、23.725665。 这两个数字指向一个确切的地方——希腊雅典的帕台农神庙。 这两个数字是 CRS 定义的坐标。

尽管地球是一个 3 维球体，但我们使用经度（南北方向的垂直线）和纬度（东西方向的水平线）的 2 维坐标系来识别地球表面的位置。 将 3D 球体（地球）转换为 2D 坐标系会引入一些扭曲。 我们将在下一节“地图投影”中探讨这些

注意：没有完美的 CRS

CRS 的任何选择都涉及到扭曲以下一项或全部内容的权衡：

• 形状
• 比例/距离
• 区域

很重要！！！ 地理空间分析中的大多数错误都源于为所需操作选择了错误的 CRS。 如果你不想花几天几夜的时间进行调试，请仔细阅读本节！

常见的 CRS 陷阱：

• 混合坐标系：组合数据集时，空间对象必须具有相同的参考系。 请务必将所有内容转换为相同的 CRS。 我们将在下面向您展示如何执行此转换。
• 计算面积：在测量形状的面积之前使用等面积 CRS。
• 计算距离：计算对象之间的距离时使用等距 CRS。

2.4 地图投影

地图投影通过将地球曲面的坐标转换为平面来使地球表面变平。 因为地球不是平坦的（我希望我们在这里同意），所以地球到二维平面的任何投影都只是现实的近似。

事实上，地球是一个大地水准面，意思是一个形状不规则的球体，并不完全是一个球体。 最著名的投影是墨卡托投影。 如上面的 gif 所示，墨卡托投影会使远离赤道的物体膨胀。

这些通货膨胀导致我们的无知令人惊讶地暴露出来，例如美国、中国、印度和欧洲如何融入非洲。

2.5 地理配准

地理配准是将坐标分配给矢量或栅格以将它们投影到地球表面模型上的过程。 它使我们能够创建地图图层。

只需在 Google 地图中点击一下，你就可以从卫星视图无缝切换到道路网络视图。 地理配准使这种转变成为可能。

2.6 地理编码

地理编码是将人类可读的地址转换为一组地理坐标的过程。 例如，希腊雅典的帕台农神庙位于纬度37.988579，经度23.7285437。

有几个库可以为你处理地理编码。 在 Python 中，geopandas 有一个地理编码实用程序，我们将在下面的文章中介绍。

3、Python 地理空间库

在本文中，我们将了解 geopandas 和 shapely，这是使用 Python 进行地理空间分析的两个最有用的库。

• Shapely - 一个允许操作和分析平面几何对象的库。
• Geopandas - 一个库，允许你处理表示表格数据的 shapefile（如 pandas），其中每一行都与一个几何图形相关联。 它提供了对许多空间功能的访问，用于应用几何、绘制地图和地理编码。 Geopandas 内部使用 shapely 来定义几何形状。

3.1 Shapely

你可以将什么放入几何中？

基本形状对象是点、线和多边形，但你也可以在同一对象中定义多个对象。 然后你就有了多点、多线和多多边形。 这些对于由各种几何形状定义的对象非常有用，例如拥有岛屿的国家。

让我们看看它是什么样子的：

import shapely
from shapely.geometry import Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygo

Shapely 通过 x、y 坐标定义一个点，如下所示：

Point(0,0)

我们可以计算形状物体之间的距离，例如两点：

a = Point(0, 0)
b = Point(1, 0)
a.distance(b)

1.0

可以将多个点放置到单个对象中：

MultiPoint([(0,0), (0,1), (1,1), (1,0)])

一系列点形成一个线对象：

line = LineString([(0,0),(1,2), (0,1)])
line

线的长度和边界可通过 length 和bounds 属性获得：

print(f'Length of line {line.length}')
print(f'Bounds of line {line.bounds}')

Length of line 3.6502815398728847
Bounds of line (0.0, 0.0, 1.0, 2.0)

多边形也由一系列点定义：

pol = Polygon([(0,0), (0,1), (1,1), (1,0)])
pol

多边形还具有有用的属性，例如面积：

pol.area

1.0

还有其他一些几何图形交互的有用函数，例如检查多边形 pol 是否与上面的线相交：

pol.intersects(line)

True

我们还可以计算交集：

pol.intersection(line)

但这个对象是什么？

print(pol.intersection(l))

GEOMETRYCOLLECTION (POINT (0 1), LINESTRING (0 0, 0.5 1))

它是一个 GeometryCollection，它是不同类型几何图形的集合。

到目前为止非常简单直观！ 你可以使用 shapely 库做更多事情，因此请务必查看文档。

3.2 Geopandas

处理地理空间数据的另一个工具是 geopandas。 众所周知，pandas DataFrame 代表表格数据集。 同样，geopandas DataFrame 表示具有两个扩展名的表格数据：

• 几何列定义与其余列关联的点、线或多边形。 此列是形状物体的集合。 无论你可以对形状对象做什么处理，也可以对几何对象做同样的处理。
• CRS 是几何列的坐标参考系，它告诉我们点、线或多边形位于地球表面的位置。 Geopandas 将几何图形映射到地球表面（例如 WGS84）。

让我们在实践中看看这一点。

# !pip install geopandas
import matplotlib
import geopandas as gpd

• 加载数据

让我们首先加载 geopandas 附带的数据集，称为“naturalearth_lowres”。 该数据集包括世界上每个国家的几何形状，并附有一些进一步的详细信息，例如人口和 GDP 估计。

world_gdf = gpd.read_file(
    gpd.datasets.get_path('naturalearth_lowres')
)
world_gdf

如果忽略几何列（形状良好的对象），这看起来就像一个常规数据框。 所有列几乎都是不言自明的。

• CRS

数据框还包括一个 CRS，它将几何列中定义的多边形映射到地球表面。

world_gdf.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

在我们的例子中，CRS 是 EPSG:4326。 该 CRS 使用纬度和经度（以度为单位）作为坐标。

CRS 的组成部分

基准面 - 参考系统，在我们的例子中定义了测量的起点（本初子午线）和地球形状的模型（椭球体）。 最常见的基准是 WGS84，但它不是唯一的。

使用区域 - 在我们的例子中，使用区域是整个世界，但是有许多 CRS 针对特定的感兴趣区域进行了优化。

轴和单位 - 通常，经度和纬度以度为单位测量。 x、y 坐标的单位通常以米为单位。

让我们看一下必须更改 CRS 的应用程序。

我们来测算一下各个国家的人口密度吧！ 我们可以测量每个几何体的面积，但请记住，我们首先需要转换为以米为单位的等面积投影。

world_gdf = world_gdf.to_crs("+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
world_gdf.crs

<Projected CRS: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +un ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Eckert IV
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

现在，我们可以通过将人口估计值除以面积来计算每个国家的人口密度。

注意：我们可以像访问常规列一样访问几何图形的区域。 尽管没有列包含几何区域，但区域是几何对象的属性。

world_gdf['pop_density'] = world_gdf.pop_est / world_gdf.area * 10**6

world_gdf.sort_values(by='pop_density', ascending=False)

请注意，几何对象现在的值的单位与以前完全不同。

只需查看上面的数据框，我们就可以快速识别异常值。 孟加拉国人口密度约为
。 南极洲的人口密度接近于零，广阔的空间里只有810人居住。

不过，将地图可视化总是更好。 那么让我们想象一下吧！

• 可视化

我们可以像 pandas 数据框一样在 world_gdf 上调用 .plot() ：

figsize = (20, 11)

world_gdf.plot('pop_density', legend=True, figsize=figsize);

上面的地图看起来不太有帮助，所以让我们通过执行以下操作来使其更好：

• 更改为墨卡托投影，因为它更熟悉。 使用参数 to_crs('epsg:4326') 可以做到这一点。
• 将颜色条转换为对数刻度，可以使用  matplotlib.colors.LogNorm(vmin=world.pop_density.min(), vmax=world.pop_density.max()) 来实现

我们可以将不同的参数传递给绘图函数，就像直接在 matplotlib 上一样。

norm = matplotlib.colors.LogNorm(vmin=world_gdf.pop_density.min(), vmax=world_gdf.pop_density.max())

world_gdf.to_crs('epsg:4326').plot("pop_density", 
                                   figsize=figsize, 
                                   legend=True,  
                                   norm=norm);

到目前为止，我们已经了解了 shapely 和 geopandas 的基础知识，但现在是时候进行完整的案例研究了。

4、案例研究：1854 年霍乱爆发的根源

我们将使用现代 Python 工具重做 John Snow 的分析，确定 1854 年伦敦布罗德街霍乱爆发的来源。 与《权力的游戏》中的对手相反，伦敦的约翰·斯诺现在做了一件事：霍乱的根源。 他在第一次地理空间分析中学会了这一点！

对于此示例，我们将使用 Robin 博客中的数据。 罗宾负责将斯诺的原始地图和数据数字化。

让我们首先检索数据并将其解压缩到当前目录中：

!wget http://www.rtwilson.com/downloads/SnowGIS_v2.zip

!unzip SnowGIS_v2.zip

让我们看看里面有什么：

!ls SnowGIS/

Cholera_Deaths.dbf
Cholera_Deaths.prj
Cholera_Deaths.sbn
Cholera_Deaths.sbx
Cholera_Deaths.shp
Cholera_Deaths.shx
OSMap.tfw
OSMap.tif
OSMap_Grayscale.tfw
OSMap_Grayscale.tif
OSMap_Grayscale.tif.aux.xml
OSMap_Grayscale.tif.ovr
Pumps.dbf
Pumps.prj
Pumps.sbx
Pumps.shp
Pumps.shx
README.txt
SnowMap.tfw
SnowMap.tif
SnowMap.tif.aux.xml
SnowMap.tif.ovr

暂时忽略文件扩展名，让我们看看这里有什么。

矢量数据：

• Cholera_Deaths ：给定空间坐标处的死亡人数
• Pumps：水泵的位置

我们可以忽略矢量数据的其他文件，只处理“.shp”文件。 .shp 称为 shapefile，是矢量对象的标准格式。

栅格数据：

• OSMap_Grayscale ：栅格 - 来自 OpenStreet Maps (OSM) 的区域的地理参考灰度图
• OSMap ：栅格 - 来自 OpenStreet Maps (OSM) 的区域地理参考地图
• SnowMap ：栅格 - 数字化和地理参考约翰·斯诺的原始地图

我们可以忽略栅格数据的其他文件，只处理“.tif”文件。 “.tif”是存储光栅和图像数据的最常见格式。

!pip install rasterio
!pip install contextily

我们已经导入了 geopandas 和 matplotlib，因此剩下的分析需要导入上下文来绘制地图。 现在让我们导入它们：

import contextily as ctx

4.1 读入数据

让我们将 Cholera_Death.shp 和 Pumps.shp 文件读入 geopandas：

deaths_df = gpd.read_file('SnowGIS/Cholera_Deaths.shp')
pumps_df = gpd.read_file('SnowGIS/Pumps.shp')

deaths_df.head()

输出看起来与 pandas 数据框完全相同。 与 geopandas 数据框的唯一区别是几何列，这是我们矢量数据集的本质。 在我们的例子中，它包括约翰·斯诺记录的死亡点坐标。

让我们看看 CRS 数据是什么样的：

deaths_df.crs

<Projected CRS: EPSG:27700>
Name: OSGB 1936 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.
- bounds: (-9.0, 49.75, 2.01, 61.01)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: OSGB 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich

另一个区别是正确定义的 shapefile 包含阐明其坐标参考系统 (CRS) 的元数据。 在本例中，它是 EPSG：27700。

现在让我们简单地看一下泵数据：

pumps_df

同样，pumps_df 持有 Broad Street 附近水泵的位置。

以下是泵的 CRS 数据：

pumps_df.crs

<Projected CRS: EPSG:27700>
Name: OSGB 1936 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.
- bounds: (-9.0, 49.75, 2.01, 61.01)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: OSGB 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich

注意：处理地理空间数据时，你应该确保所有源都具有相同的 CRS。 我怎么强调都不为过。 这可能是处理地理空间数据时最常见的错误来源。

4.2 绘制疫情爆发图

现在，我们可以在伦敦布罗德街的地图上绘制死亡人数和水泵数据。

我们将首先绘制死亡图表来开始构建情节：

ax = deaths_df.plot(column='Count', alpha=0.5, edgecolor='k', legend=True)

参考 ax，我们可以绘制泵的位置，并用红色 X 标记它们。我们还可以将图形放大：

ax = deaths_df.plot(column='Count', figsize=(15, 15), alpha=0.5, edgecolor='k', legend=True)
pumps_df.plot(ax=ax, marker='x', color='red', markersize=50)

<AxesSubplot:>

我们现在想要在数据下方显示伦敦布罗德街的地图。 这是我们可以使用上下文来读取 CRS 数据的地方：

ax = deaths_df.plot(column='Count', figsize=(15, 15), alpha=0.5, edgecolor='k', legend=True)
pumps_df.plot(ax=ax, marker='x', color='red', markersize=50)

ctx.add_basemap(
    ax,
    # CRS definition. Without the line below, the map stops making sense
    crs=deaths_df.crs.to_string(),
)

现在让我们看一下相同的数据，但是是在约翰·斯诺的原始地图上。 我们可以通过将源参数更改为 SnowMap.tif 来做到这一点，如下所示：

ax = deaths_df.plot(column='Count', figsize=(15, 15), alpha=0.5, edgecolor='k', legend=True)
pumps_df.plot(ax=ax, marker='x', color='red', markersize=50);

ctx.add_basemap(ax,
    crs=deaths_df.crs.to_string(),
    # Using the original map, hand-drawn by Snow
    source="SnowGIS/SnowMap.tif"
)

约翰·斯诺了解到，大多数霍乱死亡病例都集中在布罗德街和列克星敦街交叉口的某个特定水泵周围（地图中间附近的红色 X）。 他将此次疫情的爆发归因于该泵的供水系统受到感染。

有趣的是，自 1854 年以来，该地区几乎没有发生任何变化。

原文链接：Analyze Geospatial Data in Python: GeoPandas and Shapely

BimAnt翻译整理，转载请标明出处