期刊复现 | Python 实现多组线性回归拟合图（含误差棒、渐变色标与置信区间）

论文地址：

https://www.nature.com/articles/s41467-025-63556-2

论文题目：

Amplified local cooling effect of forestation in warming Europe

复现图片:

图 c 解读：该图展示不同下垫面对应的参数变化趋势拟合结果，三组散点及误差棒分别代表林地、开阔地与二者差值数据，搭配线性拟合曲线可直观呈现变量间相关关系、斜率与相关系数，反映数据变化规律、误差范围及相关性强弱。正式图注：图 c 不同下垫面下参数趋势的线性拟合分析。

图 d 解读：该图为变量相关性散点图，散点颜色映射另一指标数值，误差棒表征数据测量误差，结合整体线性拟合曲线，可清晰体现两组变量的关联趋势、数据离散程度与整体相关特征。正式图注：图 d 变量间相关性及线性拟合结果。

图 e 解读：该图呈现单日不同时段地表温度差值与对应指标的分布及拟合关系，不同颜色散点与曲线分别对应夜间、日均值、日间三组数据，阴影区域表征拟合误差范围，可对比不同时段数据的分布特征、变化趋势与拟合精度。正式图注：图 e 不同时段地表温度差值与指标的拟合关系。

合并总图 解读：该图整合三组分析结果，依次展示参数趋势拟合、变量相关性、分时段拟合结果，全面呈现多组数据的分布规律、相关关系与变化特征，可整体对比不同维度下的数据表现。正式图注：图 1 综合分析结果汇总图。

完整代码：

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import os
import matplotlib
from PIL import Image

# 全局配置
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.size'] = 12
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42

COLOR_SCHEMES = {1: 'RdPu'}
SELECTED_SCHEME = 1
selected_cmap = COLOR_SCHEMES.get(SELECTED_SCHEME, 'RdPu')

SAVE_BASE_DIR = os.path.join(os.getcwd(), "图表")
if not os.path.exists(SAVE_BASE_DIR):
os.makedirs(SAVE_BASE_DIR)

LABEL_FONT_SIZE = 16

# 工具函数
def add_rotated_label(ax, x1, y1, x2, y2, text, color, txt_x, offset=0.02):
p1 = ax.transData.transform((x1, y1))
p2 = ax.transData.transform((x2, y2))
angle = np.degrees(np.arctan2(p2[1] – p1[1], p2[0] – p1[0]))
txt_y = (y2 – y1) / (x2 – x1) * txt_x + y1 + offset
ax.text(txt_x, txt_y, text, color=color, fontsize=12, fontweight='bold',
rotation=angle, rotation_mode='anchor', transform=ax.transData,
va='center', ha='left', zorder=10)

def add_rotated_text(ax, line, slope, intercept, r_value, color, txt_x, offset=0.005):
x1, y1 = line.get_xydata()[0]
x2, y2 = line.get_xydata()[-1]
p1 = ax.transData.transform((x1, y1))
p2 = ax.transData.transform((x2, y2))
dx = p2[0] – p1[0]
dy = p2[1] – p1[1]
angle = np.degrees(np.arctan2(dy, dx))
txt_y = slope * txt_x + intercept + offset
ax.text(txt_x, txt_y,
f"slope = {slope * 100:.2f}, r = {r_value:.2f}",
color=color, fontsize=12, fontweight='bold',
rotation=angle, rotation_mode='anchor',
transform=ax.transData, va='bottom', ha='left', zorder=10)

# 图表C
def plot_chart_c(data_path="data1.xlsx"):
df = pd.read_excel(data_path)
x = df['Trend of SC_w (%/decade)'].values
y_forest = df['Forest α_w (unitless/decade)'].values
err_forest = df['Forest α_w Error'].values
y_openland = df['Openland α_w (unitless/decade)'].values
err_openland = df['Openland α_w Error'].values
y_delta = df['Δα_w (unitless/decade)'].values
err_delta = df['Δα_w Error'].values

fig, ax = plt.subplots(figsize=(6, 5), dpi=150)

