『matplotlib』利用basemap嵌入地图信息

import importlib
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from matplotlib.colors import LinearSegmentedColormap
from osgeo import gdal

__all__ = ["fig"]

# plt.rcParams['font.family'] = 'SimHei'
# plt.rcParams['font.style'] = 'normal'  # 或者 "italic"
plt.rcParams["font.size"] = 15
plt.rc('font', size=15, weight='bold')
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题,或者转换负号为字符串


def register_my_cmap():
    color_list = ["#FFFFFF", "#87CEFA", "#00FA9A", "#FFFF00", "#FF4500", "#FF0000"]
    my_cmap = LinearSegmentedColormap.from_list('rain', color_list)
    cm.register_cmap(cmap=my_cmap)


def fig(src: str, file: str):
    r"""
    :param src: 要读取的TIFF文件地址
    :param file: 输出的图片的存储地址
    :return: None
    """
    if isinstance(src, Path):
        src = str(src)

    img_info = gdal.Open(src)
    im_width = img_info.RasterXSize  # 获取宽度，数组第二维，左右方向元素长度，代表经度范围
    im_height = img_info.RasterYSize  # 获取高度，数组第一维，上下方向元素长度，代表纬度范围

    im_GeoTransform = img_info.GetGeoTransform()  # 获取仿射矩阵，含有 6 个元素的元组

    img = img_info.GetRasterBand(1).ReadAsArray(0, 0, im_width, im_height)  # 得到数据的值信息

    # 获取该数据的经纬度信息
    geo_info = []
    for i in range(im_height):
        row = []
        for j in range(im_width):
            xGeo = im_GeoTransform[0] + j * im_GeoTransform[1] + i * im_GeoTransform[2]
            yGeo = im_GeoTransform[3] + j * im_GeoTransform[4] + i * im_GeoTransform[5]
            col = [xGeo, yGeo]
            row.append(col)
        geo_info.append(row)

    geo_info = np.array(geo_info)
    # 至此，geo_info 的维度信息为：(im_height, im_width, 2)

    longitudes_list = geo_info[0, :, 0]  # 获取经度信息列表
    latitudes_list = geo_info[:, 0, 1]  # 获取维度信息列表

    fig = plt.figure(figsize=(16, 9))  # 创建画布

    axes = fig.add_subplot(111)  # 添加子图，这里 111 表示就一张图，返回该子图的坐标系对象

    lon, lat = np.meshgrid(longitudes_list, latitudes_list)  # 生成和 img 相同 shape 的矩阵，分别表示每个位置的经纬度

    basemap = importlib.import_module("mpl_toolkits.basemap")  # 动态导入basemap，因为它和 gdal 冲突

    m = basemap.Basemap(lat_0=0, lon_0=0, ax=axes)  # 前两个参数是定位中心，axes 是你绘图所在的 坐标轴系

    xi, yi = m(lon, lat)  # 将地理坐标替换磨人的数组下标索引

    # 将特定的值设置为mask
    mask = img <= 0
    img_mask = np.ma.array(img, mask=mask)  # 生成带掩码的矩阵

    register_my_cmap()
    # 这里强调一下 levels 参数，可以控制颜色条的分级，看最终示例图就明白了
    gci = m.contourf(x=xi, y=yi, data=np.squeeze(img_mask), levels=np.arange(0, 11, 1), cmap=plt.cm.twilight)

    r"""
    画经纬线及其标签的参数含义：
        color: 经纬线颜色，默认 black
        textcolor： 经纬度标注颜色， 默认 black
        linewidth： 经纬线宽，默认 1
        labels：控制经纬度标注显示位置，[left, right, top, bottom]，显示是 1，不显示是 0
        labelstyle：设置为 "+/-" 表示标注带正负号，否则带 N/S
        fmt： 标注显示的格式字符串
        fontsize：字体大小
        xoffset: 设置左右标注偏移图片主体的距离
        yoffset: 设置上下标注偏移图片主体的距离
        ax： 传入的Axes示实例，表明你是想在哪个坐标系上画经纬网
    """

    m.drawparallels(circles=np.arange(-90., 91., 30.), labels=[1, 0, 0, 0], fontsize=12, ax=axes)  # 画纬线
    m.drawmeridians(meridians=np.arange(-180., 185., 30.), labels=[0, 0, 1, 0], fontsize=12, ax=axes)  # 画经线

    # 添加海岸线,省州边界以及国家行政边界
    m.drawcoastlines()
    # m.drawstates()
    # m.drawcountries()

    # 添加 colorbar
    cbar = m.colorbar(gci, location='bottom', pad="5%", fig=fig, ax=axes)
    cbar.set_ticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
    colorbar_label_font_dict = {
        "family": "times new roman",
        "color": "darkred",
        "style": 'italic',
        "weight": "bold",
        "size": "16",
    }
    cbar.set_label(r'$rain(mm\cdot h^{-1})$', fontdict=colorbar_label_font_dict)

    plt.suptitle(f"Global Precipitation on {file[:4]}-{file[4:6]}-{file[6:]}", family="times new roman",
                 size=24, weight="bold", style="italic", color="orange",
                 verticalalignment="center",  # 水平对齐方式
                 horizontalalignment="center",  # 垂直对齐方式
                 y=0.96)
    # plt.show()
    plt.savefig(file)
    plt.close()