使用 Python 绘制填充地图、栅格地图 + 等降水量线

最近有小伙伴问到在栅格地图上绘制等降水量线的问题。本课程将 R 版本的代码翻译成 Python，使用 Python 生态中的 xarray、rasterio、geopandas、matplotlib 等库，实现与 R 版本完全等价的三张地图：

1. 各省份年均降水量填充地图（分段色彩）
2. 降水量栅格地图 + 等降水量线（连续渐变色）
3. 栅格地图 + 带文字标签的等降水量线

附件中提供了 2021 年各月的降水量栅格数据：pre_2021.nc，数据来源于 中国1km分辨率逐月降水量数据集（1901-2023）：https://data.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2 。

在学习本课程之前，建议预先了解：

使用 Python 绘制历年中国省市区县地图（小地图版本+长版）：https://rstata.duanshu.com/#/brief/course/72c36d788f9d4497a7e33f23f9437598

---

一、R 语言与 Python 的包对应关系

绘制降水量地图涉及的核心 R 包及其 Python 对应关系：

R 包	Python 对应	主要用途
terra / raster	rasterio + xarray	读取、裁剪、处理栅格数据
sf	geopandas + shapely	矢量数据读取和处理
ggplot2 + ggspatial	matplotlib + geopandas	地图绘制
nngeo::st_remove_holes	shapely exterior 提取	去除多边形孔洞
rasterToContour	matplotlib.contour	生成等值线
rasterToPolygons	scipy.ndimage + shapely.geometry.box	栅格转多边形
geomtextpath	matplotlib ax.text + 旋转角度	等值线文字标签

---

二、使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错。因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径。

2.1 安装 reticulate（仅首次）

# 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像）
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

# 仅在尚未安装时才安装
if (!requireNamespace("reticulate", quietly = TRUE)) {
  install.packages("reticulate")
  message("reticulate 安装完成！")
} else {
  message("reticulate 已安装，版本：", packageVersion("reticulate"))
}

2.2 虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

library(reticulate)

.venv_name   <- ".venv"
.venv_python <- virtualenv_python(.venv_name)

# 虚拟环境不存在时自动创建
if (!file.exists(.venv_python)) {
  virtualenv_create(.venv_name)
  .venv_python <- virtualenv_python(.venv_name)
}

# 通过环境变量抢先锁定 Python（优先级最高，早于任何 {python} chunk）
Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

2.3 安装 Python 包（仅首次）

本课程需要以下 Python 包：

• xarray：读取 NetCDF 格式的气象数据
• rasterio：栅格数据的裁剪、掩膜操作
• geopandas：矢量地理数据处理
• shapely：几何对象操作
• numpy / pandas：数值计算和表格处理
• matplotlib：地图绘制
• scipy：栅格降采样（zoom 函数）
• netCDF4：NetCDF 文件读取后端

py_pkgs <- c(
  "numpy", "pandas", "geopandas", "shapely",
  "rasterio", "xarray", "netCDF4", "scipy", "matplotlib"
)
installed <- py_list_packages(".venv")$package
need_install <- setdiff(py_pkgs, installed)

if (length(need_install) > 0) {
  message("正在安装缺失的包: ", paste(need_install, collapse = ", "))
  virtualenv_install(".venv", packages = need_install)
  message("安装完成！")
} else {
  message("所有 Python 包已就绪，无需安装")
}

2.4 验证激活状态

py_config()

2.5 查看已安装的包

pkgs <- py_list_packages(".venv")
key_pkgs <- c("numpy", "pandas", "geopandas", "rasterio",
              "xarray", "netCDF4", "scipy", "matplotlib")
pkgs[pkgs$package %in% key_pkgs, c("package", "version")]

2.6 虚拟环境管理常用命令

# 查看所有已创建的虚拟环境
virtualenv_list()

# 删除虚拟环境（当不再需要时）
# virtualenv_remove(".venv")

# 升级某个包
# virtualenv_install(".venv", packages = "geopandas", ignore_installed = TRUE)

---

三、全局配置与数据读取

3.1 全局参数配置

import os
import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.colors as mcolors
import matplotlib.font_manager as fm
from matplotlib.patches import Rectangle
from matplotlib.colors import LinearSegmentedColormap
from shapely.geometry import box, mapping
from shapely.ops import unary_union

# 工作目录（与 Rmd 文件所在目录一致）
BASE_DIR = os.getcwd()
print(f"工作目录: {BASE_DIR}")