ax.errorbar(x, y_forest, yerr=err_forest, fmt='none', ecolor='#7ccd7c',
elinewidth=1, capsize=2, alpha=0.6)
ax.scatter(x, y_forest, color='#5cb85c', label='Forest α$_w$',
s=60, edgecolors='none', alpha=0.8)

ax.errorbar(x, y_openland, yerr=err_openland, fmt='none', ecolor='#fff68f',
elinewidth=1, capsize=2, alpha=0.6)
ax.scatter(x, y_openland, color='#f0e64c', label='Openland α$_w$',
s=60, edgecolors='none', alpha=0.8)

ax.errorbar(x, y_delta, yerr=err_delta, fmt='none', ecolor='#e6b8af',
elinewidth=1, capsize=2, alpha=0.6)
ax.scatter(x, y_delta, color='#cc8888', label='Δα$_w$',
s=60, edgecolors='none', alpha=0.8)

slope_f, intercept_f, r_f, _, _ = stats.linregress(x, y_forest)
x_fit = np.array([-17, 1])
y_fit_f = slope_f * x_fit + intercept_f
line_f, = ax.plot(x_fit, y_fit_f, color='#2e7d32', linewidth=2.5, zorder=5)

slope_o, intercept_o, r_o, _, _ = stats.linregress(x, y_openland)
y_fit_o = slope_o * x_fit + intercept_o
line_o, = ax.plot(x_fit, y_fit_o, color='#d4ac0d', linewidth=2.5, zorder=5)

slope_d, intercept_d, r_d, _, _ = stats.linregress(x, y_delta)
y_fit_d = slope_d * x_fit + intercept_d
line_d, = ax.plot(x_fit, y_fit_d, color='#b74141', linewidth=2.5, zorder=5)

ax.axhline(0, color='gray', linestyle='–', alpha=0.7, lw=1.5)
ax.axvline(0, color='gray', linestyle='–', alpha=0.7, lw=1.5)

ax.set_xlim(-17, 1)
ax.set_ylim(-0.10, 0.06)
ax.set_xticks(np.arange(-16, 2, 4))
ax.set_xlabel('Trend of SC$_w$ (%/decade)', fontsize=12)
ax.set_ylabel('Trend of α$_w$ (unitless/decade)', fontsize=12)
ax.tick_params(axis='both', labelsize=12)
ax.legend(loc='lower right', fontsize=12)
ax.text(-0.12, 1.02, 'c', transform=ax.transAxes, fontsize=LABEL_FONT_SIZE, fontweight='bold')

fig.canvas.draw()
add_rotated_text(ax, line_f, slope_f, intercept_f, r_f, '#2e7d32', txt_x=-16, offset=0.003)
add_rotated_text(ax, line_o, slope_o, intercept_o, r_o, '#d4ac0d', txt_x=-16, offset=0.003)
add_rotated_text(ax, line_d, slope_d, intercept_d, r_d, '#b74141', txt_x=-14, offset=0.003)

plt.tight_layout()
save_path = os.path.join(SAVE_BASE_DIR, 'chart_c.png')
plt.savefig(save_path, dpi=300, bbox_inches='tight')
plt.close()

# 图表D
def plot_chart_d(data_path="data2.xlsx"):
df_data = pd.read_excel(data_path)
x_data = df_data['x_data'].values
y_data = df_data['y_data'].values
y_error = df_data['y_error'].values
color_data = df_data['color_data'].values

cbar_ticks = [0, -5, -10, -15]
vmin_val = -15
vmax_val = 0

fig, ax = plt.subplots(figsize=(6, 5))
ax.axhline(0, color='gray', linestyle='–', alpha=0.7, lw=1.5)
ax.axvline(0, color='gray', linestyle='–', alpha=0.7, lw=1.5)

scatter = ax.scatter(x_data, y_data, c=color_data, cmap=selected_cmap + '_r',
s=120, edgecolors='black', linewidth=1.5, zorder=2,
vmin=vmin_val, vmax=vmax_val)

