【气象常用】纬圈平均时间序列

效果图：

主要步骤：

1. 数据准备：我用的CMIP6数据，下载方法看这里【数据下载】CMIP6数据下载-CSDN博客

2. 数据处理：双线性插值， 计算距平， 计算标准差

3. 图像绘制：绘制纬圈平均时间序列

详细代码：着急的直接拖到最后有完整代码

步骤一：导入库包及图片存储路径并设置中文字体为宋体，西文为新罗马（没有的库包要先下好奥）

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import xarray as xr
from cnmaps import get_adm_maps
from scipy.interpolate import griddata 
from scipy.interpolate import RegularGridInterpolator  
from matplotlib.colors import ListedColormap
# import pygrib as pg
# 设置西文字体为新罗马字体，中文为宋体，字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False

datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

步骤二：读入数据（这里也可以考虑定义函数）

##########################################################
# 读入数据
f1 = xr.open_dataset(datapath + 'pr_Amon_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')
f2 = xr.open_dataset(datapath + 'pr_Amon_NESM3_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')
f3 = xr.open_dataset(datapath + 'pr_Amon_CAMS-CSM1-0_ssp245_r1i1p1f1_gn_201501-209912_yearmean.nc')
f4 = xr.open_dataset(datapath + 'pr_Amon_FGOALS-g3_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')

# print(f_h)
pr1 = f1['pr'][:-1, :, :].values *60*60*24*365/1000*1000
pr2 = f2['pr'][:-1, :, :].values *60*60*24*365/1000*1000 # mm
pr3 = f3['pr'].values *60*60*24*365/1000*1000
pr4 = f4['pr'][:-1, :, :].values *60*60*24*365/1000*1000

time = f1['time']
lons1 = f1['lon'].values  
lats1 = f1['lat'].values  
lon1, lat1 = np.meshgrid(lons1, lats1)
lons2 = f2['lon'].values  
lats2 = f2['lat'].values  
lon2, lat2 = np.meshgrid(lons2, lats2)
lons3 = f3['lon'].values  
lats3 = f3['lat'].values  
lon3, lat3 = np.meshgrid(lons3, lats3)
lons4 = f4['lon'].values  
lats4 = f4['lat'].values  
lon4, lat4 = np.meshgrid(lons4, lats4)

步骤三：对数据进行双线性插值，并计算距平和标准差

##########################################################
# 数据处理
# 插值要用的新的格点
time0 = range(2015, 2100, 1)
lons0 = np.linspace(lons1.min(), lon1.max(), 80) 
lats0 = np.linspace(lats1.min(), lat1.max(), 40)  
lon0, lat0 = np.meshgrid(lons0, lats0)  
new_points = np.array([lat0.ravel(), lon0.ravel()]).T  
# 双线性插值
def chazhi(data, lat_old, lon_old, lat_new, lon_new):
    data_new = np.empty((85, 40, 80))
    for i in range(0, 85, 1):
        interpolator = RegularGridInterpolator((lat_old.ravel(), lon_old.ravel()), data[i, :, :], 
                                           method='linear', bounds_error=False, fill_value=None)
        data_new[i, :, :] = interpolator(new_points).reshape(lat0.shape)
    return data_new

pr1_new = chazhi(pr1, lats1, lons1, lats0, lons0)
pr2_new = chazhi(pr2, lats2, lons2, lats0, lons0)
pr3_new = chazhi(pr3, lats3, lons3, lats0, lons0)
pr4_new = chazhi(pr4, lats4, lons4, lats0, lons0)

#计算距平（这里算的没什么实际意义）
pr1_mean = np.mean(pr1_new, 2) - np.mean(pr1_new[0:6], (0,2))
pr2_mean = np.mean(pr2_new, 2) - np.mean(pr2_new[0:6], (0,2))
pr3_mean = np.mean(pr3_new, 2) - np.mean(pr3_new[0:6], (0,2))
pr4_mean = np.mean(pr4_new, 2) - np.mean(pr4_new[0:6], (0,2))
# 计算方差
pr = np.stack((pr1_mean, pr2_mean, pr3_mean, pr4_mean), axis=0)
variance = np.std(pr, axis=0)