# 中国地图常用等积圆锥投影（Albers，单位：米）
# 与 R 版本完全一致
MYCRS = "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 小地图（南海诸岛）在 Albers 投影下的 bbox
# 等价于 R 的 small_bbox <- st_bbox(c(xmin=120000, xmax=1766004.1, ymax=2557786.0, ymin=320000))
SMALL_BBOX = box(120000, 320000, 1766004.1, 2557786.0)

# 主图绘制范围（Albers）
# 等价于 R 的 plotbbox <- st_bbox(c(xmin=-2725586, xmax=2982768, ymax=6000000, ymin=1800655))
PLOT_BBOX = box(-2725586, 1800655, 2982768, 6000000)
XMIN, YMIN, XMAX, YMAX = -2725586, 1800655, 2982768, 6000000

print("全局参数配置完成")
print(f"坐标系: {MYCRS[:50]}...")

3.2 中文字体设置

def setup_chinese_font():
    """尝试设置中文字体，按优先级依次尝试"""
    for fname in ["LXGWWenKai-Regular.ttf", "SimHei.ttf"]:
        fpath = os.path.join(BASE_DIR, fname)
        if os.path.exists(fpath):
            fm.fontManager.addfont(fpath)
            prop = fm.FontProperties(fname=fpath)
            plt.rcParams["font.family"] = prop.get_name()
            plt.rcParams["axes.unicode_minus"] = False
            print(f"  使用字体文件: {fpath}")
            return prop

    for sys_font in ["Heiti TC", "PingFang SC", "Noto Sans CJK SC", "STHeiti", "SimHei"]:
        try:
            prop = fm.FontProperties(family=sys_font)
            plt.rcParams["font.family"] = sys_font
            plt.rcParams["axes.unicode_minus"] = False
            print(f"  使用系统字体: {sys_font}")
            return prop
        except Exception:
            continue

    plt.rcParams["axes.unicode_minus"] = False
    print("  使用默认字体（中文可能无法正常显示）")
    return fm.FontProperties()

CNFONT = setup_chinese_font()

3.3 读取降水量栅格数据

降水量数据为 NetCDF 格式，包含 12 个月的月均降水量。我们需要计算年均值。

R 代码等价操作：rast("pre_2021.nc") -> rst; app(rst, mean, na.rm=T) -> rst

from rasterio.transform import from_bounds

def load_precipitation_raster(nc_path):
    """
    读取 NetCDF 降水量数据，计算年均值
    等价于 R: rast("pre_2021.nc") |> app(mean, na.rm=T)
    """
    ds = xr.open_dataset(nc_path, engine="netcdf4")
    print(f"数据集变量: {list(ds.data_vars)}")
    print(f"数据集坐标: {list(ds.coords)}")

    # 找到降水量变量
    var_name = [v for v in ds.data_vars if "pre" in v.lower() or "prec" in v.lower()]
    if not var_name:
        var_name = list(ds.data_vars)
    var_name = var_name[0]
    print(f"降水量变量名: {var_name}, shape: {ds[var_name].shape}")

    da = ds[var_name]
    da_float = da.astype(float)

    # 处理 nodata 值（R 的 terra 自动处理）
    nodata_val = da.encoding.get("_FillValue", None) or da.attrs.get("_FillValue", None)
    if nodata_val is not None:
        da_float = da_float.where(da != nodata_val)
        print(f"nodata 值: {nodata_val}")

    # 计算年均值（沿时间维度，等价于 app(rst, mean, na.rm=T)）
    if da_float.ndim == 3:
        annual_mean = da_float.mean(dim=da_float.dims[0], skipna=True).values
    else:
        annual_mean = da_float.values

    # 获取经纬度坐标
    lat = ds["lat"].values if "lat" in ds else ds["latitude"].values
    lon = ds["lon"].values if "lon" in ds else ds["longitude"].values

    # 创建 rasterio affine transform
    lon_min, lon_max = float(lon.min()), float(lon.max())
    lat_min, lat_max = float(lat.min()), float(lat.max())
    height, width = annual_mean.shape
    transform = from_bounds(lon_min, lat_min, lon_max, lat_max, width, height)

    ds.close()
    print(f"年均降水量范围: {np.nanmin(annual_mean):.1f} ~ {np.nanmax(annual_mean):.1f} mm")
    return annual_mean, transform, "EPSG:4326"

annual_mean, rst_transform, rst_crs = load_precipitation_raster("pre_2021.nc")
print(f"栅格大小: {annual_mean.shape}")

