本文主要介紹GeoPandas的使用要點。GeoPandas是一個Python開源專案，旨在提供豐富而簡單的地理空間資料處理介面。GeoPandas擴充套件了Pandas的資料型別，並使用matplotlib進行繪圖。GeoPandas官方倉庫地址為：GeoPandas。GeoPandas的官方檔案地址為：GeoPandas-doc。本文主要參考GeoPandas Examples Gallery。GeoPandas的基礎使用見Python繪製資料地圖1-GeoPandas入門指北。GeoPandas的視覺化入門見Python繪製資料地圖2-GeoPandas地圖視覺化。

本文所有程式碼見：Python-Study-Notes

GeoPandas推薦使用Python3.7版本及以上，執行環境最好是linux系統。GeoPandas安裝命令如下：

pip install geopandas

如果上述命令安裝出問題，則推薦使用conda安裝GeoPandas，命令如下：

conda install geopandas

或：

conda install --channel conda-forge geopandas

除了GeoPandas需要安裝，以下第三方庫也需要安裝：

pip install mapclassify
pip install matplotlib_scalebar
pip install rtree
pip install contextily

# jupyter notebook環境去除warning
import warnings
warnings.filterwarnings("ignore")

# 檢視geopandas版本
import geopandas as gpd

gpd.__version__

'0.13.2'

1 分級統計圖Choropleth

分級統計圖Choropleth是一種表示地理區域內資料分佈的視覺化圖表。它將地圖劃分為不同的區域，並使用顏色或陰影的不同程度來顯示該區域的資料值。通常，分級統計圖用於顯示人口統計、自然資源分佈等資料。分級統計圖以幫助觀察者更容易地理解資料在地理空間上的分佈情況和變化趨勢，有助於制定決策和規劃相關工作。

import geopandas as gpd
from geopandas import read_file

# pip install mapclassify
import mapclassify
mapclassify.__version__

'2.5.0'

# 讀取四川地圖資料，資料來自DataV.GeoAtlas，將其投影到EPSG:4573
data = gpd.read_file('https://geo.datav.aliyun.com/areas_v3/bound/510000_full.json').to_crs('EPSG:4573')
data.head()

簡單分級統計

以下程式碼通過scheme分級統計四川省各地級市所包含區縣數。

ax = data.plot(
    column="childrenNum",
    scheme="QUANTILES", # 設定分層設色標準
    edgecolor='lightgrey', 
    k=7, # 分級數量
    cmap="Blues",
    legend=True,
    # 通過fmt設定位數
    legend_kwds={"loc": "center left", "bbox_to_anchor": (1, 0.5),"fmt": "{:.2f}"}
)

# 顯示各地級市包含區縣數量
for index in data.index:
    x = data.iloc[index].geometry.centroid.x
    y = data.iloc[index].geometry.centroid.y
    name = data.iloc[index]["childrenNum"]
    ax.text(x, y, name, ha="center", va="center",color='red')

# 檢視label，也就是分級區間
labels = [t.get_text() for t in ax.get_legend().get_texts()]
labels

[' 3.00,  5.00',
 ' 5.00,  6.00',
 ' 6.00,  6.57',
 ' 6.57,  7.43',
 ' 7.43,  9.29',
 ' 9.29, 13.57',
 '13.57, 20.00']

# 檢視各個分級標準和區間數量，一般是左開右閉
res = mapclassify.Quantiles(data["childrenNum"], k=7)
res

Quantiles

   Interval      Count
----------------------
[ 3.00,  5.00] |     5
( 5.00,  6.00] |     4
( 6.00,  6.57] |     0
( 6.57,  7.43] |     3
( 7.43,  9.29] |     3
( 9.29, 13.57] |     3
(13.57, 20.00] |     3

間隔展示

ax = data.plot(
    column="childrenNum",
    scheme="BoxPlot", 
    edgecolor='k',
    cmap="OrRd", # 設定分層設色標準
    legend=True,
    # 通過interval設定是否展示區間間隔
    legend_kwds={"loc": "center left", "bbox_to_anchor": (1, 0.5), "interval": True}
)

# 顯示各地級市包含區縣數量
for index in data.index:
    x = data.iloc[index].geometry.centroid.x
    y = data.iloc[index].geometry.centroid.y
    name = data.iloc[index]["childrenNum"]
    ax.text(x, y, name, ha="center", va="center",color='red')

分類展示

以數值分類的方式展示資料，其中區縣數量為20的地級市為成都市。

ax = data.plot(
    column="childrenNum",
    categorical=True, # 以數值分類的方式展示
    legend=True,
    cmap="tab20",
    # 對於分類資料，fmt設定無用
    legend_kwds={"loc": "center left", "bbox_to_anchor": (1, 0.5), "fmt": "{:.0f}"},
)  

# 顯示各地級市包含區縣數量
for index in data.index:
    x = data.iloc[index].geometry.centroid.x
    y = data.iloc[index].geometry.centroid.y
    name = data.iloc[index]["childrenNum"]
    ax.text(x, y, name, ha="center", va="center",color='red')

