Python地理空间分析

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

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

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

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

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

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

约翰·斯诺了解到,大多数霍乱死亡病例都集中在 Broad Street 和 Lexington Street 交叉口的一个特定水泵周围(地图中间附近的红色 X)。 他将此次爆发归因于该水泵的受感染供水。

有趣的是,该地区自 1854 年以来几乎没有变化。


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

bBimAnt翻译整理,转载请标明出处