使用geo_plotkit高效绘制连续型栅格：Python实现中国地形图可视化

geo_plotkit的初衷，是对繁杂的Python地理绘图相关步骤进行高度封装，目标是仅用几句代码就能轻松绘制出达到顶刊水准的配图。这一设计理念的核心在于效率，对于科研与生产而言，效率的提升至关重要，并且这种效率优势会随着使用而不断累积放大。

本文将介绍如何利用geo_plotkit对如DEM、NDVI等连续型栅格数据进行快速可视化，并结合frykit包，以案例形式演示如何轻松实现我国任意行政区域的地形可视化。

重要说明：geo_plotkit是本地开发的绘图包，暂未上传至PyPI（因此无法通过pip install安装），此系列文章仅作为技术使用记录与分享。

使用说明

查看和获取DEM数据

首先，我们需要了解绘图包中包含了哪些DEM数据，并获取其文件路径。

from geo_plotkit.data import get_dem_path, list_dems

help(get_dem_path)

print('当前内置的DEM的tif数据有：',list_dems())
print('当前使用的DEM数据为：',get_dem_path('China_dem_1km'))

图1：get_dem_path函数的帮助文档，展示了其用法和内置的DEM数据列表。

efficient_tif_plot用法

efficient_tif_plot函数是快速实现连续栅格可视化的核心，它能绘制栅格并自动添加具有ArcGIS风格的可调色条。

import geo_plotkit.plot as gp

help(gp.efficient_tif_plot)

Help on function efficient_tif_plot in module geo_plotkit.plot.quick_plot:

efficient_tif_plot(ax, tif_file, cmap='terrain', loc=[0.9, 0.05, 0.02, 0.1], title='Elevation(m)', cmap_clip=None, clip_polygon=None, titlefontsize=12, titleloc='center', labelsize=10, cbar_ec='gray', pad=5)
    绘制单个 GeoTIFF 栅格（例如 NDVI 高程数据）并在指定位置添加可自定义的颜色条

    该函数特别适用于地理可视化（如高程、降水等空间分布），
    在给定的 Matplotlib Axes（通常带有 Cartopy 投影）上显示栅格影像，
    并附带一个可控制位置、颜色范围、裁剪区域和外观的颜色条（Colorbar）。

    ----------
    参数：
    ax : matplotlib.axes.Axes
        要绘制到的坐标轴对象，建议为带有地理投影的 Cartopy Axes（如 PlateCarree）。
    tif_file : str
        栅格文件路径（GeoTIFF 格式），建议为 WGS84 坐标系。
    cmap : str 或 matplotlib.colors.Colormap, 默认 'terrain'
        颜色映射，可以是 Matplotlib 内置 colormap 名称，也可以是 `cmaps` 模块提供的自定义色系。
    loc : list[float], 默认 [.9, .05, .02, .1]
        颜色条位置及大小，格式为 `[x, y, 宽度, 高度]`，相对于坐标轴的比例。
        例如 `[.9, .05, .02, .1]` 表示在右下角放一个竖直颜色条。
    title : str, 默认 'Elevation(m)'
        颜色条标题（用于说明数据单位或含义）。
    cmap_clip : tuple[float, float], 可选
        颜色映射裁剪范围（0~1 之间的比例值），用于增强视觉对比度。
        例如 `(0.1, 0.98)` 表示从 colormap 的 10% 到 98% 范围取色。
    clip_polygon : list[shapely.geometry.Polygon], 可选
        裁剪区域多边形列表（EPSG:4326 坐标系），如果提供则仅显示该范围内的数据。
    titlefontsize : int, 默认 12
        颜色条标题字体大小。
    titleloc : {'left', 'center', 'right'}, 默认 'center'
        颜色条标题位置。
    labelsize : int, 默认 10
        颜色条刻度标签的字体大小。
    cbar_ec : str, 默认 'gray'
        颜色条边框颜色（edgecolor）。
    pad : float, 默认 5
        颜色条标题与刻度轴之间的间距。

    ----------
    返回：
    ax : matplotlib.axes.Axes
        绘制完成的坐标轴对象。
    cbar : matplotlib.colorbar.Colorbar
        创建的颜色条对象，可用于进一步调整格式。

    ----------
    注意事项：
    1. `tif_file` 必须是地理参考的 GeoTIFF（WGS84 坐标系），否则显示可能错位。
    2. 如果需要绘制某个行政区或不规则区域，应先准备好 `clip_polygon`（如 shapefile 转换成 shapely Polygon）。
    3. `loc` 控制颜色条位置时，需要根据绘图布局自行调整，防止遮挡主体地图。
    4. `cmap_clip` 可以在数据分布极不均匀时提高视觉效果（避免高值/低值过度主导色彩）。

    ----------
    用例：
    >>> fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.PlateCarree()})
    >>> china = fshp.get_cn_border()
    >>> efficient_tif_plot(
    ...     ax, 'China_dem_1km.tif',
    ...     cmap='terrain',
    ...     loc=[.88, .05, .03, .12],
    ...     title='Elevation(m)',
    ...     cmap_clip=(0.1, 0.98),
    ...     clip_polygon=[china],
    ...     titlefontsize=9,
    ...     titleloc='center',
    ...     labelsize=8,
    ...     cbar_ec='k',
    ...     pad=6
    ... )

搭配 fshp 实现任意行政区地形可视化

