论文发表必备技能：空间数据（tif图像、shp边界）多子图可视化-基于cartopy与matplotlib的python包_matplotlib 绘制 tiff

零、背景

众所周知，我正在写毕业论文，刚画完一张2x2子图布局的偏相关系数空间分布图。

我的变量简介：具体地，我有三个因变量分别是SOS, EOS, LOS，四个自变量分别是tem, pre, IET, IEP，因此，我需要画三张大图，每张图内包含四张子图，即三个因变量与四个自变量分别的偏相关系数分布图。

结构：我的代码总共分为三部分，包括①导入必备包以及函数准备、②数据准备、以及③画图。其内容涉及到，可视化tif图像、shp边界、绘制colorbar、在图框上显示经纬度等优雅操作。

现新鲜出炉，分享给大家，热忱期盼交流指正~

一、导入必备包以及函数准备

from osgeo import gdal
import os
from matplotlib.ticker import FuncFormatter
import matplotlib as mpl
import numpy as np
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from matplotlib.ticker import MultipleLocator
from cartopy.io.shapereader import Reader
from cartopy import crs as ccrs
from matplotlib.ticker import MaxNLocator
from matplotlib.colors import Normalize
import matplotlib.patches as patches
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
 
'''
目标：获取所有.tif图像路径
输入：tif所在文件夹绝对路径
输出：包含该文件夹中所有tif文件绝对路径的list
'''
def get_paths(input_path):
    file_list = [os.path.join(input_path, file) for file in os.listdir(input_path) if file.endswith(".tif")]
    return file_list
 
'''
读写遥感数据
'''
class Grid:
    # 读取图像文件
    def read_img(self, filename):
        data = gdal.Open(filename)  # 打开文件
        im_width = data.RasterXSize  # 读取图像行数
        im_height = data.RasterYSize  # 读取图像列数
 
        im_geotrans = data.GetGeoTransform()
        # 仿射矩阵，左上角像素的大地坐标和像素分辨率。
        # 共有六个参数，分表代表左上角x坐标；东西方向上图像的分辨率；如果北边朝上，地图的旋转角度，0表示图像的行与x轴平行；左上角y坐标；
        # 如果北边朝上，地图的旋转角度，0表示图像的列与y轴平行；南北方向上地图的分辨率。
        im_proj = data.GetProjection()  # 地图投影信息
        im_data = data.ReadAsArray(0, 0, im_width, im_height)  # 此处读取整张图像
        # ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)
        # 读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。
        del data  # 释放内存，如果不释放，在arcgis，envi中打开该图像时会显示文件被占用
 
        return im_proj, im_geotrans, im_data
 
    def write_img(self, filename, im_proj, im_geotrans, im_data):
        # filename-创建的新影像
        # im_geotrans,im_proj该影像的参数，im_data，被写的影像
        # 写文件，以写成tiff为例
        # gdal数据类型包括
        # gdal.GDT_Byte,
        # gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        # gdal.GDT_Float32, gdal.GDT_Float64
 
        # 判断栅格数据的类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
        if len(im_data.shape) == 3:  # len(im_data.shape)表示矩阵的维数
            im_bands, im_height, im_width = im_data.shape  # （维数，行数，列数）
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape  # 一维矩阵
 
        #         创建文件
        driver = gdal.GetDriverByName('GTiff')  # 数据类型必须有，因为要计算需要多大内存空间
        data = driver.Create(filename, im_width, im_height, im_bands, datatype)
        data.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        data.SetProjection(im_proj)  # 写入投影
        if im_bands == 1:
            data.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                data.GetRasterBand(i + 1).WriteArray(im_data[i])
        del data
a = Grid()
 
'''
目标：寻找映射最大值、最小值
输入：四个数组（子图中绘制的四个tif）
输出：四个数组全体的最大值、最小值
'''
def find_global_max_min(arr1, arr2, arr3, arr4):
    # 初始化最大值和最小值
    global_max = None
    global_min = None
 
    # 将四个数组放入列表中
    arrays = [arr1, arr2, arr3, arr4]
 
    # 遍历每个数组并查找最大值和最小值
    for arr in arrays:
        local_max = np.nanmax(arr)
        local_min = np.nanmin(arr)
 
        # 更新全局最大值和最小值
        if global_max is None or local_max > global_max:
            global_max = local_max
        if global_min is None or local_min < global_min:
            global_min = local_min
 
    return global_max, global_min
 