ax.errorbar(x_data, y_data, yerr=y_error, fmt='none',
ecolor='darkgray', elinewidth=1.5, capsize=3, alpha=0.8, zorder=1)

slope, intercept, r_value, p_value, std_err = stats.linregress(x_data, y_data)
x_fit = np.array([-0.005, 0.065])
y_fit = slope * x_fit + intercept
ax.plot(x_fit, y_fit, color='black', linewidth=2.5, zorder=3)

annotation_text = f"slope = {slope:.2f}\\nr = {r_value:.2f}"
ax.text(0.95, 0.95, annotation_text, transform=ax.transAxes,
fontsize=12, ha='right', va='top')

cbaxes = inset_axes(ax, width="40%", height="4%", loc='lower center',
bbox_to_anchor=(-0.25, 0.2, 1, 1), bbox_transform=ax.transAxes)
cbar = fig.colorbar(scatter, cax=cbaxes, orientation='horizontal', extend='neither')
cbar.outline.set_edgecolor('black')
cbar.outline.set_linewidth(1.5)
cbar.set_label('Trend of SC$_w$ (%/decade)', fontsize=12, labelpad=8)
cbar.ax.xaxis.set_label_position('top')
cbar.ax.tick_params(labelsize=12)
cbar.set_ticks(cbar_ticks)
cbar.ax.invert_xaxis()

ax.set_xlim(-0.01, 0.068)
ax.set_ylim(-0.5, 0.2)
ax.set_xlabel('Trend of Δα$_w$ (unitless/decade)', fontsize=12)
ax.set_ylabel('Trend of daytime ΔLST$_w$ (K/decade)', fontsize=12)
ax.tick_params(axis='both', which='major', labelsize=12)
ax.text(-0.1, 1.02, 'd', transform=ax.transAxes, fontsize=LABEL_FONT_SIZE, fontweight='bold')

save_path = os.path.join(SAVE_BASE_DIR, 'chart_d.png')
plt.savefig(save_path, dpi=300, bbox_inches='tight')
plt.close()

# 图表E
def plot_chart_e(data_path="data3.xlsx", save_filename="chart_e.png"):
if not os.path.exists(data_path):
raise FileNotFoundError(f"未找到Excel文件：{data_path}")

df = pd.read_excel(data_path)
required_columns = ['SC_w (%)', 'ΔLST_w (Nighttime)', 'ΔLST_w (Daily Mean)', 'ΔLST_w (Daytime)']
for col in required_columns:
if col not in df.columns:
raise ValueError(f"缺少列：{col}")

sc_w = df['SC_w (%)'].values
night_y = df['ΔLST_w (Nighttime)'].values
daily_y = df['ΔLST_w (Daily Mean)'].values
day_y = df['ΔLST_w (Daytime)'].values

night_reg = stats.linregress(sc_w / 100, night_y)
night_slope, night_intercept, night_slope_err = night_reg.slope, night_reg.intercept, night_reg.stderr

daily_reg = stats.linregress(sc_w / 100, daily_y)
daily_slope, daily_intercept, daily_slope_err = daily_reg.slope, daily_reg.intercept, daily_reg.stderr

day_reg = stats.linregress(sc_w / 100, day_y)
day_slope, day_intercept, day_slope_err = day_reg.slope, day_reg.intercept, day_reg.stderr

fig, ax = plt.subplots(figsize=(6, 5), dpi=150)

ax.scatter(sc_w, night_y, color='#ff9999', label='nighttime', s=120, alpha=0.7, edgecolors='none')
ax.scatter(sc_w, daily_y, color='#b3b3b3', label='daily mean', s=120, alpha=0.7, edgecolors='none')
ax.scatter(sc_w, day_y, color='#9999ff', label='daytime', s=120, alpha=0.7, edgecolors='none')

x_fit = np.linspace(sc_w.min(), sc_w.max(), 100)

