地理空間資料分析入門【Shapely+Geopandas】

地理空間資料描述了地球表面上的任何物體或特征。 常見的例子包括：

• 品牌應該在哪裡開設下一家門店？
• 天氣如何影響區域銷售？
• 乘車的最佳路線是什麼？
• 哪個地區受飓風影響最嚴重？
• 冰蓋融化與碳排放有何關系？
• 哪些地區的火災風險最高？

這些問題的答案很有價值，使空間資料技能成為任何資料科學家工具集的重要補充。

讓我們從學習地理空間資料的語言開始。 在本節結束時，我們将了解：

• 矢量與栅格資料
• 地理參考系統 (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"
)           

約翰·斯諾了解到，大多數霍亂死亡病例都集中在 Broad Street 和 Lexington Street 交叉口的一個特定水泵周圍（地圖中間附近的紅色 X）。 他将此次爆發歸因于該水泵的受感染供水。

有趣的是，該地區自 1854 年以來幾乎沒有變化。

原文連結：http://www.bimant.com/blog/python-geospatial-analysis/