3.4 裁剪到中国范围

R 代码等价操作：st_union(prov) |> nngeo::st_remove_holes() -> cn; terra::crop(vect(cn)) |> terra::mask(vect(cn))

from rasterio.mask import mask as rasterio_mask
from rasterio.io import MemoryFile
import rasterio.crs as rcrs
import rasterio.features

def clip_raster_to_china(data, transform, crs_str, prov_shp_path):
    """
    将栅格裁剪到中国范围（合并省级行政区，去除内部孔洞）
    """
    prov = gpd.read_file(prov_shp_path)
    cn = unary_union(prov.geometry)

    # 去除内部孔洞（等价于 nngeo::st_remove_holes）
    if cn.geom_type == "Polygon":
        from shapely.geometry import Polygon
        cn = Polygon(cn.exterior)
    elif cn.geom_type == "MultiPolygon":
        from shapely.geometry import MultiPolygon, Polygon
        cn = MultiPolygon([Polygon(p.exterior) for p in cn.geoms])

    height, width = data.shape
    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "width": width,
        "height": height,
        "count": 1,
        "crs": rcrs.CRS.from_string(crs_str),
        "transform": transform,
        "nodata": float("nan"),
    }

    with MemoryFile() as memfile:
        with memfile.open(**profile) as dataset:
            dataset.write(data.astype("float32"), 1)
        with memfile.open() as dataset:
            masked_data, masked_transform = rasterio_mask(
                dataset, [mapping(cn)], crop=True, all_touched=True, nodata=float("nan")
            )

    result = masked_data[0].astype(float)
    result[result == float("nan")] = np.nan
    print(f"裁剪后栅格大小: {result.shape}")
    return result, masked_transform, rcrs.CRS.from_string(crs_str)

data_cn, transform_cn, crs_cn = clip_raster_to_china(
    annual_mean, rst_transform, rst_crs, "2021行政区划/省.shp"
)

---

四、分区域提取均值（各省降水量）

R 代码等价操作：terra::extract(rst, vect(prov), fun=mean, na.rm=T)

此处代码需下载讲义材料查看～

年均降水量直方图

# 等价于 R 的 hist(provmap2$年度降水量_mm)
fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(df_prov["年度降水量_mm"].dropna(), bins=15, color="#3288bd",
        edgecolor="white", linewidth=0.5)
ax.set_xlabel("年均降水量 (mm)", fontproperties=CNFONT)
ax.set_ylabel("省份数量", fontproperties=CNFONT)
ax.set_title("各省份年均降水量分布", fontproperties=CNFONT)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout()
plt.show()

五、读取地图数据

# 读取小地图版本的中国省级地图
# 等价于 R 的 read_sf("chinaprov2019mini/chinaprov2019mini.shp")
provmap = gpd.read_file("chinaprov2019mini/chinaprov2019mini.shp")
provmap = provmap[provmap["省代码"].notna()].copy()
print(f"省级地图要素数: {len(provmap)}")
print(provmap.columns.tolist())

# 线条要素：九段线、海岸线、小地图框格
# 等价于 R 的 read_sf("chinaprov2019mini/chinaprov2019mini_line.shp")
linemap = gpd.read_file("chinaprov2019mini/chinaprov2019mini_line.shp")
linemap = linemap[linemap["class"].isin(["九段线", "海岸线", "小地图框格"])].copy()
print(f"线条要素数: {len(linemap)}")

---

六、绘制填充地图（pic1）

各省份年均降水量填充地图，等价于 R 版本的 pic1.png。

6.1 合并省级数据与地图

# 找合并键
def find_join_key(left_df, right_df):
    left_cols = set(left_df.columns) - {"geometry"}
    right_cols = set(right_df.columns) - {"geometry"}
    for c in ["省", "省代码", "省份", "NAME", "PNAME"]:
        if c in left_cols and c in right_cols:
            return c
    common = left_cols & right_cols
    return next(iter(common)) if common else None

join_key = find_join_key(provmap, df_prov)
print(f"合并键: {join_key}")

provmap2 = provmap.merge(df_prov, how="left", on=join_key)
print(f"合并后: {len(provmap2)} 行")
print(provmap2[["省", "年度降水量_mm"]].head(10))

6.2 数据分组

