Python地理空间分析快速入门

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

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

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

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

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

在这里插入图片描述

推荐：用 NSDT设计器快速搭建可编程3D场景。

1、矢量数据

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

注意：矢量是数学对象。与光栅不同，你可以在不损失分辨率的情况下放大矢量。
矢量数据主要分为三种类型：

点
线。连接点创建一条线。
多边形。连接具有封闭区域的线会生成一个多边形。

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

矢量通常附带有一些属性。例如，建筑物的属性（例如，其名称、地址、价格、建造日期）可以伴随多边形。

在这里插入图片描述

2、栅格数据

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

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

在这里插入图片描述

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

肉眼看不见的东西，只吸收一小部分电磁波谱，在其他电磁波频率下就可以显露出来。
在这里插入图片描述

3、坐标参考系统 (CRS)

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

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

尽管地球是一个 3 维球体，但我们使用经度（南北垂直线）和纬度（东西水平线）的二维坐标系来确定地球表面上的位置。将 3D 球体（地球）转换为 2D 坐标系会引入一些变形。我们将在下一节地图投影中探讨这些扭曲。
在这里插入图片描述

注意：没有 CRS 是完美的

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

形状
尺度/距离
区域

重要提示！！！地理空间分析中的大多数错误来自于为所需的操作选择了错误的 CRS。如果不想花几天几夜调试，请通读本节！

常见的 CRS 陷阱包括：

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

4、地图投影

地图投影通过将坐标从地球曲面转换为平面来使地球表面变平。因为地球不是平的，地球在二维平面上的任何投影都只是对现实的近似。
在这里插入图片描述

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

5、地理配准

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

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

6、地理编码

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

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

7、Python 地理空间库

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

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

8、Shapely入门

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

让我们看看它看起来如何：

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

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 库做更多的事情，一定要查看文档。

9、Geopandas 入门

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

几何列定义与其余列关联的点、线或多边形。此列是一组Shapely的对象。
CRS 是几何列的坐标参考系统，它告诉我们点、线或多边形在地球表面的位置。 Geopandas 将几何图形映射到地球表面（例如 WGS84）。

让我们在实践中看看这个。

9.1 加载数据

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

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

在这里插入图片描述

如果您忽略几何列（一个Shapely对象），这看起来就像一个普通的数据框。所有的列都是不言自明的。

9.2 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 的组成部分包括：

Datum - 参考系统，在我们的例子中定义了测量起点（本初子午线）和地球形状模型（椭球体）。最常见的 Datum 是 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 人生活在广阔的空间中。

不过，可视化地图总是更好。所以让我们想象一下！

9.3 可视化

我们可以像 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 上一样将不同的参数传递给 plot 函数。

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);

在这里插入图片描述

10、案例研究：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 ：数字化和地理参考 John Snow 的原始地图

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

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

import contextily as ctx

10.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 数据框的唯一区别是几何列，这是我们矢量数据集的本质。在我们的例子中，它包括 John Snow 记录的死亡点坐标。

让我们看看 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。我怎么强调都不过分。在处理地理空间数据时，它可能是所有错误中最常见的来源。

10.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(),
)

在这里插入图片描述

现在让我们看一下相同的数据，但是是在 John Snow 的原始地图上。我们可以通过将源参数更改为 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"
)

在这里插入图片描述