y_n = night_slope * (x_fit / 100) + night_intercept
band_n = 0.03 + 0.001 * x_fit
ax.plot(x_fit, y_n, color='#ff3333', linewidth=2.5, zorder=5)
ax.fill_between(x_fit, y_n – band_n, y_n + band_n, color='#ffcccc', alpha=0.4)

y_d = daily_slope * (x_fit / 100) + daily_intercept
band_d = 0.04 + 0.0015 * x_fit
ax.plot(x_fit, y_d, color='#333333', linewidth=2.5, zorder=5)
ax.fill_between(x_fit, y_d – band_d, y_d + band_d, color='#dddddd', alpha=0.4)

y_t = day_slope * (x_fit / 100) + day_intercept
band_t = 0.06 + 0.002 * x_fit
ax.plot(x_fit, y_t, color='#3366ff', linewidth=2.5, zorder=5)
ax.fill_between(x_fit, y_t – band_t, y_t + band_t, color='#ccccff', alpha=0.4)

fig.canvas.draw()
add_rotated_label(ax, sc_w.min(), night_intercept, sc_w.max(), night_slope * (sc_w.max() / 100) + night_intercept,
f"slope = {night_slope:.2f} ± {night_slope_err:.2f}", '#ff3333', txt_x=sc_w.min() + 5,
offset=0.02)
add_rotated_label(ax, sc_w.min(), daily_intercept, sc_w.max(), daily_slope * (sc_w.max() / 100) + daily_intercept,
f"slope = {daily_slope:.2f} ± {daily_slope_err:.2f}", '#333333', txt_x=sc_w.min() + 5,
offset=0.02)
add_rotated_label(ax, sc_w.min(), day_intercept, sc_w.max(), day_slope * (sc_w.max() / 100) + day_intercept,
f"slope = {day_slope:.2f} ± {day_slope_err:.2f}", '#3366ff', txt_x=sc_w.min() + 5, offset=0.02)

ax.set_xlim(sc_w.min() – 1, sc_w.max() + 1)
ax.set_ylim(-0.3, 0.4)
ax.set_xlabel('SC$_w$ (%)', fontsize=12)
ax.set_ylabel('ΔLST$_w$ (K)', fontsize=12)
ax.tick_params(axis='both', labelsize=12)
ax.legend(loc='lower right', fontsize=12, frameon=True)
ax.text(-0.1, 1.02, 'e', transform=ax.transAxes, fontsize=LABEL_FONT_SIZE, fontweight='bold')

plt.tight_layout()
save_path = os.path.join(SAVE_BASE_DIR, save_filename)
plt.savefig(save_path, dpi=300, bbox_inches='tight')
plt.close()

# 拼接总图
def merge_three_plots():
p1 = os.path.join(SAVE_BASE_DIR, "chart_c.png")
p2 = os.path.join(SAVE_BASE_DIR, "chart_d.png")
p3 = os.path.join(SAVE_BASE_DIR, "chart_e.png")

im1 = Image.open(p1)
im2 = Image.open(p2)
im3 = Image.open(p3)

height = im1.height
im2 = im2.resize((int(im2.width * height / im2.height), height))
im3 = im3.resize((int(im3.width * height / im3.height), height))

total_width = im1.width + im2.width + im3.width
new_im = Image.new('RGB', (total_width, height), (255, 255, 255))

new_im.paste(im1, (0, 0))
new_im.paste(im2, (im1.width, 0))
new_im.paste(im3, (im1.width + im2.width, 0))

merge_path = os.path.join(SAVE_BASE_DIR, "合并总图.png")
new_im.save(merge_path, dpi=(300, 300))

# 主入口
if __name__ == "__main__":
DATA1_PATH = "data1.xlsx"
DATA2_PATH = "data2.xlsx"
DATA3_PATH = "data3.xlsx"

try:
plot_chart_c(DATA1_PATH)
plot_chart_d(DATA2_PATH)
plot_chart_e(DATA3_PATH)
merge_three_plots()
except Exception as e:
pass

数据获取

评论+私信获取

VX获取：期刊复现 | Python 实现多组线性回归拟合图（含误差棒、渐变色标与置信区间）