'''
输入：sos_data_lists, eos_data_lists, los_data_lists
输出：数组正值与负值占比
'''
def calculate_percentage(data_list):
    result = []
 
    for data in data_list:
        positive_count = 0
        negative_count = 0
        total_count = 0
 
        for row in data:
            for value in row:
                if -1 <= value <= 1:
                    total_count += 1
                    if value > 0:
                        positive_count += 1
                    elif value < 0:
                        negative_count += 1
 
        if total_count == 0:
            result.append("No valid data in the array")
        else:
            positive_percentage = (positive_count / total_count) * 100
            negative_percentage = (negative_count / total_count) * 100
            result.append(positive_percentage)
            result.append(negative_percentage)
 
    return result

二、数据准备

#tif数据所在文件夹
sos_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\sos')
eos_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\eos')
los_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\los')
 
#获取数组
#sos
proj, geotrans,sos_tem_corr = a.read_img(sos_paths[-1])
_,_,sos_pre_corr = a.read_img(sos_paths[-2])
_,_,sos_iet_corr = a.read_img(sos_paths[-3])
_,_,sos_iep_corr = a.read_img(sos_paths[-4])
#eos
_,_,eos_tem_corr = a.read_img(eos_paths[-1])
_,_,eos_pre_corr = a.read_img(eos_paths[-2])
_,_,eos_iet_corr = a.read_img(eos_paths[-3])
_,_,eos_iep_corr = a.read_img(eos_paths[-4])
#los
_,_,los_tem_corr = a.read_img(los_paths[-1])
_,_,los_pre_corr = a.read_img(los_paths[-2])
_,_,los_iet_corr = a.read_img(los_paths[-3])
_,_,los_iep_corr = a.read_img(los_paths[-4])
 
#按制图顺序，整理数据顺序
#sos
sos_data_lists = [sos_tem_corr, sos_pre_corr,sos_iet_corr,sos_iep_corr]
sos_data_lists = np.array(sos_data_lists)
#eos
eos_data_lists = [eos_tem_corr, eos_pre_corr,eos_iet_corr,eos_iep_corr]
eos_data_lists = np.array(eos_data_lists)
#los
los_data_lists = [los_tem_corr, los_pre_corr,los_iet_corr,los_iep_corr]
los_data_lists = np.array(los_data_lists)
 
#寻找sos/eos/los偏相关系数最大最小值
sos_max, sos_min = find_global_max_min(sos_tem_corr, sos_pre_corr, sos_iet_corr,sos_iep_corr)
eos_max, eos_min = find_global_max_min(eos_tem_corr, eos_pre_corr, eos_iet_corr,eos_iep_corr)
los_max, los_min = find_global_max_min(los_tem_corr, los_pre_corr, los_iet_corr,los_iep_corr)
print("SOS 数据集全局最大值:", sos_max, "\nSOS 数据集全局最小值:", sos_min, "\nEOS 数据集全局最大值:", eos_max, "\nEOS 数据集全局最小值:", eos_min, "\nLOS 数据集全局最大值:", los_max, "\nLOS 数据集全局最小值:", los_min)
 
#边界数据
shp_path = r'D:\Hydro_essay\bound\bound_without_proj.shp'
 
#获取各子图正值、负值占比
#sos-tem/pre/iet/iep
sos_perc = calculate_percentage(sos_data_lists)
sos_p_perc = []#4个值
sos_p_perc.extend([sos_perc[0], sos_perc[2], sos_perc[4], sos_perc[6]])
sos_n_perc = []#4个值
sos_n_perc.extend([sos_perc[1], sos_perc[3], sos_perc[5], sos_perc[7]])
#eos-tem/pre/iet/iep
eos_perc = calculate_percentage(eos_data_lists)
eos_p_perc = []#4个值
eos_p_perc.extend([eos_perc[0], eos_perc[2], eos_perc[4], eos_perc[6]])
eos_n_perc = []#4个值
eos_n_perc.extend([eos_perc[1], eos_perc[3], eos_perc[5], eos_perc[7]])
#los-tem/pre/iet/iep
los_perc = calculate_percentage(los_data_lists)
los_p_perc = []#4个值
los_p_perc.extend([los_perc[0], los_perc[2], los_perc[4], los_perc[6]])
los_n_perc = []#4个值
los_n_perc.extend([los_perc[1], los_perc[3], los_perc[5], los_perc[7]])

