【论文复刻】温升阶梯图+拟合线

复刻了一下这篇论文里的fig1a：Mitigation of China’s carbon neutrality to global warming | Nature Communications

相较于原图没有绘制检验区间，没有加文本

效果图：

主要步骤：

1. 数据准备：我这里用的CMIP6数据，大家可以替换为自己的数据（CMIP6数据下载方法可以看这里【数据下载】CMIP6数据下载-CSDN博客）

2. 数据处理：对原数据进行了全球平均、年平均，然后计算了较工业革命前的温升，并对温升进行了二次拟合

3. 图像绘制：绘制了温升的阶梯图、拟合线图，对图例进行了精细化设置，并设置了1.5/2℃的水平线、阴影区间及坐标轴相关信息

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

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

###############################################################################
# 导入库及文件
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt

import netCDF4 as nc

from scipy.optimize import curve_fit

# 设置西文字体为新罗马字体，中文为宋体，字号为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/'

# 定义读入数据函数
def read_file(filename, var):
    file = Dataset(datapath+filename, mode='r')
    time = file.variables['time']
    time_new = nc.num2date(time[:], time.units)
    data = file.variables[var][:]
    file.close()
    return data

# 定义全球平均函数
def glo_mean(data):
    data_mean = np.mean(data, (1, 2))
    return data_mean

# 定义年平均函数
def year_mean(data,nyear):
    data_mean = np.empty(nyear)
    for i in range(0, nyear, 1):
        data_mean[i] = np.mean(data[i * 12:(i + 1) * 12], 0)
    return data_mean

# 这里是我的文件名称哈
path_his = 'tas_Amon_CAMS-CSM1-0_historical_r1i1p1f1_gn_185001-201412.nc'
path_585 = 'tas_Amon_CAMS-CSM1-0_ssp585_r1i1p1f1_gn_201501-209912.nc'
path_370 = 'tas_Amon_CAMS-CSM1-0_ssp370_r1i1p1f1_gn_201501-209912.nc'
path_245 = 'tas_Amon_CAMS-CSM1-0_ssp245_r1i1p1f1_gn_201501-209912.nc'
path_126 = 'tas_Amon_CAMS-CSM1-0_ssp126_r1i1p1f1_gn_201501-209912.nc'

# 计算工业革命前全球平均地表温度
t_his = np.mean(read_file(path_his, 'tas')[0:51*12])

# 利用 read_file 读入数据
# 利用 gal_mean 计算全球平均
# 利用 year_mean 计算年平均
t_585 = year_mean(glo_mean(read_file(path_585, 'tas')), 85)
t_370 = year_mean(glo_mean(read_file(path_370, 'tas')), 85)
t_245 = year_mean(glo_mean(read_file(path_245, 'tas')), 85)
t_126 = year_mean(glo_mean(read_file(path_126, 'tas')), 85)

# 计算升温情况
dt_585 = t_585 - t_his
dt_370 = t_370 - t_his
dt_245 = t_245 - t_his
dt_126 = t_126 - t_his

步骤三：计算拟合函数

x = np.arange(0, 85, 1)

# 定义拟合函数
def f(x, a, b, c):
    return a * x**2 + b * x + c
def fit(data):
    params, covariance = curve_fit(f, x, data)
    a, b, c = params
    y = f(x, a, b, c)  # 计算对应的y值
    return y

# 计算拟合线
dt_585_fit = fit(dt_585)
dt_370_fit = fit(dt_370)
dt_245_fit = fit(dt_245)
dt_126_fit = fit(dt_126)

步骤四：绘制阶梯图、拟合线图，设置图例、1.5/2℃的水平线、阴影区间及坐标轴相关信息

###############################################################################
# 绘制图像
colors = ['#F77D75', '#7CAD00', '#00BEC3', '#C67CFF']
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.1, 0.6, 0.6, 0.4])

# 绘制阶梯线
plt.step(x, dt_585, color=colors[0], linewidth=1, label='SSP 5-8.5')
plt.step(x, dt_370, color=colors[1], linewidth=1, label='SSP 3-7.0')
plt.step(x, dt_245, color=colors[2], linewidth=1, label='SSP 2-4.5')
plt.step(x, dt_126, color=colors[3], linewidth=1, label='SSP 1-2.6')

# 绘制拟合线
plt.plot(x, dt_585_fit, color=colors[0], linewidth=2, label='SSP 5-8.5_fit')
plt.plot(x, dt_370_fit, color=colors[1], linewidth=2, label='SSP 3-7.0_fit')
plt.plot(x, dt_245_fit, color=colors[2], linewidth=2, label='SSP 2-4.5_fit')
plt.plot(x, dt_126_fit, color=colors[3], linewidth=2, label='SSP 1-2.6_fit')

# 设置横纵坐标轴
xlabels = np.arange(2015, 2100, 5)
ylabels = np.arange(0, 6, 1)
ax.set(xlim=(0, 84),xticks=x[::5], xticklabels=xlabels,
       ylim=(0, 5), yticklabels=ylabels)

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('ΔT(℃)', fontsize=12)

# 设置图例
ax.legend(frameon=False, fontsize=12, loc=2, 
          ncol=2, columnspacing=5.5,
          bbox_to_anchor=(0.06, 0.99))

# 绘制检验线
plt.axhline(y=1.5, color='g', linestyle='--')
plt.axhline(y=2, color='r', linestyle='--')