步骤四：绘制图像

##########################################################
#绘制图像
fig = plt.figure(figsize=(15, 6))
level = range(-100, 101, 5)
x = range(0, 85, 1)
y = range(-90, 90, 1)

def plot_cont(location, data, text, textposition):
    
    ax = fig.add_axes(location)

    cs = ax.contourf(x, lats0, data.T, cmap='YlGn', levels=level, extend='both')
    
    ax.set(xticks=x[::20], xticklabels=time0[::20],
           yticks=y[::30], yticklabels=y[::30] )
    title = ax.set_title(text, fontsize=12)
    title.set_position(textposition)  # 调整位置
    return cs
    
# 调用上面的函数
c1 = plot_cont([0.1, 0.38, 0.3, 0.2], pr1_mean, '(b) BCC-CSM2-MR', [0.1, 1.05])
c2 = plot_cont([0.43, 0.38, 0.3, 0.2], pr2_mean, '(c) NESM3', [0.08, 1.05])
c3 = plot_cont([0.1, 0.1, 0.3, 0.2], pr3_mean, '(d) CAMS-CSM1-0',  [0.1, 1.05])
c4 = plot_cont([0.43, 0.1, 0.3, 0.2], pr4_mean, '(e) FGOALS-g3', [0.1, 1.05])
c5 = plot_cont([0.26, 0.66, 0.3, 0.2], 2 * variance, '(a) Model spread', [0.1, 1.05])
# 色标设置
rect = [0.23, 0.015, 0.35, 0.02]
cbar_ax = fig.add_axes(rect)

cb = fig.colorbar(c1, drawedges=True,
                  cax=cbar_ax, orientation='horizontal')
cb.set_label('mm', fontsize=12)
cb.ax.tick_params(length=0)

步骤五：保存图像

###############################################################################
# 存图
plt.savefig(figpath+'209 纬圈平均的时间序列', dpi=600, bbox_inches = 'tight')
plt.show()

完整代码在这里：

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import xarray as xr
from cnmaps import get_adm_maps
from scipy.interpolate import griddata 
from scipy.interpolate import RegularGridInterpolator  
from matplotlib.colors import ListedColormap
# import pygrib as pg
# 设置西文字体为新罗马字体，中文为宋体，字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False

datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'
##########################################################
# 读入数据
f1 = xr.open_dataset(datapath + 'pr_Amon_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')
f2 = xr.open_dataset(datapath + 'pr_Amon_NESM3_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')
f3 = xr.open_dataset(datapath + 'pr_Amon_CAMS-CSM1-0_ssp245_r1i1p1f1_gn_201501-209912_yearmean.nc')
f4 = xr.open_dataset(datapath + 'pr_Amon_FGOALS-g3_ssp245_r1i1p1f1_gn_201501-210012_yearmean.nc')

# print(f_h)
pr1 = f1['pr'][:-1, :, :].values *60*60*24*365/1000*1000
pr2 = f2['pr'][:-1, :, :].values *60*60*24*365/1000*1000 # mm
pr3 = f3['pr'].values *60*60*24*365/1000*1000
pr4 = f4['pr'][:-1, :, :].values *60*60*24*365/1000*1000