三、开始画图

这部分我将colorbar的title写为中文的（可能是因为要毕业了吧）。如果你需要英文出图，只需要替换一下fontdict即可。

# 字体和样式设置
size1 = 4
fontdict1 = {'weight': 'bold', 'size': 4, 'color': 'k', 'family': 'Times New Roman'}
fontdict2 = {'size': 12, 'color': 'k', 'family': 'Times New Roman'}
fontdict3 = {'size': 6, 'color': 'k', 'family': 'Times New Roman'}
fontdict4 = {'weight': 'bold', 'size': 3.9, 'color': 'k', 'family': 'Times New Roman'}
fontdict5 = {'weight': '1000', 'size': 4.5, 'color': 'k', 'family': "SimSun"}#colorbar
fontdict6 = {'weight': '1000', 'size': 3.9, 'color': 'k', 'family': "SimSun"}#text
#全体字体
mpl.rcParams.update({
    'text.usetex': False,
    'font.family': 'serif',
    'font.serif': ['Times New Roman'],
    'font.size': size1,
    'mathtext.fontset': 'stix',
})
# 投影设置
proj = ccrs.PlateCarree()
# 图框经纬度范围
extent = [102.7, 104.7, 27.2, 29.2]#注意四个数字的排序！
#tif图像经纬度设置
extent_tif = [102.84521586, 104.43331875, 27.37916667,  28.90833333]
# 图布大小设置
fig = plt.figure(figsize=(25, 16), dpi=300)
# 子图布局参数
left_margin_adjusted = 0.345 #左列子图左侧图框距图布边儿距离
bottom_margin_adjusted = 0.098 #下列子图下侧图框距图布边儿距离
colorbar_left_adjusted = 0.82  # 色标的左边界位置调整
# 建立轴
ax_a = fig.add_axes([0.05, 0.52, 0.39, 0.39], projection=proj)
ax_b = fig.add_axes([left_margin_adjusted, 0.52, 0.39, 0.39], projection=proj)
ax_c = fig.add_axes([0.05, bottom_margin_adjusted, 0.39, 0.39], projection=proj)
ax_d = fig.add_axes([left_margin_adjusted, bottom_margin_adjusted, 0.39, 0.39], projection=proj)
# 设置轴上tick的长度和宽度，以及经纬度数字与轴的距离
tick_length = 1.2  # 刻度线的长度
tick_width = 0.4  # 刻度线的宽度
label_pad = 1  # 刻度标签与轴的距离
axis_line_width = 0.35  # 轴线（ticks）宽度
# 设置归一化对象，该对象用于统一不同图像的色标-修改xos_min;xos_max替换数据1/4！！！
norm = Normalize(vmin=sos_min, vmax=sos_max)
#坐标设置
#我希望左上子图只有左侧和上侧显示坐标；右上子图只有右侧和上侧显示坐标等...四个子图各不相同
# 左上子图设置
ax_a.tick_params(axis='x', which='both', bottom=False, top=True, labelbottom=False, labeltop=True,length=tick_length, width=tick_width, pad=label_pad)
ax_a.tick_params(axis='y', which='both', left=True, right=False, labelleft=True, labelright=False,length=tick_length, width=tick_width, pad=label_pad)
ax_a.xaxis.set_ticks_position('top')
ax_a.yaxis.set_ticks_position('left')
# 右上子图设置
ax_b.tick_params(axis='x', which='both', bottom=False, top=True, labelbottom=False, labeltop=True,length=tick_length, width=tick_width, pad=label_pad)
ax_b.tick_params(axis='y', which='both', left=False, right=True, labelleft=False, labelright=True,length=tick_length, width=tick_width, pad=label_pad)
ax_b.xaxis.set_ticks_position('top')
ax_b.yaxis.set_ticks_position('right')
# 左下子图设置
ax_c.tick_params(axis='x', which='both', bottom=True, top=False, labelbottom=True, labeltop=False,length=tick_length, width=tick_width, pad=label_pad)
ax_c.tick_params(axis='y', which='both', left=True, right=False, labelleft=True, labelright=False,length=tick_length, width=tick_width, pad=label_pad)
ax_c.xaxis.set_ticks_position('bottom')
ax_c.yaxis.set_ticks_position('left')
# 右下子图设置
ax_d.tick_params(axis='x', which='both', bottom=True, top=False, labelbottom=True, labeltop=False,length=tick_length, width=tick_width, pad=label_pad)
ax_d.tick_params(axis='y', which='both', left=False, right=True, labelleft=False, labelright=True,length=tick_length, width=tick_width, pad=label_pad)
ax_d.xaxis.set_ticks_position('bottom')
ax_d.yaxis.set_ticks_position('right')
#子图title（我要画三张图，所以三个list）
sos_titles = ['(a) SOS_tem','(b) SOS_pre','(c) SOS_IET','(d) SOS_IEP']
eos_titles = ['(a) EOS_tem','(b) EOS_pre','(c) EOS_IET','(d) EOS_IEP']
los_titles = ['(a) LOS_tem','(b) LOS_pre','(c) LOS_IET','(d) LOS_IEP']
 
