Python | 为什么Cartopy绘制的地图不能用?

起因

之前开组会，遇到一些绘制中国地图的图，然后老师提醒，cartopy或者NCL绘制的中国的地图不能用。

为什么不能用，却没有明说。

知其然，不知其所以然，于是便去查了一下。

Cartopy和NCL中的问题

可以看出，Cartopy和NCL的地图这里缺一点，那里缺一点。

以后如果有人问为什么不能用Cartopy或者NCL中的数据来绘制地图，直接发这个图给他。

正确的绘图

在这里，可以使用cnmap来绘制，cnmaps 是一个可以让中国地图画起来更丝滑的地图类 python 扩展包

cnmaps 最简单也最快的安装方法是使用 pip 来安装 cnmaps：

pip install -U cnmaps

绘制边界

用最简单直接的方式，来绘制你的第一张中国地图。

import cartopy.crs as ccrsimport matplotlib.pyplot as pltfrom cnmaps import get_adm_maps, draw_mapsfig = plt.figure(figsize=(10,10))ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_maps(get_adm_maps(level='国')) plt.show()

审图号

当然，使用上面的代码绘制也会有个问题，带中国地图的稿件可能会遇到审图号的问题，这个可以去网上去下载。

这里我们使用气象家园一位用户提供的数据：

有需要可以去该网站下载。

可以看出审图号是：GS（2024）0650号

这里使用Python代码来绘制：