time = f1['time']
lons1 = f1['lon'].values  
lats1 = f1['lat'].values  
lon1, lat1 = np.meshgrid(lons1, lats1)
lons2 = f2['lon'].values  
lats2 = f2['lat'].values  
lon2, lat2 = np.meshgrid(lons2, lats2)
lons3 = f3['lon'].values  
lats3 = f3['lat'].values  
lon3, lat3 = np.meshgrid(lons3, lats3)
lons4 = f4['lon'].values  
lats4 = f4['lat'].values  
lon4, lat4 = np.meshgrid(lons4, lats4)
##########################################################
# 数据处理
# 插值要用的新的格点
time0 = range(2015, 2100, 1)
lons0 = np.linspace(lons1.min(), lon1.max(), 80) 
lats0 = np.linspace(lats1.min(), lat1.max(), 40)  
lon0, lat0 = np.meshgrid(lons0, lats0)  
new_points = np.array([lat0.ravel(), lon0.ravel()]).T  
# 双线性插值
def chazhi(data, lat_old, lon_old, lat_new, lon_new):
    data_new = np.empty((85, 40, 80))
    for i in range(0, 85, 1):
        interpolator = RegularGridInterpolator((lat_old.ravel(), lon_old.ravel()), data[i, :, :], 
                                           method='linear', bounds_error=False, fill_value=None)
        data_new[i, :, :] = interpolator(new_points).reshape(lat0.shape)
    return data_new

pr1_new = chazhi(pr1, lats1, lons1, lats0, lons0)
pr2_new = chazhi(pr2, lats2, lons2, lats0, lons0)
pr3_new = chazhi(pr3, lats3, lons3, lats0, lons0)
pr4_new = chazhi(pr4, lats4, lons4, lats0, lons0)

#计算距平（这里算的没什么实际意义）
pr1_mean = np.mean(pr1_new, 2) - np.mean(pr1_new[0:6], (0,2))
pr2_mean = np.mean(pr2_new, 2) - np.mean(pr2_new[0:6], (0,2))
pr3_mean = np.mean(pr3_new, 2) - np.mean(pr3_new[0:6], (0,2))
pr4_mean = np.mean(pr4_new, 2) - np.mean(pr4_new[0:6], (0,2))
# 计算标准差
pr = np.stack((pr1_mean, pr2_mean, pr3_mean, pr4_mean), axis=0)
variance = np.std(pr, axis=0)
##########################################################
#绘制图像
fig = plt.figure(figsize=(15, 6))
level = range(-100, 101, 5)
x = range(0, 85, 1)
y = range(-90, 90, 1)

def plot_cont(location, data, text, textposition):
    
    ax = fig.add_axes(location)

    cs = ax.contourf(x, lats0, data.T, cmap='YlGn', levels=level, extend='both')
    
    ax.set(xticks=x[::20], xticklabels=time0[::20],
           yticks=y[::30], yticklabels=y[::30] )
    title = ax.set_title(text, fontsize=12)
    title.set_position(textposition)  # 调整位置
    return cs
    
# 调用上面的函数
c1 = plot_cont([0.1, 0.38, 0.3, 0.2], pr1_mean, '(b) BCC-CSM2-MR', [0.1, 1.05])
c2 = plot_cont([0.43, 0.38, 0.3, 0.2], pr2_mean, '(c) NESM3', [0.08, 1.05])
c3 = plot_cont([0.1, 0.1, 0.3, 0.2], pr3_mean, '(d) CAMS-CSM1-0',  [0.1, 1.05])
c4 = plot_cont([0.43, 0.1, 0.3, 0.2], pr4_mean, '(e) FGOALS-g3', [0.1, 1.05])
c5 = plot_cont([0.26, 0.66, 0.3, 0.2], 2 * variance, '(a) Model spread', [0.1, 1.05])
# 色标设置
rect = [0.23, 0.015, 0.35, 0.02]
cbar_ax = fig.add_axes(rect)

cb = fig.colorbar(c1, drawedges=True,
                  cax=cbar_ax, orientation='horizontal')
cb.set_label('mm', fontsize=12)
cb.ax.tick_params(length=0)
# ###############################################################################
# 存图
plt.savefig(figpath+'209 纬圈平均的时间序列', dpi=600, bbox_inches = 'tight')
plt.show()