# 四个子图设置坐标线、经纬度、以及绘制shp文件和tif文件
for idx,ax in enumerate([ax_a, ax_b, ax_c, ax_d]):
    print(idx)
    ax.set_extent(extent, crs=proj)  # 设置地图显示的范围
    ax.set_xticks(np.arange(int(np.floor(extent[0])), int(np.ceil(extent[1])) + 1, 1), crs=proj)
    ax.set_yticks(np.arange(int(np.floor(extent[2])), int(np.ceil(extent[3])) + 1, 1), crs=proj)
    ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    #ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,linewidth=.8, linestyle=':', color='k', alpha=0.5)
    # 设置刻度只显示整数
    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
    # 设置轴线宽度
    for spine in ax.spines.values():
        spine.set_linewidth(axis_line_width)
    #小刻度线
    xminorticks = MultipleLocator(1)
    yminorticks = MultipleLocator(1)
    ax.xaxis.set_minor_locator(xminorticks)
    ax.yaxis.set_minor_locator(yminorticks)
    # 绘制矢量边界
    ax.add_geometries(Reader(shp_path).geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.25)
    #绘制tif-修改xos_data_lists替换数据2/4！！！
    ax.imshow(sos_data_lists[idx], extent=extent_tif, cmap='coolwarm', norm=norm, transform=proj, interpolation='nearest')
    #注：在使用imshow或contourf函数时，指定的extent必须准确匹配tif数据的实际地理范围,因此指定新的extent
    # 设置标题，并调整位置-修改xos_titles替换数据3/4!!!
    ax.text(0.5, 0.92, sos_titles[idx], transform=ax.transAxes, fontdict=fontdict1, ha='center')
    #添加正值负值比例-修改xos_p_perc；xos_n_perc[idx]替换数据4/4！！！
    # Red rectangle for > 0
    red_patch = patches.Rectangle((0.5, 0.18), 0.075, 0.035, transform=ax.transAxes, color='red', zorder=5)
    ax.add_patch(red_patch)
    red_ratio = sos_p_perc[idx]
    blue_ratio = sos_n_perc[idx]
    ax.text(0.6, 0.19, f'{red_ratio:.2f}% ', transform=ax.transAxes, fontdict=fontdict4, va='center')
    ax.text(0.8, 0.19, f'(正值)', transform=ax.transAxes, fontdict=fontdict6, va='center')
    # Blue rectangle for < 0
    blue_patch = patches.Rectangle((0.5, 0.10), 0.075, 0.035, transform=ax.transAxes, color='blue', zorder=5)
    ax.add_patch(blue_patch)
    ax.text(0.6, 0.11, f'{blue_ratio:.2f}% ', transform=ax.transAxes, fontdict=fontdict4, va='center')
    ax.text(0.8, 0.11, f'(负值)', transform=ax.transAxes, fontdict=fontdict6, va='center')
 
print('drawing colorbar')
# 设置色标位置和大小
colorbar_width = 0.015
colorbar_height = 0.8
ax_colorbar = fig.add_axes([colorbar_left_adjusted-0.08, bottom_margin_adjusted, colorbar_width, colorbar_height])
# 创建色标
sm = mpl.cm.ScalarMappable(cmap='coolwarm', norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, cax=ax_colorbar)
# 设置色标刻度
cbar.set_ticks([-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8])
# 设置色标刻度标签格式为一位小数
cbar.ax.yaxis.set_major_formatter(FuncFormatter(lambda x, _: f'{x:.1f}'))
# 设置色标的刻度参数和移除边框
cbar.ax.tick_params(axis='y', which='both', length=tick_length, width=tick_width, pad=label_pad)
cbar.outline.set_linewidth(0)  # 移除色标的边框
# 设置色标标签
cbar.set_label('偏相关系数', fontdict=fontdict5)
 