技巧：轻松获取任意行政区的经纬度边界范围

借助 frykit.shp（简称 fshp）可以非常方便地获取各级行政区的矢量边界。

import frykit.shp as fshp

fshp.get_cn_province('湖北省').bounds

图2：使用 fshp 获取湖北省的边界坐标范围。

案例一：绘制全国地形图

使用 fshp.get_cn_border() 获取全国边界，并作为 clip_polygon 参数传入（注意需要放在列表中）。

import geo_plotkit.plot as gp
from   geo_plotkit import bold_border, set_tick_style, set_journal_style, set_map_ticks
from geo_plotkit.data import get_dem_path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import frykit.plot as fplt
import frykit.shp as fshp

set_journal_style('ieee',fontsize=13)
fig, ax = plt.subplots(figsize=(12,9),subplot_kw={'projection': ccrs.PlateCarree()})
extent = [73.5, 135.5, 12, 54]
set_map_ticks(ax, extent, grid=True, display_direction=['top','left'])
im = gp.add_basemap(ax, fname='eAtlas', clip_bounds=extent)
gp.add_admin(ax, level="国界", linewidth=1)

bold_border(ax, 2)
set_tick_style(ax, major_length=6, major_width=1.2, direction='out')

ax, cbar = gp.efficient_tif_plot(
    ax,
    get_dem_path('China_dem_1km'),
    cmap='terrain',
    clip_polygon=[fshp.get_cn_border()]
)

图3：使用 geo_plotkit 结合全国边界绘制的中国地形图。

案例二：绘制四川省地形图

方法与全国类似，只需将 clip_polygon 参数替换为四川省的边界。

import geo_plotkit.plot as gp
from   geo_plotkit import bold_border, set_tick_style, set_journal_style, set_map_ticks
from geo_plotkit.data import get_dem_path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import frykit.plot as fplt
import frykit.shp as fshp

set_journal_style('ieee',fontsize=13)
fig, ax = plt.subplots(figsize=(12,9),subplot_kw={'projection': ccrs.PlateCarree()})
extent = [97, 109, 26, 34.5]
set_map_ticks(ax, extent, grid=True, display_direction=['top','left'])
im = gp.add_basemap(ax, fname='clean', clip_bounds=extent)

bold_border(ax, 2)
set_tick_style(ax, major_length=6, major_width=1.2, direction='out')

ax, cbar = gp.efficient_tif_plot(
    ax,
    get_dem_path('China_dem_1km'),
    cmap='OceanLakeLandSnow',
    clip_polygon=[fshp.get_cn_province('四川省')],
    title='DEM(m)',
    titlefontsize=15,
    loc=[.9, .05, .03, .12],
    labelsize=12,
    pad=8
)
fplt.add_cn_province(ax, '四川省', lw=1.2)

图4：四川省的地形可视化效果图，使用了不同的配色方案。

案例三：绘制多省份合并区域地形图

有时我们需要绘制多个省份合并区域的地形。首先，可以使用 shapely.ops.unary_union 来合并多个矢量并获取其外接边界范围。

from shapely.ops import unary_union

unary_union(fshp.get_cn_province(['湖北省','湖南省','江西省'])).bounds
# (108.36248799999998, 24.48894, 118.482124, 33.277547999999996)

当 clip_polygon 参数直接传入一个由多个多边形组成的列表（如 fshp.get_cn_province([...]) 的返回结果）时，无需再额外嵌套成列表。

import geo_plotkit.plot as gp
from   geo_plotkit import bold_border, set_tick_style, set_journal_style, set_map_ticks
from geo_plotkit.data import get_dem_path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import frykit.plot as fplt
import frykit.shp as fshp

set_journal_style('ieee',fontsize=13)
fig, ax = plt.subplots(figsize=(12,9),subplot_kw={'projection': ccrs.PlateCarree()})
extent = [108, 120, 24, 34]
set_map_ticks(ax, extent, grid=True, display_direction=['top','left'])
im = gp.add_basemap(ax, fname='clean', clip_bounds=extent)
gp.add_admin(ax, level="省", linewidth=1)

bold_border(ax, 2)
set_tick_style(ax, major_length=6, major_width=1.2, direction='out')

ax, cbar = gp.efficient_tif_plot(
    ax,
    get_dem_path('China_dem_1km'),
    cmap='OceanLakeLandSnow',
    cmap_clip=(0.1,0.9), # 切片色系
    clip_polygon=fshp.get_cn_province(['湖北省','湖南省','江西省']),
    title='DEM(m)',
    titlefontsize=15,
    loc=[.85, .05, .03, .12],
    labelsize=12,
    pad=8
)

图5：湖北、湖南、江西三省合并区域的地形图。

通用参数说明

除了搭配 frykit 包使用，efficient_tif_plot 函数也完全支持自定义的栅格数据和矢量边界。例如，你可以用它来绘制NDVI栅格图，或者裁剪显示某个特定流域的范围。

1. tif_file：直接传入你的栅格文件路径（要求为WGS84坐标系，.tif后缀）。
2. clip_polygon：可以使用 geopandas 读取Shapefile，并通过 .geometry 属性获取矢量边界后传入。

希望通过本文，你能掌握使用 geo_plotkit 进行高效地理可视化的方法。如果你对 Python 地理数据处理和绘图有更多兴趣，欢迎在云栈社区进行更深入的交流与探讨。

---