自定義分級

# 自定義分級標準
def custom(value):
    # 設定ABC三個等級
    level = None
    if value > 15:
       level = 'A'
    elif value > 8:
        level = 'B'
    else:
        level = 'C'
    return level

# 根據自定義函數對映為新的列
data['level'] = data['childrenNum'].apply(custom)
data.head()

ax = data.plot(
    column="level",
    categorical=True, # 以數值分類的方式展示
    legend=True,
    cmap="coolwarm",
    # 對於分類資料，fmt設定無用
    legend_kwds={"loc": "center left", "bbox_to_anchor": (1, 0.5), "fmt": "{:.0f}"},
)  

# 顯示各地級市包含區縣數量
for index in data.index:
    x = data.iloc[index].geometry.centroid.x
    y = data.iloc[index].geometry.centroid.y
    name = data.iloc[index]["childrenNum"]
    ax.text(x, y, name, ha="center", va="center",color='red')

2 通過DataFrame建立GeoDataFrame

基於經緯度資料

GeoDataFrame有一個geometry列，我們可以通過經緯度資料Latitude和Longitude建立該列。

import pandas as pd
# 生成關於南美城市的dataframe資料
df = pd.DataFrame(
    {
        "City": ["Buenos Aires", "Brasilia", "Santiago", "Bogota", "Caracas"],
        "Country": ["Argentina", "Brazil", "Chile", "Colombia", "Venezuela"],
        "Latitude": [-34.58, -15.78, -33.45, 4.60, 10.48],
        "Longitude": [-58.66, -47.91, -70.66, -74.08, -66.86],
    }
)

df

# 將dataframe轉換為GeoDataFrame
import geopandas as gpd
gdf = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude), crs="EPSG:4326"
)
gdf

# 在南美地圖上展示
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
# 定位到南美
ax = world.cx[-90:-55, -25:15].plot(color="white", edgecolor="black")
# 在ax區域上繪製地圖
gdf.plot(ax=ax, color="red")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed57a460>

基於WTK格式資料

WKT (Well-Known Text) 是一種用於描述地理位置的資料格式。WTK格式的資料包含點、線、多邊形等地理位置資訊。WTK格式的資料可以被許多GIS軟體和地理位置分析工具所讀取和處理。我們可以將帶有WKT資料的DataFrame轉換為GeoDataframe。

df = pd.DataFrame(
    {
        "City": ["Buenos Aires", "Brasilia", "Santiago", "Bogota", "Caracas"],
        "Country": ["Argentina", "Brazil", "Chile", "Colombia", "Venezuela"],
        "Coordinates": [
            "POINT(-58.66 -34.58)",
            "POINT(-47.91 -15.78)",
            "POINT(-70.66 -33.45)",
            "POINT(-74.08 4.60)",
            "POINT(-66.86 10.48)",
        ],
    }
)
df

# 建立新列然後資料轉換
df["Coordinates"] = gpd.GeoSeries.from_wkt(df["Coordinates"])
gdf = gpd.GeoDataFrame(df, geometry="Coordinates", crs="EPSG:4326")

print(gdf.head())

           City    Country                  Coordinates
0  Buenos Aires  Argentina  POINT (-58.66000 -34.58000)
1      Brasilia     Brazil  POINT (-47.91000 -15.78000)
2      Santiago      Chile  POINT (-70.66000 -33.45000)
3        Bogota   Colombia    POINT (-74.08000 4.60000)
4       Caracas  Venezuela   POINT (-66.86000 10.48000)

# 在南美地圖上展示
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
# 定位到南美
ax = world.cx[-90:-55, -25:15].plot(color="white", edgecolor="black")
# 在ax區域上繪製地圖
gdf.plot(ax=ax, color="red")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed445100>

3 新增比例尺

在地理空間資料分析和視覺化過程中，比例尺可以幫助我們瞭解地圖上的距離和大小關係。基於matplotlib進行視覺化時，可以利用matplotlib-scalebar庫新增比例尺。

簡單比例尺

import geopandas as gpd
# pip install matplotlib_scalebar安裝
from matplotlib_scalebar.scalebar import ScaleBar

# nybb為紐約五大地圖
nybb = gpd.read_file(gpd.datasets.get_path("nybb"))
# 北美地區常見座標系統，座標以米為單位
nybb = nybb.to_crs(32619)  
nybb.head()

如下所示，建立ScaleBar物件所需的唯一引數是dx。dx表示輸入圖片每一個畫素代表的長度，units為dx的單位。此引數的值取決於座標參考系的單位。在前面nybb資料集已經使用epsg:32619座標系統，該座標系以單位米為單位，如下所示，可以看到nybb.crs輸出結果中Axis Info項標識了該參考系以metre米為單位。

nybb.crs