plt.show()

四、结果展示

再来1个英文版本的吧（证明我都做啦哈哈）

五、全部代码集合

#%%
from osgeo import gdal
import os
from matplotlib.ticker import FuncFormatter
import matplotlib as mpl
import numpy as np
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from matplotlib.ticker import MultipleLocator
from cartopy.io.shapereader import Reader
from cartopy import crs as ccrs
from matplotlib.ticker import MaxNLocator
from matplotlib.colors import Normalize
import matplotlib.patches as patches
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
 
'''
目标：获取所有.tif图像路径
输入：tif所在文件夹绝对路径
输出：包含该文件夹中所有tif文件绝对路径的list
'''
def get_paths(input_path):
    file_list = [os.path.join(input_path, file) for file in os.listdir(input_path) if file.endswith(".tif")]
    return file_list
 
'''
读写遥感数据
'''
class Grid:
    # 读取图像文件
    def read_img(self, filename):
        data = gdal.Open(filename)  # 打开文件
        im_width = data.RasterXSize  # 读取图像行数
        im_height = data.RasterYSize  # 读取图像列数
 
        im_geotrans = data.GetGeoTransform()
        # 仿射矩阵，左上角像素的大地坐标和像素分辨率。
        # 共有六个参数，分表代表左上角x坐标；东西方向上图像的分辨率；如果北边朝上，地图的旋转角度，0表示图像的行与x轴平行；左上角y坐标；
        # 如果北边朝上，地图的旋转角度，0表示图像的列与y轴平行；南北方向上地图的分辨率。
        im_proj = data.GetProjection()  # 地图投影信息
        im_data = data.ReadAsArray(0, 0, im_width, im_height)  # 此处读取整张图像
        # ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)
        # 读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。
        del data  # 释放内存，如果不释放，在arcgis，envi中打开该图像时会显示文件被占用
 
        return im_proj, im_geotrans, im_data
 
    def write_img(self, filename, im_proj, im_geotrans, im_data):
        # filename-创建的新影像
        # im_geotrans,im_proj该影像的参数，im_data，被写的影像
        # 写文件，以写成tiff为例
        # gdal数据类型包括
        # gdal.GDT_Byte,
        # gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        # gdal.GDT_Float32, gdal.GDT_Float64
 
        # 判断栅格数据的类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
        if len(im_data.shape) == 3:  # len(im_data.shape)表示矩阵的维数
            im_bands, im_height, im_width = im_data.shape  # （维数，行数，列数）
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape  # 一维矩阵
 
        #         创建文件
        driver = gdal.GetDriverByName('GTiff')  # 数据类型必须有，因为要计算需要多大内存空间
        data = driver.Create(filename, im_width, im_height, im_bands, datatype)
        data.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        data.SetProjection(im_proj)  # 写入投影
        if im_bands == 1:
            data.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                data.GetRasterBand(i + 1).WriteArray(im_data[i])
        del data
a = Grid()
 
'''
目标：寻找映射最大值、最小值
输入：四个数组（子图中绘制的四个tif）
输出：四个数组全体的最大值、最小值
'''
def find_global_max_min(arr1, arr2, arr3, arr4):
    # 初始化最大值和最小值
    global_max = None
    global_min = None
 
    # 将四个数组放入列表中
    arrays = [arr1, arr2, arr3, arr4]
 
    # 遍历每个数组并查找最大值和最小值
    for arr in arrays:
        local_max = np.nanmax(arr)
        local_min = np.nanmin(arr)
 
        # 更新全局最大值和最小值
        if global_max is None or local_max > global_max:
            global_max = local_max
        if global_min is None or local_min < global_min:
            global_min = local_min
 
    return global_max, global_min
 
'''
输入：sos_data_lists, eos_data_lists, los_data_lists
输出：数组正值与负值占比
'''
def calculate_percentage(data_list):
    result = []
 
    for data in data_list:
        positive_count = 0
        negative_count = 0
        total_count = 0
 
        for row in data:
            for value in row:
                if -1 <= value <= 1:
                    total_count += 1
                    if value > 0:
                        positive_count += 1
                    elif value < 0:
                        negative_count += 1
 
        if total_count == 0:
            result.append("No valid data in the array")
        else:
            positive_percentage = (positive_count / total_count) * 100
            negative_percentage = (negative_count / total_count) * 100
            result.append(positive_percentage)
            result.append(negative_percentage)
 
    return result
 
