geopandas：Python绘制数据地图

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__

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

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

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

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>

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>

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

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>

# 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(
    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的变换
)

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>

geopandas.overlay(layer1, layer2, how)

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>

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

gdf.plot(cmap="tab10")

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

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

gdf.plot(cmap="tab10")

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

# 需要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>

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

gdf.plot(cmap="tab10")

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

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

gdf.plot(cmap="tab10")

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

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

# 创建点 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>

# 左连接
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

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

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

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

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

s.area

0    0.5
1    NaN
2    0.0
dtype: float64

# 判断缺失值
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

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)

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)

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

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)

# cx.providers

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)

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)

ax = data.plot(alpha=0.5, edgecolor="k")
# crs告诉数据集用的坐标系统，这里data.crs为WGS 84(经纬度)
cx.add_basemap(ax,
                crs=data.crs,
                source="demo.tif",
               )