<Projected CRS: EPSG:32619>
Name: WGS 84 / UTM zone 19N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 72°W and 66°W, northern hemisphere between equator and 84°N, onshore and offshore. Aruba. Bahamas. Brazil. Canada - New Brunswick (NB); Labrador; Nunavut; Nova Scotia (NS); Quebec. Colombia. Dominican Republic. Greenland. Netherlands Antilles. Puerto Rico. Turks and Caicos Islands. United States. Venezuela.
- bounds: (-72.0, 0.0, -66.0, 84.0)
Coordinate Operation:
- name: UTM zone 19N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

在下面程式碼中新增了比例尺和畫素尺寸，該比例尺採用的是線段式表示方式，即在地圖上繪製一條線段並註明該地圖上該線段所代表的實際距離。

ax = nybb.plot()
# 在地圖中新增比例尺和畫素尺寸
scalebar =ScaleBar(dx=1,units="m")
ax.add_artist(scalebar)

<matplotlib_scalebar.scalebar.ScaleBar at 0x7f75ed426a30>

確定比例尺基準長度

如下所示，以經緯度為單位的epsg:4326座標系，其單位尺度為度（經緯度）。

# nybb為紐約五大區地圖
nybb = gpd.read_file(gpd.datasets.get_path("nybb"))
nybb = nybb.to_crs(4326) 
nybb.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed3a9cd0>

nybb.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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

可以通過計算該地圖中相鄰兩點之間的距離長度來確定比例尺基準長度，要注意的是這兩點應該位於待繪製的地圖中。

from shapely.geometry.point import Point

points = gpd.GeoSeries(
    [Point(-73.9, 40.7), Point(-74.9, 40.7)], crs=4326
)  
# 將兩點轉換到以米為單位的座標系
points = points.to_crs(32619)  
# 計算點之間的距離，距離單位為座標系的單位
distance_meters = points[0].distance(points[1])
distance_meters

84698.53985065906

nybb = nybb.to_crs(4326) 
ax = nybb.plot()
ax.add_artist(ScaleBar(distance_meters,"m"))

<matplotlib_scalebar.scalebar.ScaleBar at 0x7f75ed113400>

比例尺自定義

通過更改 ScaleBar 引數能夠調整比例尺的顯示效果，ScaleBar具體引數如下所示。這些引數的使用可以自行嘗試。

scalebar = ScaleBar(
    dx, # 畫素和長度之間的比例尺。例如，如果一個畫素代表1毫米，則dx=0.001。
    units="m", # 長度單位
    dimension="si-length", # 維度
    label=None, # 刻度尺標籤
    length_fraction=None, # 刻度尺長度佔比
    height_fraction=None, # 刻度尺高度佔比
    width_fraction=None, # 刻度尺寬度佔比
    location=None, # 刻度尺的位置
    pad=None, # 刻度尺和邊框之間的間距
    border_pad=None, # 刻度尺和邊框之間的邊距
    sep=None, # 刻度尺標籤和刻度之間的距離
    frameon=None, # 是否顯示邊框
    color=None, # 刻度尺和標籤的顏色
    box_color=None, # 邊框的顏色
    box_alpha=None, # 邊框的透明度
    scale_loc=None, # 刻度線的位置
    label_loc=None, # 刻度尺標籤的位置
    font_properties=None, # 標籤和刻度線的字型屬性
    label_formatter=None, # 標籤的格式化函數
    scale_formatter=None, # 刻度線的格式化函數
    fixed_value=None, # 固定的數值
    fixed_units=None, # 固定的單位
    animated=False, # 是否允許動畫
    rotation=None, # 刻度尺的旋轉角度
    bbox_to_anchor=None, # bbox的錨點
    bbox_transform=None, # bbox的變換
)

此外，也可以更改一畫素代表的長度單位，如ScaleBar(2, dimension="si-length", units="km")表示圖中1畫素代表實際si-length（國際單位制）中的2km。所支援的長度單位引數如下表所示：

dimension	units
si-length	km, m, cm, um
imperial-length	in, ft, yd, mi
si-length-reciprocal	1/m, 1/cm
angle	deg

一些比例尺引數調整的範例如下

nybb = gpd.read_file(gpd.datasets.get_path("nybb")).to_crs(32619)
ax = nybb.plot()

# 改變位置和方向
scale1 = ScaleBar(
    dx=1,
    label="Scale 1",
    location="lower left",  # 位置
    label_loc="left",
    scale_loc="top",  # 註釋文字相對於橫線方向位置
)

# 改變顏色
scale2 = ScaleBar(
    dx=1,
    label="Scale 2",
    location="center",
    color="#b32400",
    box_color="yellow",
    box_alpha=0.8,  # 透明度
)

# 改變文字
scale3 = ScaleBar(
    dx=1,
    label="Scale 3",
    font_properties={
        "size": "large",
    },  
    location="lower right",  # 位置
    scale_formatter=lambda value, unit: f"> {value} {unit} <",
)


# 改變長度
scale4 = ScaleBar(
    dx=1,
    label="Scale 4",
    length_fraction=0.5, # 表示刻度線佔繪圖區域的50%
    scale_loc="top",
    label_loc="left",
    border_pad=1,
    pad=0.25,
)