# 等价于 R 的 cut(年度降水量_mm, breaks=seq(-400,2400,by=400))
def classify_precip(val):
    if pd.isna(val):
        return "No data"
    elif val <= 0:
        return "No data"
    elif val <= 400:
        return "0 ~ 400"
    elif val <= 800:
        return "400 ~ 800"
    elif val <= 1200:
        return "800 ~ 1200"
    elif val <= 1600:
        return "1200 ~ 1600"
    elif val <= 2000:
        return "1600 ~ 2000"
    else:
        return "2000 ~ 2400"

provmap2["group"] = provmap2["年度降水量_mm"].apply(classify_precip)
print(provmap2["group"].value_counts())

6.3 绘图

此处代码需下载讲义材料查看～

七、栅格地图 + 等降水量线

7.1 将栅格重投影到 Albers（关键步骤）

原始栅格是 WGS84 坐标，绘制等降水量线时需要与地图投影（Albers）一致。正确做法是先把栅格重投影到 Albers，然后直接在地图坐标系的 axes 上调用 ax.contour 绘制等值线，无需生成 GeoDataFrame 中转。

为什么之前的做法（WGS84 → GeoDataFrame → 转 Albers → plot）效果不好？
裁剪后的栅格边界有大量 NaN，contour 会把连续等值线切成成千上万条碎片线段，这些碎片拼接后虽然位置正确，但密集的重叠会让线条看起来很粗且凌乱。ax.clabel 也无法正确找到适合放标签的位置。
直接在重投影栅格上调用 contour，则等值线是连续的，标签也能自动沿线排布。

等价于 R 中对 raster 对象做 projectRaster 再 rasterToContour。

from rasterio.warp import calculate_default_transform, reproject, Resampling
import rasterio.transform

def reproject_raster_to_albers(data, transform, src_crs_str="EPSG:4326"):
    """
    将 WGS84 栅格重投影到 Albers 等积圆锥投影。
    等价于 R 的 projectRaster(rst, crs=mycrs)
    返回 (data_albers, albers_transform, x_arr, y_arr)
    """
    h, w = data.shape
    src_crs = rasterio.crs.CRS.from_string(src_crs_str)
    dst_crs = rasterio.crs.CRS.from_string(MYCRS)

    bounds = rasterio.transform.array_bounds(h, w, transform)
    dst_transform, dst_width, dst_height = calculate_default_transform(
        src_crs, dst_crs, w, h, *bounds
    )

    data_albers = np.full((dst_height, dst_width), np.nan, dtype=np.float32)

    with MemoryFile() as memfile:
        with memfile.open(driver="GTiff", dtype="float32",
                          width=w, height=h, count=1,
                          crs=src_crs, transform=transform, nodata=np.nan) as src_ds:
            src_ds.write(data.astype("float32"), 1)
        with memfile.open() as src_ds:
            reproject(
                source=rasterio.band(src_ds, 1),
                destination=data_albers,
                src_transform=transform, src_crs=src_crs,
                dst_transform=dst_transform, dst_crs=dst_crs,
                resampling=Resampling.bilinear,
                src_nodata=np.nan, dst_nodata=np.nan,
            )

    # 生成像元中心坐标（x 递增，y 递减=北到南）
    x_arr = np.array([dst_transform.c + dst_transform.a * (j + 0.5)
                      for j in range(dst_width)])
    y_arr = np.array([dst_transform.f + dst_transform.e * (i + 0.5)
                      for i in range(dst_height)])
    return data_albers, dst_transform, x_arr, y_arr

print("正在重投影栅格到 Albers...")
data_albers, albers_transform, albers_x, albers_y = reproject_raster_to_albers(
    data_cn, transform_cn
)
print(f"Albers 栅格大小: {data_albers.shape}")
print(f"值范围: {np.nanmin(data_albers):.1f} ~ {np.nanmax(data_albers):.1f} mm")

7.2 栅格数据转多边形（降采样）

等价于 R 的 rst |> aggregate(fact=10) |> rasterToPolygons() |> st_as_sf()

from scipy.ndimage import zoom