#step1:数据准备
#tif数据所在文件夹
sos_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\sos')
eos_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\eos')
los_paths = get_paths(r'D:\code\python\phenology\gee\partial_corr\result\corr\los')
 
#获取数组
#sos
proj, geotrans,sos_tem_corr = a.read_img(sos_paths[-1])
_,_,sos_pre_corr = a.read_img(sos_paths[-2])
_,_,sos_iet_corr = a.read_img(sos_paths[-3])
_,_,sos_iep_corr = a.read_img(sos_paths[-4])
#eos
_,_,eos_tem_corr = a.read_img(eos_paths[-1])
_,_,eos_pre_corr = a.read_img(eos_paths[-2])
_,_,eos_iet_corr = a.read_img(eos_paths[-3])
_,_,eos_iep_corr = a.read_img(eos_paths[-4])
#los
_,_,los_tem_corr = a.read_img(los_paths[-1])
_,_,los_pre_corr = a.read_img(los_paths[-2])
_,_,los_iet_corr = a.read_img(los_paths[-3])
_,_,los_iep_corr = a.read_img(los_paths[-4])
 
#按制图顺序，整理数据顺序
#sos
sos_data_lists = [sos_tem_corr, sos_pre_corr,sos_iet_corr,sos_iep_corr]
sos_data_lists = np.array(sos_data_lists)
#eos
eos_data_lists = [eos_tem_corr, eos_pre_corr,eos_iet_corr,eos_iep_corr]
eos_data_lists = np.array(eos_data_lists)
#los
los_data_lists = [los_tem_corr, los_pre_corr,los_iet_corr,los_iep_corr]
los_data_lists = np.array(los_data_lists)
 
#寻找sos/eos/los偏相关系数最大最小值
sos_max, sos_min = find_global_max_min(sos_tem_corr, sos_pre_corr, sos_iet_corr,sos_iep_corr)
eos_max, eos_min = find_global_max_min(eos_tem_corr, eos_pre_corr, eos_iet_corr,eos_iep_corr)
los_max, los_min = find_global_max_min(los_tem_corr, los_pre_corr, los_iet_corr,los_iep_corr)
print("SOS 数据集全局最大值:", sos_max, "\nSOS 数据集全局最小值:", sos_min, "\nEOS 数据集全局最大值:", eos_max, "\nEOS 数据集全局最小值:", eos_min, "\nLOS 数据集全局最大值:", los_max, "\nLOS 数据集全局最小值:", los_min)
 
#边界数据
shp_path = r'D:\Hydro_essay\bound\bound_without_proj.shp'
 