ax.add_artist(scale1)
ax.add_artist(scale2)
ax.add_artist(scale3)
ax.add_artist(scale4)

<matplotlib_scalebar.scalebar.ScaleBar at 0x7f753de7d850>

4 圖層操作與幾何運算

4.1 圖層疊加

在geopandas中，overlay()函數是用於將兩個地理圖層進行疊加分析的函數。它可以進行求交集、並集、差集和對稱差集等操作。overlay()函數的基本語法如下：

geopandas.overlay(layer1, layer2, how)

其中，layer1和layer2是兩個geopandas地理圖層物件，how是一個字串，指定要進行的疊加操作。how引數有以下取值：

• intersection：交集
• union：並集
• difference：差集
• symmetric_difference：對稱差集
• identity

在下面的範例中展示overlay函數的使用方式。

準備資料

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon, Point

# 畫一個圓
center = Point(2, 2)  # 圓心座標
radius = 1  # 圓的半徑
circle = center.buffer(radius)
gdf1 = gpd.GeoDataFrame({'geometry': circle, 'circle':[0]})
gdf1.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f74a2e0fc10>

gdf1

# 畫兩個正方形
square = gpd.GeoSeries([Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
                        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])])
gdf2 = gpd.GeoDataFrame({'geometry': square, 'square':[0,1]})
gdf2.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f753db6ab20>

gdf2

# 展示共同繪圖結果
ax = gdf1.plot()
gdf2.plot(ax=ax)

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed2480a0>

交集intersection

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='intersection')
gdf

gdf.plot(cmap="tab10")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed10acd0>

並集union

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='union')
gdf

gdf.plot(cmap="tab10")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed15deb0>

差集difference

# 需要pip install rtree
# 提取在gdf1中，但不在gdf2中的區域
gdf = gpd.overlay(gdf1, gdf2, how='difference')
# 也可以用以下寫法更加直觀
# gdf = gdf1.overlay(gdf2, how='difference')
gdf

gdf.plot(cmap="tab10")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed0d27f0>

對稱差集symmetric_difference

# 需要pip install rtree
# 提取不在gdf1和pdf2交集的區域
gdf = gpd.overlay(gdf1, gdf2, how='symmetric_difference')
gdf

gdf.plot(cmap="tab10")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed259070>

identity

identity是ArcGIS中常用的操作。意思是將源地理圖層與參考圖層進行比較，以在源圖層中標識與參考圖層中相交的區域。使用identity的一個典型場景是當需要分析兩個圖層交集的時候。例如，可能有一個圖層包含了所有的道路，另一個圖層包含了所有的建築。通過使用identity可以找到所有的建築物位於哪些道路上。

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='identity')
gdf

gdf.plot(cmap="tab10")

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed256d60>

4.2 空間連線

空間連線允許將兩個或多個空間資料集合併成一個新的資料集。例如，我們有兩個資料集，一個包含所有城市的邊界，另一個包含所有的人口資料。通過空間連線，我們可以將這兩個資料集合併成一個新的資料集，其中每個城市都會有相應的人口資料。GeoPandas提供sjoin函數將兩個GeoDataFrame資料集基於空間關係進行連線。sjoin函數常用引數如下：

sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuffix='right')

其中，引數含義如下：

• left_df：左側的GeoDataFrame資料集。
• right_df：右側的GeoDataFrame資料集。
• how：連線方式，可選項如下：
  • inner （預設選項）：返回兩個GeoDataFrame中具有共同空間索引的幾何體的交集。
  • left：返回左側GeoDataFrame中的所有幾何體，以及右側GeoDataFrame中與之相交的幾何體。如果右側GeoDataFrame中沒有與左側相交的幾何體，則右側資料中的所有列都將為null。
  • right：與left相反，返回右側GeoDataFrame中的所有幾何體，以及左側GeoDataFrame中與之相交的幾何體。如果左側GeoDataFrame中沒有與右側相交的幾何體，則左側資料中的所有列都將為null。
• predicate：連線的空間關係，常用選項如下：
  • intersects （預設選項）：返回兩個幾何體相交的所有幾何體。
  • contains：返回左側GeoDataFrame中包含於右側GeoDataFrame中的所有幾何體。
  • within：返回右側GeoDataFrame中包含於左側GeoDataFrame中的所有幾何體。
  • touches：返回兩個幾何體相切的所有幾何體。
  • crosses：返回兩個幾何體相交但不相切的所有幾何體。
  • overlaps：返回兩個幾何體部分重疊的所有幾何體。
• lsuffix：組合後左側資料集中幾何物件列的字尾，預設為left。
• rsuffix：組合後右側資料集中幾何物件列的字尾，預設為right。

以下範例展示瞭如何使用sjoin函數進行空間連線。

準備資料

# 建立點 GeoDataFrame
points = gpd.GeoDataFrame(
    [
        {'id': 'p1', 'geometry': Point(0, 0)},
        {'id': 'p2', 'geometry': Point(1, 1)},
        {'id': 'p3', 'geometry': Point(2, 2)},
        {'id': 'p4', 'geometry': Point(3, 3)}
    ],
    crs='EPSG:4326'
)
points

