Python GIS 地理信息数据分析入门：GeoPandas 和 Shapely

翻译自： Analyze Geospatial Data in Python: GeoPandas and Shapely
原作者： Ioannis Prapas

更好的排版以及代码输出，可点击下方和鲸社区链接一键跑通

地理空间数据介绍

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

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

这些问题的回答是有价值的，这使得空间数据分析技能成为数据科学家技能的重要补充。

基础知识

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

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

矢量

矢量数据表示世界上的几何图形。当你打开一张导航地图时，你会看到矢量数据。道路网络、建筑物、餐馆和自动取款机都是矢量，并有其相关的属性。

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

有三种主要的矢量数据类型：

• 点
• 线：将点连接起来就形成了一条线
• 多边形：将线连起来形成的封闭的区域

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

定义物体的具体属性一般会伴随着向量。例如，建筑物的属性（如名称、地址、价格、建造日期）可以伴随着一个多边形。

光栅

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

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

这些不是你常见的图像。它们包含我们的眼睛可以看到的RGB数据，以及来自可见电磁波谱以外的多光谱甚至高光谱信息。我们不再局限于只有3个通道/颜色（RGB），而是可以得到许多通道的图像。

肉眼看不见的东西只吸收了电磁波谱的一小部分，他们可以在其他电磁频率下显示出来。

光栅VS矢量表

坐标参考系统（CRS）

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

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

尽管地球是一个三维球体，但我们使用经度（南北走向的垂直线）和纬度（东西走向的水平线）的二维坐标系统来确定地球表面的位置。将一个三维球体（地球仪）转换成一个二维坐标系会带来一些扭曲。我们将在下一节 "地图投影 "中探讨这些失真问题。

注意：没有一个CRS是完美的

对CRS的任何选择都涉及到对以下一个或所有方面的扭曲的权衡：

• 形状
• 比例/距离
• 面积

非常重要 ！! 地理空间分析中的大多数错误来自于为所需操作选择了错误的CRS。如果你不想花费几天几夜的时间进行调试，请彻底阅读本节内容

常见的CRS误区 ：

• 混合坐标系：当结合数据集时， 空间对象必须具有相同的参考系统 。请确保将所有的东西转换为相同的CRS。我们在下面向你展示如何进行这种转换。
• 计算面积：在测量一个形状的面积之前，使用一个等面积的CRS。
• 计算距离：在计算物体之间的距离时，使用等距CRS。

地图投影

地图投影通过将坐标从地球的弯曲表面转化为一个平面，使地球的表面变得平坦。因为地球不是平的（我希望我们在这里达成共识），任何将地球投影到一个二维平面的做法都只是对现实的一种近似。

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

这些膨胀导致我们很容易忽略一些事实，比如美国、中国、印度和欧洲都比非洲小。

地理参考

地理参考是为矢量或光栅分配坐标的过程，以便将它们投射到地球表面的模型上。它使我们能够创建地图的图层。

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

地理编码

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

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

Python地理空间库

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

Shapely - 一个允许操作和分析平面几何对象的库。

pip install shapely

Geopandas - 一个允许你处理代表表格数据的shapefiles的库（如pandas），其中每一行都与一个几何体相关。它提供了对许多空间函数的访问，用于应用几何图形、绘制地图和地理编码。Geopandas内部使用shapely来定义几何图形。

Shapely

你可以把什么放进几何体（geometry ）？

基本的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)

多个点可以被放置在一个物体中：

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

一系列的点构成了一个线对象：

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

一条线的长度和边界可通过长度和边界属性获得：

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

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

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

多边形也有一些属性，如面积：

pol.area

还有其他一些有用的功能，在这些功能中，几何体是相互作用的，例如检查多边形plo是否与上面的line相交：

pol.intersects(line)

我们也可以计算出交叉点：

pol.intersection(line)

那么它的数据类型是什么呢？

print(pol.intersection(line))

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

到目前为止，它是非常直接和直观的！你可以用shapely库做更多的事情，所以一定要查看文档。

Geopandas基础知识

另一个用于处理地理空间数据的工具是geopandas。正如我们所知，pandas 的 DataFrames 表示表格数据集。同样地，geopandas 的 DataFrames 表示的是有两个扩展功能的表格数据集：

• geometry列定义了一个与其他列相关的点、线或多边形。这一列是一个shapely对象的集合。你能对shapley的形状对象做什么，你也能对geometry
  对象做什么。
• CRS列是geometry
  列的坐标参考系统，它告诉我们一个点、线或多边形在地球表面的位置。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

177 rows × 6 columns

如果你忽略geometry列（一个shapely 对象），这看起来就像一个普通的DataFrames，所有列的意义都很容易明白。

CRS

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

world_gdf.crs

在我们的案例中，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

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

注意 ：我们可以像访问普通列一样访问几何体的面积。虽然没有列包含几何体的面积，但面积是几何体对象的一个属性。

world_gdf['pop_density'] = world_gdf.pop_est / world_gdf.area * 10**6
world_gdf.sort_values(by='pop_density', ascending=False)

177 rows × 7 columns

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

仅仅看一下上面的数据框架，我们就可以迅速识别出异常值。孟加拉国的人口密度约为$1175 persons/km^2$. 南极洲的人口密度几乎为零，只有810人生活在一个广阔的空间里。

不过，将地图可视化总是更好。所以让我们可视化吧!

可视化

我们可以在world_gdf上调用.plot()，就像调用pandas的dataframe一样。

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的基本知识，但现在是时候进入一个完整的案例研究。

案例研究：1854年霍乱爆发的来源

我们将使用现代Python工具重做约翰-斯诺的分析，确定1854年伦敦布罗德街霍乱爆发的来源。

在这个例子中，我们将使用罗宾的博客中的数据。罗宾做了将斯诺的原始地图和数据数字化的工作。

让我们首先检索数据并在我们的当前目录中解压：

# !wget http://www.rtwilson.com/downloads/SnowGIS_v2.zip
# !unzip SnowGIS_v2.zip
!ls SnowGIS/

暂时不考虑文件的扩展名，让我们看看我们这里有什么。

矢量数据

Cholera_Deaths: 某一空间坐标上的死亡人数
Pumps：水泵的位置 我们可以忽略矢量数据的其他文件，只处理'.shp'文件。.shp被称为shapefiles，是矢量对象的标准格式。

光栅数据

OSMap_Grayscale: 光栅--来自OpenStreet Maps (OSM)的该地区的地理参照灰度地图
OSMap: 光栅--来自OpenStreet Maps (OSM)的该地区的地理参考图
SnowMap: 光栅 -- 数字化和地理参照的约翰-斯诺的原始地图

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

我们已经导入了geopandas和matplotlib，所以在剩下的分析中我们只需要导入用于绘制地图的contextily。现在让我们来导入这些：

import contextily as ctx

读入数据

让我们把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

另一个区别是，正确定义的形状文件包括阐明其坐标参考系统（CRS）的元数据。在这种情况下，它是EPSG:27700。

现在让我们简要地看一下水泵的数据：

pumps_df

同样，pump_df保存了Broad Street附近的水泵的位置。

以下是水泵的CRS数据：

pumps_df.crs

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

绘制疫情图

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

我们将首先通过绘制死亡人数来建立该图：

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)

我们现在想在数据下面显示一张伦敦布罗德街的地图。这时我们可以使用contextily来读取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(
    # CRS definition. Without the line below, the map stops making sense
    crs=deaths_df.crs.to_string(),
)

现在我们在约翰-斯诺的原始地图上看看同样的数据。我们可以通过改变source参数为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);