#获取各子图正值、负值占比
#sos-tem/pre/iet/iep
sos_perc = calculate_percentage(sos_data_lists)
sos_p_perc = []#4个值
sos_p_perc.extend([sos_perc[0], sos_perc[2], sos_perc[4], sos_perc[6]])
sos_n_perc = []#4个值
sos_n_perc.extend([sos_perc[1], sos_perc[3], sos_perc[5], sos_perc[7]])
#eos-tem/pre/iet/iep
eos_perc = calculate_percentage(eos_data_lists)
eos_p_perc = []#4个值
eos_p_perc.extend([eos_perc[0], eos_perc[2], eos_perc[4], eos_perc[6]])
eos_n_perc = []#4个值
eos_n_perc.extend([eos_perc[1], eos_perc[3], eos_perc[5], eos_perc[7]])
#los-tem/pre/iet/iep
los_perc = calculate_percentage(los_data_lists)
los_p_perc = []#4个值
los_p_perc.extend([los_perc[0], los_perc[2], los_perc[4], los_perc[6]])
los_n_perc = []#4个值
los_n_perc.extend([los_perc[1], los_perc[3], los_perc[5], los_perc[7]])
#%%
#step2:开始画图
# 字体和样式设置
size1 = 4
fontdict1 = {'weight': 'bold', 'size': 4, 'color': 'k', 'family': 'Times New Roman'}
fontdict2 = {'size': 12, 'color': 'k', 'family': 'Times New Roman'}
fontdict3 = {'size': 6, 'color': 'k', 'family': 'Times New Roman'}
fontdict4 = {'weight': 'bold', 'size': 3.9, 'color': 'k', 'family': 'Times New Roman'}
fontdict5 = {'weight': '1000', 'size': 4.5, 'color': 'k', 'family': "SimSun"}#colorbar
fontdict6 = {'weight': '1000', 'size': 3.9, 'color': 'k', 'family': "SimSun"}#text
#全体字体
mpl.rcParams.update({
    'text.usetex': False,
    'font.family': 'serif',
    'font.serif': ['Times New Roman'],
    'font.size': size1,
    'mathtext.fontset': 'stix',
})
# 投影设置
proj = ccrs.PlateCarree()
# 图框经纬度范围
extent = [102.7, 104.7, 27.2, 29.2]#注意四个数字的排序！
#tif图像经纬度设置
extent_tif = [102.84521586, 104.43331875, 27.37916667,  28.90833333]
# 图布大小设置
fig = plt.figure(figsize=(25, 16), dpi=300)
# 子图布局参数
left_margin_adjusted = 0.345 #左列子图左侧图框距图布边儿距离
bottom_margin_adjusted = 0.098 #下列子图下侧图框距图布边儿距离
colorbar_left_adjusted = 0.82  # 色标的左边界位置调整
# 建立轴
ax_a = fig.add_axes([0.05, 0.52, 0.39, 0.39], projection=proj)
ax_b = fig.add_axes([left_margin_adjusted, 0.52, 0.39, 0.39], projection=proj)
ax_c = fig.add_axes([0.05, bottom_margin_adjusted, 0.39, 0.39], projection=proj)
ax_d = fig.add_axes([left_margin_adjusted, bottom_margin_adjusted, 0.39, 0.39], projection=proj)
# 设置轴上tick的长度和宽度，以及经纬度数字与轴的距离
tick_length = 1.2  # 刻度线的长度
tick_width = 0.4  # 刻度线的宽度
label_pad = 1  # 刻度标签与轴的距离
axis_line_width = 0.35  # 轴线（ticks）宽度
# 设置归一化对象，该对象用于统一不同图像的色标-修改xos_min;xos_max替换数据1/4！！！
norm = Normalize(vmin=sos_min, vmax=sos_max)
#坐标设置
#我希望左上子图只有左侧和上侧显示坐标；右上子图只有右侧和上侧显示坐标等...四个子图各不相同
# 左上子图设置
ax_a.tick_params(axis='x', which='both', bottom=False, top=True, labelbottom=False, labeltop=True,length=tick_length, width=tick_width, pad=label_pad)
ax_a.tick_params(axis='y', which='both', left=True, right=False, labelleft=True, labelright=False,length=tick_length, width=tick_width, pad=label_pad)
ax_a.xaxis.set_ticks_position('top')
ax_a.yaxis.set_ticks_position('left')
# 右上子图设置
ax_b.tick_params(axis='x', which='both', bottom=False, top=True, labelbottom=False, labeltop=True,length=tick_length, width=tick_width, pad=label_pad)
ax_b.tick_params(axis='y', which='both', left=False, right=True, labelleft=False, labelright=True,length=tick_length, width=tick_width, pad=label_pad)
ax_b.xaxis.set_ticks_position('top')
ax_b.yaxis.set_ticks_position('right')
# 左下子图设置
ax_c.tick_params(axis='x', which='both', bottom=True, top=False, labelbottom=True, labeltop=False,length=tick_length, width=tick_width, pad=label_pad)
ax_c.tick_params(axis='y', which='both', left=True, right=False, labelleft=True, labelright=False,length=tick_length, width=tick_width, pad=label_pad)
ax_c.xaxis.set_ticks_position('bottom')
ax_c.yaxis.set_ticks_position('left')
# 右下子图设置
ax_d.tick_params(axis='x', which='both', bottom=True, top=False, labelbottom=True, labeltop=False,length=tick_length, width=tick_width, pad=label_pad)
ax_d.tick_params(axis='y', which='both', left=False, right=True, labelleft=False, labelright=True,length=tick_length, width=tick_width, pad=label_pad)
ax_d.xaxis.set_ticks_position('bottom')
ax_d.yaxis.set_ticks_position('right')
#子图title（我要画三张图，所以三个list）
sos_titles = ['(a) SOS_tem','(b) SOS_pre','(c) SOS_IET','(d) SOS_IEP']
eos_titles = ['(a) EOS_tem','(b) EOS_pre','(c) EOS_IET','(d) EOS_IEP']
los_titles = ['(a) LOS_tem','(b) LOS_pre','(c) LOS_IET','(d) LOS_IEP']
 