points.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed360bb0>

# 建立多邊形 GeoDataFrame
polygons = gpd.GeoDataFrame(
    [
        {'id': 'P1', 'geometry': Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])},
        {'id': 'P2', 'geometry': Polygon([(1, 1), (1, 3), (3, 3), (3, 1)])}
    ],
    crs='EPSG:4326'
)
polygons

polygons.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed49deb0>

sjoin函數

# 左連線
join_left_df = points.sjoin(polygons, how="left")
# 輸出結果中的每一行都表示左側GeoDataFrame中的一個幾何物件與右側GeoDataFrame中的一個幾何物件進行了連線後得到的結果。
# index_right表示右側GeoDataFrame中的行索引
# id_left：表示左側GeoDataFrame中的幾何物件的ID
# id_right：表示右側GeoDataFrame中的幾何物件的ID
# geometry：表示連線後的幾何物件
join_left_df

# 右連線
join_right_df = points.sjoin(polygons, how="right")
join_right_df

# 設定predicate
join_right_within_df = points.sjoin(polygons, how="left", predicate="contains")
join_right_within_df

4.3 幾何操作

GeoPandas提供了多種用於幾何操作的函數，具體如下：

• 構造方法
  • buffer(distance, resolution=16)：返回一個GeoSeries，其中包含與每個幾何物件距離在給定距離內的所有點的幾何形狀。
  • boundary：返回一個GeoSeries，其中包含每個幾何形狀的集合理論邊界的低維物件。
  • centroid：返回一個GeoSeries，其中包含每個幾何質心的點。
  • convex_hull：返回一個GeoSeries，其中包含表示包含每個物件中所有點的最小凸多邊形的幾何形狀，除非物件中的點數小於三個。對於兩個點，凸包會摺疊成一個線串；對於一個點，凸包是一個點。
  • envelope：返回一個GeoSeries，其中包含包含每個物件的點或最小矩形多邊形（其邊與座標軸平行）的幾何形狀。
  • simplify(tolerance, preserve_topology=True)：返回一個GeoSeries，其中包含每個物件的簡化表示。在geopandas中，simplify函數可以用來簡化多邊形的形狀，以減少地圖資料的大小，同時也可以提高繪圖的效率。當繪圖資料特別大時，該函數很有用。tolerance：簡化容差值，代表簡化幾何物件的形狀後的最大允許誤差。當 tolerance 值越小時，簡化後的幾何物件的形狀越接近原始幾何物件的形狀。preserve_topology：是否保持拓撲結構，預設值為True，表示保持拓撲結構。
  • unary_union：返回一個幾何形狀，其中包含GeoSeries中所有幾何形狀的聯合。
• 幾何變化方法
  • affine_transform(self, matrix)：使用仿射變換矩陣來變換 GeoSeries 的幾何形狀。matrix 為一個包含6、12個元素的列表或元組（2d情況、3d情況）的仿射變換矩陣。關於 matrix 引數的使用需要有仿射變換的知識。
  • rotate(ngle, origin='center', use_radians=False)：旋轉 GeoSeries 的座標系。
  • scale(xfact=1.0, yfact=1.0, zfact=1.0, origin='center')：沿著(x, y, z)三個方向縮放 GeoSeries 的幾何形狀。
  • skew(xs=0.0, ys=0.0, origin='center', use_radians=False)：基於原點origin，沿著 x 和 y 兩個方向傾斜/扭曲 GeoSeries 的幾何形狀。
  • translate(xoff=0.0, yoff=0.0, zoff=0.0)：平移 GeoSeries 的座標系。

構造方法使用範例

import geopandas as gpd

# 載入資料集
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# 展示結果
world.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed587670>

# buffer函數
buffered = world.geometry.buffer(distance=5)

# 顯示結果
buffered.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f74a32e7a30>

# 獲取幾何形狀邊界
boundary = world.geometry.boundary

# 顯示結果
boundary.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f74a3183fd0>

# 獲取幾何質心
centroids = world.geometry.centroid

# 顯示結果
centroids.plot(marker='*', color='green', markersize=5)

<matplotlib.axes._subplots.AxesSubplot at 0x7f74a31560d0>

# 獲取幾何形狀的凸包
convex_hulls = world.geometry.convex_hull

# 顯示結果
convex_hulls.plot(alpha=0.5, edgecolor='k')

<matplotlib.axes._subplots.AxesSubplot at 0x7f75ed012eb0>

# 獲取幾何形狀的外接矩形
envelopes = world.geometry.envelope

# 顯示結果
envelopes.plot(alpha=0.5, edgecolor='k')

<matplotlib.axes._subplots.AxesSubplot at 0x7f753dee4940>

# 對幾何物件進行簡化處理
simplified = world.geometry.simplify(tolerance=0.1)

# 顯示結果
simplified.plot(alpha=0.5, edgecolor='k')