def raster_to_polygons(data, transform, resample_factor=10):
    """
    将栅格数据降采样后转成多边形
    等价于 raster::aggregate(fact=10) + rasterToPolygons()
    """
    h, w = data.shape
    new_h = max(1, h // resample_factor)
    new_w = max(1, w // resample_factor)

    # 降采样（等价于 aggregate(fact=10)）
    data_f = data.astype(float)
    scale_h = new_h / h
    scale_w = new_w / w
    data_small = zoom(data_f, (scale_h, scale_w), order=1, prefilter=False)

    # 新的 transform
    lon_min = transform.c
    lat_max = transform.f
    lon_max = transform.c + transform.a * w
    lat_min = transform.f + transform.e * h
    new_transform = from_bounds(lon_min, lat_min, lon_max, lat_max, new_w, new_h)

    # 每个像元转成一个矩形多边形
    rows_list = []
    for i in range(new_h):
        for j in range(new_w):
            val = data_small[i, j]
            if np.isnan(val):
                continue
            x0 = new_transform.c + j * new_transform.a
            y1 = new_transform.f + i * new_transform.e
            x1 = x0 + new_transform.a
            y0 = y1 + new_transform.e
            geom = box(min(x0, x1), min(y0, y1), max(x0, x1), max(y0, y1))
            rows_list.append({"mean": val, "geometry": geom})

    gdf = gpd.GeoDataFrame(rows_list, crs="EPSG:4326")
    return gdf

print("正在生成栅格多边形（约需1分钟）...")
rstpoints = raster_to_polygons(data_cn, transform_cn, resample_factor=10)
print(f"多边形数量: {len(rstpoints)}")

# 转到 Albers 投影
rstpoints = rstpoints.to_crs(MYCRS)

7.3 分离主图和小地图区域

等价于 R 中用 st_intersection 分离南海小地图和主图的操作：

# 小地图内的多边形
small_bbox_gdf = gpd.GeoDataFrame(geometry=[SMALL_BBOX], crs=MYCRS)
rstpoints_small = gpd.overlay(rstpoints, small_bbox_gdf, how="intersection")
print(f"小地图多边形数: {len(rstpoints_small)}")

# 主图内的多边形
plot_bbox_gdf = gpd.GeoDataFrame(geometry=[PLOT_BBOX], crs=MYCRS)
rstpoints_main = gpd.overlay(rstpoints, plot_bbox_gdf, how="intersection")
print(f"主图多边形数: {len(rstpoints_main)}")

7.4 移动小地图到恰当位置

等价于 R 的 mutate(geometry = geometry * 0.5 + c(2100000, 1665139))

from shapely.affinity import scale as affine_scale, translate

def scale_translate_geom(geom, scale=0.5, tx=2100000, ty=1665139):
    """
    等价于 R 的 geometry * 0.5 + c(2100000, 1665139)
    先以原点缩放 0.5 倍，再平移
    """
    scaled = affine_scale(geom, xfact=scale, yfact=scale, origin=(0, 0))
    return translate(scaled, xoff=tx, yoff=ty)

rstpoints_small = rstpoints_small.copy()
rstpoints_small["geometry"] = rstpoints_small["geometry"].apply(scale_translate_geom)
rstpoints_small = rstpoints_small.set_crs(MYCRS, allow_override=True)

# 合并主图和小地图的多边形
rstpoints_all = pd.concat([rstpoints_main, rstpoints_small], ignore_index=True)
rstpoints_all = gpd.GeoDataFrame(rstpoints_all, crs=MYCRS)
print(f"合并后多边形数: {len(rstpoints_all)}")

7.5 绘制栅格地图 + 等降水量线（pic2）

等降水量线直接在 Albers 坐标系的 axes 上调用 ax.contour 绘制，与地图坐标系完全对齐。
这等价于 R 的 geom_sf(data=line, color="white", linewidth=0.5)。

此处代码需下载讲义材料查看～

这样我们就把等降水量线添加到了地图上。

---

八、添加等降水量线文字标签（pic3）

R 版本使用 geomtextpath::geom_textsf 在线上自动添加文字标签。
Python 中使用 ax.clabel 实现相同效果——它会在等值线的合适位置插入标签，并将线切断，防止文字与线条重叠。

# 等价于 R 的:
# line %>% st_cast("LINESTRING") %>% mutate(level = paste0(level, " mm")) -> line2
# geom_textsf(data=line2, aes(label=level), linecolour="white", color="gray30", ...)

fig, ax = plt.subplots(figsize=(10, 8.5))
fig.patch.set_facecolor("white")

# 绘制栅格多边形
rstpoints_all.plot(ax=ax, column="mean", cmap=cmap, norm=norm,
                   linewidth=0, edgecolor="none", zorder=1)

# 绘制省级轮廓
provmap.plot(ax=ax, facecolor="none", edgecolor="gray",
             linewidth=0.03, zorder=3)

# 绘制线条
for cls, style in LINE_STYLES.items():
    subset = linemap[linemap["class"] == cls]
    if len(subset) > 0:
        subset.plot(ax=ax, color=style["color"],
                    linewidth=style["linewidth"], zorder=5)

# 等降水量线 + 文字标签（等价于 geomtextpath::geom_textsf）
# 同样直接在 Albers 坐标系中绘制
cs3 = ax.contour(
    albers_x,
    albers_y[::-1],
    np.flipud(data_albers),
    levels=contour_levels,
    colors="white",
    linewidths=0.5,
    zorder=6,
    alpha=0.9,
)

# clabel：只标注 4 个 level，用 FontProperties 控制字体大小
labeled_levels = [800, 1200, 1600, 2000]
labels = ax.clabel(
    cs3,
    levels=labeled_levels,
    inline=True,
    inline_spacing=2,
    fmt=lambda v: f"{int(v)} mm",
    colors="dimgray",
    use_clabeltext=True,
    zorder=8,
)
# 强制设置标签字体和大小（clabel 的 fontsize 参数不总生效，需逐个设置）
_font_path = os.path.join(BASE_DIR, "LXGWWenKai-Regular.ttf")
label_font = fm.FontProperties(fname=_font_path, size=4) if os.path.exists(_font_path) else fm.FontProperties(size=4)
for lbl in labels:
    lbl.set_fontproperties(label_font)
    lbl.set_fontsize(4)
    lbl.set_bbox(dict(facecolor="white", edgecolor="none",
                      alpha=0.75, pad=0.3, boxstyle="round,pad=0.08"))

# 省份名称
for _, row in provmap.iterrows():
    try:
        pt = row.geometry.representative_point()
        prov_name = row.get("省", "")
        if prov_name:
            ax.text(pt.x, pt.y, prov_name,
                    fontsize=5.5, color="gray",
                    ha="center", va="center",
                    fontproperties=CNFONT, zorder=10)
    except Exception:
        pass

# 此处代码需下载讲义材料查看～

# 等比例保存
xlim = ax.get_xlim(); ylim = ax.get_ylim()
dw = xlim[1] - xlim[0]; dh = ylim[1] - ylim[0]
fig.set_size_inches(10, 10 * dh / dw)
ax.set_xlim(xlim); ax.set_ylim(ylim)

fig.savefig("pic3.png", dpi=400, bbox_inches="tight",
            facecolor="white", edgecolor="none")
print("pic3.png 已保存")
plt.show()

九、直接运行 Python 脚本

除了在 Rmd 文档中逐步演示，也可以直接运行 main.py 脚本一次性生成所有地图：

# 在虚拟环境中运行 main.py
py_run_file("main.py")

或者在终端中：

# 激活虚拟环境后运行
.venv/bin/python main.py

---

附录：R 与 Python 函数对照表

R 代码	Python 等价	说明
rast("pre_2021.nc")	xr.open_dataset(nc_path)	读取 NetCDF 栅格
app(rst, mean, na.rm=T)	.mean(dim=..., skipna=True)	计算多波段均值
st_union(prov)	unary_union(prov.geometry)	合并多边形
nngeo::st_remove_holes()	Polygon(geom.exterior)	去除多边形孔洞
terra::crop(vect(cn))	rasterio.mask.mask(crop=True)	裁剪栅格
terra::mask(vect(cn))	rasterio.mask.mask(all_touched=True)	掩膜裁剪
terra::extract(fun=mean)	geometry_mask + np.nanmean	区域均值提取
rasterToContour	ax.contour（直接在 Albers axes 上）	生成等值线
aggregate(fact=10)	scipy.ndimage.zoom	降采样
rasterToPolygons	shapely.geometry.box 逐像元	栅格转多边形
st_transform(mycrs)	.to_crs(MYCRS)	坐标系转换
st_intersection	gpd.overlay(how="intersection")	空间相交
geometry * 0.5 + c(tx, ty)	affine_scale + translate	几何变换
geom_textsf	ax.clabel(inline=True, fmt=...)	等值线文字标签
annotation_scale	手动绘制矩形比例尺	比例尺
annotation_north_arrow	ax.annotate + ax.text	指北针
scale_fill_gradientn	LinearSegmentedColormap.from_list	渐变色板
ggsave(dpi=400)	fig.savefig(dpi=400)	保存高分辨率图

# 为生成的图片添加水印（如需要）
# lapply(fs::dir_ls(regexp = "png"), rstatatools::addrswm)

点击这里跳转到 RStata 短书平台获取附件：使用 Python 绘制填充地图、栅格地图 + 等降水量线