from pathlib import Pathimport matplotlib.pyplot as pltimport numpy as npimport shapefileDEFAULT_DATA_DIR = Path(r"中国标准行政区划数据GS（2024）0650号_派生数据")DEFAULT_OUTPUT = Path(r"/output/china_map.png")DEFAULT_LANGUAGE = "zh"DEFAULT_PROVINCE_BOUNDARIES = TrueDEFAULT_DPI = 600MAIN_EXTENT = (70.0, 149.0, 13.0, 55.0)INSET_EXTENT = (106.0, 122.0, 1.0, 23.5)FONT_CANDIDATES = ["Microsoft YaHei","SimHei","Arial Unicode MS","DejaVu Sans",]defmiller_project(lon, lat):    lon_arr = np.asarray(lon, dtype=float)    lat_arr = np.clip(np.asarray(lat, dtype=float), -89.999, 89.999)    lon_rad = np.deg2rad(lon_arr)    lat_rad = np.deg2rad(lat_arr)    x = lon_rad    y = 1.25 * np.log(np.tan(np.pi / 4.0 + 0.4 * lat_rad))return x, ydefproject_extent(extent):    lon_min, lon_max, lat_min, lat_max = extent    x0, _ = miller_project([lon_min], [lat_min])    x1, _ = miller_project([lon_max], [lat_min])    _, y0 = miller_project([lon_min], [lat_min])    _, y1 = miller_project([lon_min], [lat_max])return float(x0[0]), float(x1[0]), float(y0[0]), float(y1[0])defconfigure_matplotlib():    plt.rcParams["font.sans-serif"] = FONT_CANDIDATES    plt.rcParams["axes.unicode_minus"] = False    plt.rcParams["savefig.facecolor"] = "white"    plt.rcParams["figure.facecolor"] = "white"defresolve_data_dir(data_dir):if data_dir:        candidate = Path(data_dir).expanduser()if candidate.exists():return candidateraise FileNotFoundError(f"数据目录不存在: {candidate}")if DEFAULT_DATA_DIR.exists():return DEFAULT_DATA_DIRdefbuild_layer_paths(data_dir):    layers = {"border": data_dir / "国界线+海岸线" / "中国_陆地国界线(不含未定国界).shp","coast": data_dir / "国界线+海岸线" / "中国_海岸线new.shp","nine_dash": data_dir / "缓冲区(国界线九段线)" / "15-30km-九段线.shp","undetermined": data_dir / "未定国界" / "新疆_未定国界(虚线)-版本2.shp","province": data_dir / "省界线.shp",    }    missing = [path for path in layers.values() ifnot path.exists()]if missing:        missing_text = "\n".join(str(path) for path in missing)raise FileNotFoundError(f"以下 shp 文件缺失:\n{missing_text}")return layersdefshape_parts(shape):    points = np.asarray(shape.points, dtype=float)if points.size == 0:return    indices = list(shape.parts) + [len(points)]for start, end in zip(indices[:-1], indices[1:]):        segment = points[start:end]if len(segment) >= 2:yield segmentdefplot_line_shapefile(    ax,    shp_path,    *,    color="black",    linewidth=1.0,    linestyle="-",    zorder=3,):    reader = shapefile.Reader(str(shp_path))for shape in reader.iterShapes():for segment in shape_parts(shape):            x, y = miller_project(segment[:, 0], segment[:, 1])            ax.plot(                x,                y,                color=color,                linewidth=linewidth,                linestyle=linestyle,                solid_capstyle="round",                zorder=zorder,            )defplot_point_shapefile(    ax,    shp_path,    *,    color="black",    size=0.35,    zorder=4,):    reader = shapefile.Reader(str(shp_path))    xs = []    ys = []for shape in reader.iterShapes():        points = np.asarray(shape.points, dtype=float)if len(points) == 0:continue        x, y = miller_project(points[:, 0], points[:, 1])        xs.append(x)        ys.append(y)ifnot xs:return    ax.scatter(        np.concatenate(xs),        np.concatenate(ys),        s=size,        c=color,        linewidths=0,        zorder=zorder,    )defformat_lon_label(value):returnf"{int(value)}°E"defformat_lat_label(value):returnf"{int(value)}°N"defsetup_axis(    ax,    extent,    *,    show_tick_labels,    tick_fontsize,):    lon_min, lon_max, lat_min, lat_max = extent    x0, x1, y0, y1 = project_extent(extent)    ax.set_xlim(x0, x1)    ax.set_ylim(y0, y1)    ax.set_aspect("equal", adjustable="box")    xticks = np.array([80, 100, 120, 140], dtype=float)    yticks = np.array([10, 20, 30, 40, 50], dtype=float)    xticks = xticks[(xticks >= lon_min) & (xticks <= lon_max)]    yticks = yticks[(yticks >= lat_min - 5) & (yticks <= lat_max)]    xtick_pos, _ = miller_project(xticks, np.full_like(xticks, lat_min))    _, ytick_pos = miller_project(np.full_like(yticks, lon_min), yticks)    ax.set_xticks(xtick_pos)    ax.set_yticks(ytick_pos)if show_tick_labels:        ax.set_xticklabels([format_lon_label(val) for val in xticks], fontsize=tick_fontsize)        ax.set_yticklabels([format_lat_label(val) for val in yticks], fontsize=tick_fontsize)else:        ax.set_xticklabels([])        ax.set_yticklabels([])    ax.grid(True, color="#808080", linewidth=0.45, alpha=0.55)for spine in ax.spines.values():        spine.set_linewidth(0.8)        spine.set_color("black")defdraw_layers(    ax,    layers,    *,    province_boundaries,    border_linewidth,    nine_dash_linewidth,    coast_size,    undetermined_linewidth,):if province_boundaries:        plot_line_shapefile(            ax,            layers["province"],            color="#5A5A5A",            linewidth=0.35,            linestyle="-",            zorder=2,        )    plot_line_shapefile(        ax,        layers["border"],        color="black",        linewidth=border_linewidth,        linestyle="-",        zorder=4,    )    plot_line_shapefile(        ax,        layers["nine_dash"],        color="black",        linewidth=nine_dash_linewidth,        linestyle="-",        zorder=4,    )# 海岸线用散点而不是连续折线，能避免群岛区域在导出时出现跨岛连线和十字伪影。    plot_point_shapefile(        ax,        layers["coast"],        color="black",        size=coast_size,        zorder=5,    )    plot_line_shapefile(        ax,        layers["undetermined"],        color="black",        linewidth=undetermined_linewidth,        linestyle=":",        zorder=4,    )defadd_inset(    fig,    layers,    *,    language,    province_boundaries,):    inset = fig.add_axes([0.70, 0.11, 0.21, 0.18])    setup_axis(inset, INSET_EXTENT, show_tick_labels=False, tick_fontsize=8)    draw_layers(        inset,        layers,        province_boundaries=province_boundaries,        border_linewidth=0.55,        nine_dash_linewidth=0.55,        coast_size=0.22,        undetermined_linewidth=0.45,    )    label = "南海诸岛"if language == "zh"else"Nanhai\nZhudao"    label_x = 116.0if language == "zh"else115.5    label_y = 2.8    tx, ty = miller_project([label_x], [label_y])    inset.text(        tx[0],        ty[0],        label,        fontsize=7if language == "zh"else8,        fontweight="bold",        ha="center",        va="bottom",    )return insetdefdraw_china_map(    data_dir,    output,    *,    language,    province_boundaries,    dpi,):    configure_matplotlib()    layers = build_layer_paths(data_dir)    fig = plt.figure(figsize=(23.159861111111113 / 2.54, 18.39736111111111 / 2.54))    ax = fig.add_axes([0.08, 0.08, 0.84, 0.84])    setup_axis(ax, MAIN_EXTENT, show_tick_labels=True, tick_fontsize=10)    draw_layers(        ax,        layers,        province_boundaries=province_boundaries,        border_linewidth=1.0,        nine_dash_linewidth=1.0,        coast_size=0.35,        undetermined_linewidth=0.55,    )    add_inset(fig, layers, language=language, province_boundaries=province_boundaries)    output.parent.mkdir(parents=True, exist_ok=True)    fig.savefig(output, dpi=dpi, bbox_inches="tight", pad_inches=0.06)    plt.close(fig)return outputdefmain():    data_dir = resolve_data_dir(None)    output = DEFAULT_OUTPUT.expanduser()    result = draw_china_map(        data_dir,        output,        language=DEFAULT_LANGUAGE,        province_boundaries=DEFAULT_PROVINCE_BOUNDARIES,        dpi=DEFAULT_DPI,    )    print(f"图片已保存: {result}")if __name__ == "__main__":    main()

可以根据上述代码修改到自己想要的位置。

References

• https://github.com/cnmetlab/cnmaps
• http://bbs.06climate.com/forum.php?mod=viewthread&tid=111711
• https://github.com/huangynj/NCL-Chinamap
• https://www.mdpi.com/2073-4433/17/3/327

往期回顾

• Python | 数据压缩 | 八面体高斯网格规则化
• Python | 压缩近3倍netcdf文件大小
• Python | Skyborn | 涌现约束法
• Python | 球面谐波 | 流函数 | 势函数
• Python | Skyborn | 对流层顶高度的计算
• Python | Skyborn | 地转风 | 地转平
• Python | Skyborn | 缺测值填充
• Python | Skyborn | 空间趋势 | 快速计算
• Python | Skyborn | 涌现约束法
• Python | 计算台风与海岸线距离