<matplotlib.axes._subplots.AxesSubplot at 0x7f753dd89a60>

merged = world.geometry.unary_union
# 將合併後的幾何物件轉換為GeoDataFrame
gdf_merged = gpd.GeoDataFrame(geometry=[merged])
# 列印後只有一行
print(gdf_merged)
gdf_merged.plot()

                                            geometry
0  MULTIPOLYGON (((-162.440 -79.281, -163.027 -78...





<matplotlib.axes._subplots.AxesSubplot at 0x7f753dd36d60>

幾何變化方法使用範例

# 讀取資料集
import geopandas as gpd

nybb = gpd.read_file(gpd.datasets.get_path('nybb'))
ax = nybb.plot()
# 啟用科學計數法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

from shapely.affinity import affine_transform


# 仿射變換
# 定義仿射變換引數
a, b, d, e, xoff, yoff = 1.5, 0.5, 0, 0.5, 1.5, 0

tmp = nybb.copy()
# 對nybb資料集中的幾何物件進行仿射變換
tmp['geometry'] = tmp['geometry'].apply(lambda x: affine_transform(x, [a, b, d, e, xoff, yoff]))

# 顯示變換後的nybb資料集
tmp.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f753dee4f10>

# 旋轉
nybb_rotate = nybb.geometry.rotate(angle=45)
nybb_rotate.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f753d983d60>

# 縮放
nybb_scale = nybb.geometry.scale(xfact=2, yfact=2, zfact=1)
nybb_scale.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f753b71aee0>

# 傾斜/扭曲
nybb_skew = nybb.geometry.skew(xs=0.1, ys=0.2, use_radians=True)
ax = nybb_skew.plot()
# 啟用科學計數法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

# 平移
nybb_translated = nybb.geometry.translate(xoff=100000, yoff=100000, zoff=0)
ax = nybb_translated.plot()
# 啟用科學計數法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

4.4 彙總

在geopandas中，dissolve函數可以對具有相同屬性值的幾何物件進行合併，從而生成新的幾何物件。在彙總過程中，可以選擇保留某些欄位的資訊，也可以對其他欄位進行統計計算。dissolve函數如下：

geopandas.GeoDataFrame.dissolve(by=None, aggfunc='first', as_index=True, **kwargs)

函數引數介紹：

• by: 可以是一個欄位名，也可以是一列欄位名的列表。表示按照哪些欄位進行彙總。預設為None，即將所有要素合併成一個要素。
• aggfunc: 統計函數，用於對其他欄位進行計算，可以是以下函數之一：
  • 'first': 返回第一個非空值。
  • 'last': 返回最後一個非空值。
  • 'mean': 返回平均值。
  • 'sum': 返回總和。
  • 'min': 返回最小值。
  • 'max': 返回最大值。
  • 自定義函數：可以傳入自定義的聚合函數。
• as_index: 是否將by引數指定的欄位作為行索引，預設為True。
• *kwargs: 其他引數。

下面範例程式碼演示了dissolve函數的使用。

import geopandas as gpd

# 讀取湖北省地圖資料
data = gpd.read_file('https://geo.datav.aliyun.com/areas_v3/bound/420000_full.json')
data.head()

data.plot(cmap='tab20')

<matplotlib.axes._subplots.AxesSubplot at 0x7f753d7ad850>

# 使用dissolve函數合併幾何體，根據地級市的區縣數分組
dissolve_data = data.dissolve(by='childrenNum')
dissolve_data.head()

dissolve_data.plot(cmap='tab20')

<matplotlib.axes._subplots.AxesSubplot at 0x7f753da46ee0>

# 使用dissolve函數合併幾何體，根據地級市的區縣數分組，其他列求均值
dissolve_data = data.dissolve(by='childrenNum', aggfunc='mean')
dissolve_data.head()

4.5 缺失值與空值處理

在使用geopandas處理地理空間資料時，經常會遇到None和Empty這兩個概念。雖然它們都表示缺失值，但它們之間有著一些區別。

• None：表示屬性或者列的值不存在，或者沒有被填充。在geopandas中，如果一個geometry列的值為None，那意味著這個幾何物件不存在。
• Empty：表示屬性或者列的值存在，但是值為空。在geopandas中，如果一個geometry列的值為空，那意味著這個幾何物件是存在的，但是它沒有任何形狀或者座標資訊。

以下為具有一個多邊形、一個缺失值和一個空多邊形的GeoSeries範例：

from shapely.geometry import Polygon
s = gpd.GeoSeries([Polygon([(0, 0), (1, 1), (0, 1)]), None, Polygon([])])
s

0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
1                                                 None
2                                        POLYGON EMPTY
dtype: geometry

在geopandas空間運算中，缺失的幾何圖形通常會傳播。在結果中，這些缺失的幾何圖形也會缺失。另一方面，空的幾何圖形被視為幾何圖形。結果將取決於所進行的運算。如下所示：

s.area

0    0.5
1    NaN
2    0.0
dtype: float64

我們可以通過isna函數和is_empty屬性判斷是否為缺失值或者空值：

# 判斷缺失值
s.isna()

0    False
1     True
2    False
dtype: bool

# 判斷空值
s.is_empty

0    False
1    False
2     True
dtype: bool

# 判斷缺失或為空
s.is_empty | s.isna()

0    False
1     True
2     True
dtype: bool

# 提取既不缺失也不為空的值
s[~(s.is_empty | s.isna())]

0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
dtype: geometry

5 背景地圖疊加

contextily是一個Python庫，它提供了一種簡單的方法將背景地圖（通常是Web瓦片地圖，如OpenStreetMap、Stamen Maps、Mapbox等）新增到地理空間資料視覺化中。使用contextily庫可以使地理空間資料視覺化更加生動、直觀，同時可以提供更多的地理資訊。瓦片地圖是一種基於網格的地圖顯示方式，將地圖劃分為多個小塊，每個小塊稱為「瓦片」，每個瓦片都有自己的座標和編號。這些瓦片可以按需載入，使使用者能夠快速地瀏覽地圖，同時減少了載入時間和資源消耗。瓦片地圖常用於線上地圖應用程式，例如谷歌地圖和百度地圖。
contextily支援使用WGS84 (EPSG:4326)和Spheric Mercator (EPSG:3857)座標系，在Web地圖應用程式中，一般使用EPSG:3857（以米為單位）來顯示瓦片地圖，並使用EPSG:4326（以經緯度為單位）來標記瓦片地圖上的位置。

contextily庫的主要功能包括：

• 從Web地圖提供商獲取地圖圖層
• 將地圖圖層與地理空間資料集合並
• 使用Matplotlib或Bokeh繪製地圖

本文主要介紹contextily簡單使用，contextily具體使用可參考其官方檔案：contextily-doc。contextily庫中基於add_basemap函數在地圖上新增背景地圖。下面是該函數常用可用引數的介紹：

• ax: matplotlib axes物件，用於繪製地圖
• crs: 輸出地圖的座標系，預設為'EPSG:3857'
• source: 底圖的來源，支援多種來源，如OpenStreetMap、Stamen Terrain、Stamen Toner等等，預設為OpenStreetMap
• zoom: 底圖的縮放級別，預設為None，自動根據ax的extent和crs計算。zoom值越高，底圖的縮放級別就越大，地圖顯示的範圍也就越小，細節也會越來越清晰。
• url: 底圖的url地址，預設為None，自動根據source和zoom計算。
• attribution: 底圖的版權資訊，預設為None
• alpha: 底圖的透明度，預設為1.0
• *kwargs: 其他matplotlib.image()函數的可選引數，如cmap、vmin、vmax等等

⚠⚠source引數選擇不同底圖的來源，可能需要大量時間或者特定網路，如果失敗多重試執行程式碼。

5.1 簡單背景地圖疊加

import geopandas as gpd
# 讀取深圳市地圖資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/440300_full.json")
# 簡單繪圖
ax = data.plot(alpha=0.5, edgecolor="k")

# 確定資料所使用的座標系
data.crs
# 將資料集所使用座標系轉為EPSG:3857
data_wm = data.to_crs(epsg=3857)

import contextily as cx
import matplotlib.pyplot as plt 
fig, ax = plt.subplots(figsize=(10, 10))
ax = data_wm.plot(ax = ax, alpha=0.5, edgecolor="k")
# 將自動下載瓦片地圖
cx.add_basemap(ax)

# 儲存地圖
fig.savefig('save.jpg', pad_inches=0.5, bbox_inches='tight', dpi=300)

在上面的程式碼中，如果僅使用經緯度資料疊加瓦片地圖，需要在add_basemap函數中設定crs引數，如下所示：

import contextily as cx
import matplotlib.pyplot as plt 
fig, ax = plt.subplots(figsize=(10, 10))
ax = data.plot(ax = ax, alpha=0.5, edgecolor="k")
# 將自動下載瓦片地圖
cx.add_basemap(ax, crs=data.crs)

可以通過調整zoom引數改變背景瓦片地圖的細節程度，建議zoom值不要過大，下載速度太慢。此外可以通過設定attribution=""去除繪圖水印。

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, zoom=12, attribution="")

5.2 客製化化背景地圖

通過設定add_basemap的source引數能夠指定不同的資料來源，以在地圖上新增不同型別的底圖。如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite)
ax.set_axis_off()

當然也可以疊加多個背景地圖，如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels)