# 绘制阴影
plt.axvspan(6, 25, color='purple', alpha=0.1)
plt.axvspan(26, 45, color='purple', alpha=0.1)
plt.axvspan(46, 65, color='purple', alpha=0.1)
plt.axvspan(66, 85, color='purple', alpha=0.1)

步骤五：保存图像

###############################################################################
# 保存图像
plt.savefig(figpath+'302 温升阶梯+拟合线.png', bbox_inches = 'tight', dpi=600, format='png')
plt.show()

完整代码在这里：

###############################################################################
# 导入库及文件
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt

import netCDF4 as nc

from scipy.optimize import curve_fit

# 设置西文字体为新罗马字体，中文为宋体，字号为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/'

# 定义读入数据函数
def read_file(filename, var):
    file = Dataset(datapath+filename, mode='r')
    time = file.variables['time']
    time_new = nc.num2date(time[:], time.units)
    data = file.variables[var][:]
    file.close()
    return data

# 定义全球平均函数
def glo_mean(data):
    data_mean = np.mean(data, (1, 2))
    return data_mean

# 定义年平均函数
def year_mean(data,nyear):
    data_mean = np.empty(nyear)
    for i in range(0, nyear, 1):
        data_mean[i] = np.mean(data[i * 12:(i + 1) * 12], 0)
    return data_mean

# 这里是我的文件名称哈
path_his = 'tas_Amon_CAMS-CSM1-0_historical_r1i1p1f1_gn_185001-201412.nc'
path_585 = 'tas_Amon_CAMS-CSM1-0_ssp585_r1i1p1f1_gn_201501-209912.nc'
path_370 = 'tas_Amon_CAMS-CSM1-0_ssp370_r1i1p1f1_gn_201501-209912.nc'
path_245 = 'tas_Amon_CAMS-CSM1-0_ssp245_r1i1p1f1_gn_201501-209912.nc'
path_126 = 'tas_Amon_CAMS-CSM1-0_ssp126_r1i1p1f1_gn_201501-209912.nc'

# 计算工业革命前全球平均地表温度
t_his = np.mean(read_file(path_his, 'tas')[0:51*12])

# 利用 read_file 读入数据
# 利用 gal_mean 计算全球平均
# 利用 year_mean 计算年平均
t_585 = year_mean(glo_mean(read_file(path_585, 'tas')), 85)
t_370 = year_mean(glo_mean(read_file(path_370, 'tas')), 85)
t_245 = year_mean(glo_mean(read_file(path_245, 'tas')), 85)
t_126 = year_mean(glo_mean(read_file(path_126, 'tas')), 85)

# 计算升温情况
dt_585 = t_585 - t_his
dt_370 = t_370 - t_his
dt_245 = t_245 - t_his
dt_126 = t_126 - t_his

x = np.arange(0, 85, 1)

# 定义拟合函数
def f(x, a, b, c):
    return a * x**2 + b * x + c
def fit(data):
    params, covariance = curve_fit(f, x, data)
    a, b, c = params
    y = f(x, a, b, c)  # 计算对应的y值
    return y

# 计算拟合线
dt_585_fit = fit(dt_585)
dt_370_fit = fit(dt_370)
dt_245_fit = fit(dt_245)
dt_126_fit = fit(dt_126)

###############################################################################
# 绘制图像
colors = ['#F77D75', '#7CAD00', '#00BEC3', '#C67CFF']
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.1, 0.6, 0.6, 0.4])

# 绘制阶梯线
plt.step(x, dt_585, color=colors[0], linewidth=1, label='SSP 5-8.5')
plt.step(x, dt_370, color=colors[1], linewidth=1, label='SSP 3-7.0')
plt.step(x, dt_245, color=colors[2], linewidth=1, label='SSP 2-4.5')
plt.step(x, dt_126, color=colors[3], linewidth=1, label='SSP 1-2.6')

# 绘制拟合线
plt.plot(x, dt_585_fit, color=colors[0], linewidth=2, label='SSP 5-8.5_fit')
plt.plot(x, dt_370_fit, color=colors[1], linewidth=2, label='SSP 3-7.0_fit')
plt.plot(x, dt_245_fit, color=colors[2], linewidth=2, label='SSP 2-4.5_fit')
plt.plot(x, dt_126_fit, color=colors[3], linewidth=2, label='SSP 1-2.6_fit')

# 设置横纵坐标轴
xlabels = np.arange(2015, 2100, 5)
ylabels = np.arange(0, 6, 1)
ax.set(xlim=(0, 84),xticks=x[::5], xticklabels=xlabels,
       ylim=(0, 5), yticklabels=ylabels)

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('ΔT(℃)', fontsize=12)

# 设置图例
ax.legend(frameon=False, fontsize=12, loc=2, 
          ncol=2, columnspacing=5.5,
          bbox_to_anchor=(0.06, 0.99))

# 绘制检验线
plt.axhline(y=1.5, color='g', linestyle='--')
plt.axhline(y=2, color='r', linestyle='--')

# 绘制阴影
plt.axvspan(6, 25, color='purple', alpha=0.1)
plt.axvspan(26, 45, color='purple', alpha=0.1)
plt.axvspan(46, 65, color='purple', alpha=0.1)
plt.axvspan(66, 85, color='purple', alpha=0.1)

###############################################################################
# 保存图像
plt.savefig(figpath+'302 温升阶梯+拟合线.png', bbox_inches = 'tight', dpi=600, format='png')
plt.show()