# 四个子图设置坐标线、经纬度、以及绘制shp文件和tif文件
for idx,ax in enumerate([ax_a, ax_b, ax_c, ax_d]):
    print(idx)
    ax.set_extent(extent, crs=proj)  # 设置地图显示的范围
    ax.set_xticks(np.arange(int(np.floor(extent[0])), int(np.ceil(extent[1])) + 1, 1), crs=proj)
    ax.set_yticks(np.arange(int(np.floor(extent[2])), int(np.ceil(extent[3])) + 1, 1), crs=proj)
    ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    #ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,linewidth=.8, linestyle=':', color='k', alpha=0.5)
    # 设置刻度只显示整数
    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
    ax.yaxis.set_major_locator(MaxNLocator(integer=True))
    # 设置轴线宽度
    for spine in ax.spines.values():
        spine.set_linewidth(axis_line_width)
    #小刻度线
    xminorticks = MultipleLocator(1)
    yminorticks = MultipleLocator(1)
    ax.xaxis.set_minor_locator(xminorticks)
    ax.yaxis.set_minor_locator(yminorticks)
    # 绘制矢量边界
    ax.add_geometries(Reader(shp_path).geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.25)
    #绘制tif-修改xos_data_lists替换数据2/4！！！
    ax.imshow(sos_data_lists[idx], extent=extent_tif, cmap='coolwarm', norm=norm, transform=proj, interpolation='nearest')
    #注：在使用imshow或contourf函数时，指定的extent必须准确匹配tif数据的实际地理范围,因此指定新的extent
    # 设置标题，并调整位置-修改xos_titles替换数据3/4!!!
    ax.text(0.5, 0.92, sos_titles[idx], transform=ax.transAxes, fontdict=fontdict1, ha='center')
    #添加正值负值比例-修改xos_p_perc；xos_n_perc[idx]替换数据4/4！！！
    # Red rectangle for > 0
    red_patch = patches.Rectangle((0.5, 0.18), 0.075, 0.035, transform=ax.transAxes, color='red', zorder=5)
    ax.add_patch(red_patch)
    red_ratio = sos_p_perc[idx]
    blue_ratio = sos_n_perc[idx]
    ax.text(0.6, 0.19, f'{red_ratio:.2f}% ', transform=ax.transAxes, fontdict=fontdict4, va='center')
    ax.text(0.8, 0.19, f'(正值)', transform=ax.transAxes, fontdict=fontdict6, va='center')
    # Blue rectangle for < 0
    blue_patch = patches.Rectangle((0.5, 0.10), 0.075, 0.035, transform=ax.transAxes, color='blue', zorder=5)
    ax.add_patch(blue_patch)
    ax.text(0.6, 0.11, f'{blue_ratio:.2f}% ', transform=ax.transAxes, fontdict=fontdict4, va='center')
    ax.text(0.8, 0.11, f'(负值)', transform=ax.transAxes, fontdict=fontdict6, va='center')
 
print('drawing colorbar')
# 设置色标位置和大小
colorbar_width = 0.015
colorbar_height = 0.8
ax_colorbar = fig.add_axes([colorbar_left_adjusted-0.08, bottom_margin_adjusted, colorbar_width, colorbar_height])
# 创建色标
sm = mpl.cm.ScalarMappable(cmap='coolwarm', norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, cax=ax_colorbar)
# 设置色标刻度
cbar.set_ticks([-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8])
# 设置色标刻度标签格式为一位小数
cbar.ax.yaxis.set_major_formatter(FuncFormatter(lambda x, _: f'{x:.1f}'))
# 设置色标的刻度参数和移除边框
cbar.ax.tick_params(axis='y', which='both', length=tick_length, width=tick_width, pad=label_pad)
cbar.outline.set_linewidth(0)  # 移除色标的边框
# 设置色标标签
cbar.set_label('偏相关系数', fontdict=fontdict5)
 
plt.show()

完结撒花~欢迎交流和指正！

PS:图中的操作用Arcgis好像很快/(ㄒoㄒ)/~~