此外也可以疊加多個不同來源的背景圖層。如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.Watercolor, zoom=12)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels, zoom=10)

從上面的案例我們可以看到，contextily通過Provider預置提供商的名稱來獲取相應的Web瓦片地圖。contextily所有預置的地圖提供商可以通過以下cx.providers命令獲取。可以嘗試根據這些提供商客製化瓦片地圖格式。

# cx.providers

除了使用contextily預置的地圖提供商，可以通過source設定給定瓦片地圖地址來指定需要新增的底圖。例如可以新增天地圖，高德地圖，騰訊地圖的瓦片地圖的地址。一些範例的瓦片地圖地址可見：高德谷歌騰訊天地圖地圖瓦片url和在geopandas中疊加線上地圖。
一般地圖服務提供XYZ瓦片地圖連結，其中的xyz代表了地圖的座標系。如下所示：

• x：表示在地圖水平方向上的位置，從左到右遞增，即經度值。
• y：表示在地圖豎直方向上的位置，從上到下遞增，即緯度值。
• z：表示地圖的縮放級別，從0開始遞增，數值越大，地圖顯示的範圍越小，細節越豐富。

在瓦片地圖中，地圖被分成了許多小塊，每個小塊都有一個唯一的編號，也就是xyz座標系。當我們使用地圖服務時，通過改變xyz的值，就可以獲取到不同位置、不同縮放級別下的地圖瓦片，從而達到展示不同地圖的目的。直接通過url設定瓦片地圖範例如下：

fig, ax = plt.subplots(figsize=(10, 10))

ax = data_wm.plot(ax=ax, alpha=0.5, edgecolor="k")

cx.add_basemap(ax, 
                source='https://webst01.is.autonavi.com/appmaptile?style=6&x={x}&y={y}&z={z}',
                zoom=12)

fig.savefig('save.jpg', pad_inches=0, bbox_inches='tight', dpi=300)

fig, ax = plt.subplots(figsize=(10, 10))

ax = data_wm.plot(ax=ax, alpha=0.5, edgecolor="k")

cx.add_basemap(ax, 
                source='https://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
                zoom=12)

fig.savefig('save.jpg', pad_inches=0, bbox_inches='tight', dpi=300)
ax.set_xlim(data_wm.total_bounds[0], data_wm.total_bounds[2])
ax.set_ylim(data_wm.total_bounds[1], data_wm.total_bounds[3])

(2559177.946084248, 2615308.057854809)

5.3 離線背景地圖

在有些時候我們需要離線使用背景瓦片地圖，contextly提供bounds2raster函數用於根據給定的空間範圍和地圖縮放級別，將線上地圖服務中的柵格資料下載為本地檔案。bounds2raster函數常用引數如下：

• w：float型別，表示空間範圍的最小值。
• s：float型別，表示空間範圍的最小值。
• e：float型別，表示空間範圍的最大值。
• n：float型別，表示空間範圍的最大值。
• path：str型別，表示下載的柵格資料檔案的儲存路徑。
• zoom：int或者字串型別，表示地圖縮放級別。如果為字串型別，可以設定為'auto'，表示自動確定最佳的縮放級別。
• source：str型別，表示地圖服務的地址。
• ll：bool型別，表示w、s、e、n是否使用經緯度座標系，預設為False。
• wait：int型別，表示兩次下載之間的等待時間，單位為秒。預設為0。
• max_retries：int型別，表示下載失敗後最大的重試次數。預設為2次。

bounds2raster函數返回RGB影象陣列和瓦片影象邊界框[minX，maxX，minY，maxY]，此外由於網路地圖總是基於WGS84 Web Mercator（EPSG:3857）座標系，因此bounds2raster函數返回和儲存的圖片都是基於EPSG:3857座標系。

import geopandas as gpd
# 讀取鄭州市地圖資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/410100_full.json")
# 簡單繪圖
ax = data.plot(alpha=0.5, edgecolor="k")

# 疊加地圖
ax = data.plot(alpha=0.5, edgecolor="k")
# crs告訴資料集用的座標系統，這裡data.crs為WGS 84(經緯度)
cx.add_basemap(ax,
                crs=data.crs,
                source=cx.providers.Stamen.Watercolor
               )

提取待繪圖區域的邊界資訊

west, south, east, north = bbox = data.total_bounds
bbox

array([112.721178,  34.262109, 114.220962,  34.989506])

根據邊界資訊下載資料

import contextily as cx
import matplotlib.pyplot as plt 
img, ext = cx.bounds2raster(west,
                             south,
                             east,
                             north,
                             "demo.tif",
                             source=cx.providers.Stamen.Watercolor,
                             ll=True
                            )

# 展示下載的資料
plt.axis('off')
plt.imshow(img)
# 邊界範圍
ext

(12523442.714243276, 12719121.506653327, 4030983.1236470537, 4187526.157575096)

有了背景地圖，add_basemap中的source函數設定檔案路徑地址就可以離線疊加地圖。

ax = data.plot(alpha=0.5, edgecolor="k")
# crs告訴資料集用的座標系統，這裡data.crs為WGS 84(經緯度)
cx.add_basemap(ax,
                crs=data.crs,
                source="demo.tif",
               )

6 參考

• GeoPandas
• GeoPandas-doc
• GeoPandas Examples Gallery
• Python繪製資料地圖1-GeoPandas入門指北
• Python繪製資料地圖2-GeoPandas地圖視覺化
• matplotlib-scalebar
• contextily
• contextily-doc
• 高德谷歌騰訊天地圖地圖瓦片url
• 在geopandas中疊加線上地圖