【气象常用】时间序列的线性拟合

就是99

于 2024-05-16 10:00:00 发布

阅读量745

点赞数 21

分类专栏：气象常用文章标签： python 数据分析

本文链接：https://blog.csdn.net/weixin_46921265/article/details/138805928

版权

气象常用专栏收录该内容

16 篇文章 1 订阅

订阅专栏

效果图：

主要步骤：

1. 数据准备：下载Hadley Centre observations datasets的HadSST数据

可参考【气象常用】时间序列图-CSDN博客

2. 数据处理：计算线性拟合

3. 图像绘制：绘制折线及拟合线，并添加文本

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

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

###############################################################################
# 导入库并设置字体
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import linregress
# 设置西文字体为新罗马字体，中文为宋体，字号为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

步骤二：从csv文件读入数据

###############################################################################
datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'

# 读取csv文件，并选择指定的列
datan = pd.read_csv(datapath+'HadSST.4.0.1.0_annual_NHEM.csv', 
                   usecols=['year', 'anomaly', 'lower_bound_95pct_bias_uncertainty_range', 'upper_bound_95pct_bias_uncertainty_range'])

xn = datan['year']
yn = datan['anomaly']

步骤三：计算拟合线

###############################################################################
# 拟合线计算
k1, b1, r_value1, p_value1, std_err1 = linregress(xn[0:51], yn[0:51])
k2, b2, r_value2, p_value2, std_err2 = linregress(xn[51:121], yn[51:121])
k3, b3, r_value3, p_value3, std_err3 = linregress(xn[121:], yn[121:])

# 根据回归系数生成拟合直线的x值
x_fit1 = np.linspace(1850, 1901, 51)
x_fit2 = np.linspace(1901, 1970, 69)
x_fit3 = np.linspace(1971, 2025, 54)

# 计算拟合直线的y值
y_fit1 = k1 * x_fit1 + b1
y_fit2 = k2 * x_fit2 + b2
y_fit3 = k3 * x_fit3 + b3

步骤四：绘制图像并添加文本

###############################################################################
#绘制图像
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.1, 0.6, 0.6, 0.4])

# 绘制海温折线
ax.plot(xn, yn, 'k-', linewidth=1,  )
ax.plot(x_fit1, y_fit1, 'k--', linewidth=1, )
ax.plot(x_fit2, y_fit2, 'k--', linewidth=1, )
ax.plot(x_fit3, y_fit3, 'k--', linewidth=1, )

# 添加数值
ax.text(1850, -0.6 , 'y = %7.4f * x + %7.4f' %(k1, abs(b1)), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.6 , 'y = %7.4f * x - %7.4f' %(k2, abs(b2)), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.6 , 'y = %7.4f * x - %7.4f' %(k3, abs(b3)), color='k', fontsize=10, rotation=0, horizontalalignment="left")

ax.text(1850, -0.7 , 'r_value：%7.4f' %(r_value1), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.7 , 'r_value：%7.4f' %(r_value2), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.7 , 'r_value：%7.4f' %(r_value3), color='k', fontsize=10, rotation=0, horizontalalignment="left")

ax.text(1850, -0.8 , 'std_err：%7.4f' %(std_err1), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.8 , 'std_err：%7.4f' %(std_err2), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.8 , 'std_err：%7.4f' %(std_err3), color='k', fontsize=10, rotation=0, horizontalalignment="left")

# 设置坐标轴信息
# xlim，ylim：设置x轴y轴显示得最大最小值
# xticks，yticks： 设置坐标轴刻度，可控制隔几个标注，本图为两个一标
# xticklabels，yticklabels：设置坐标轴刻度名称，维度要和xticks，yticks一致奥
# xlabel，ylabel：坐标轴名称
# title：图名
xlabels = np.arange(1850, 2025, 20)
ylabels = np.arange(-12, 13, 4)/10
ax.set_xticks(xlabels)
ax.set_xticklabels(xlabels, rotation = 0, fontsize = 12)
ax.set_yticks(ylabels)
ax.set_yticklabels(ylabels, fontsize = 12)
ax.set_xlabel('Year',fontsize=12)
ax.set_ylabel('Temperature Anomaly(℃)',fontsize=12)
ax.set_title('Sea Surface Temperature (HadSST4) Anomalies Over the Last 100 Years', fontsize=14)


# 设置水平检验线
plt.axhline(0, color='grey', linestyle='--')

步骤五：保存图像

###############################################################################
# 保存图像
plt.savefig(figpath+'208 时间序列线性拟合.png', bbox_inches = 'tight', dpi=600, format='png')
plt.show()

完整代码在这里：

###############################################################################
# 导入库并设置字体
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import linregress
# 设置西文字体为新罗马字体，中文为宋体，字号为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/'

# 读取csv文件，并选择指定的列
datan = pd.read_csv(datapath+'HadSST.4.0.1.0_annual_NHEM.csv', 
                   usecols=['year', 'anomaly', 'lower_bound_95pct_bias_uncertainty_range', 'upper_bound_95pct_bias_uncertainty_range'])

xn = datan['year']
yn = datan['anomaly']
###############################################################################
# 拟合线计算
k1, b1, r_value1, p_value1, std_err1 = linregress(xn[0:51], yn[0:51])
k2, b2, r_value2, p_value2, std_err2 = linregress(xn[51:121], yn[51:121])
k3, b3, r_value3, p_value3, std_err3 = linregress(xn[121:], yn[121:])

# 根据回归系数生成拟合直线的x值
x_fit1 = np.linspace(1850, 1901, 51)
x_fit2 = np.linspace(1901, 1970, 69)
x_fit3 = np.linspace(1971, 2025, 54)

# 计算拟合直线的y值
y_fit1 = k1 * x_fit1 + b1
y_fit2 = k2 * x_fit2 + b2
y_fit3 = k3 * x_fit3 + b3
###############################################################################
#绘制图像
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.1, 0.6, 0.6, 0.4])

# 绘制海温折线
ax.plot(xn, yn, 'k-', linewidth=1,  )
ax.plot(x_fit1, y_fit1, 'k--', linewidth=1, )
ax.plot(x_fit2, y_fit2, 'k--', linewidth=1, )
ax.plot(x_fit3, y_fit3, 'k--', linewidth=1, )

# 添加数值
ax.text(1850, -0.6 , 'y = %7.4f * x + %7.4f' %(k1, abs(b1)), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.6 , 'y = %7.4f * x - %7.4f' %(k2, abs(b2)), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.6 , 'y = %7.4f * x - %7.4f' %(k3, abs(b3)), color='k', fontsize=10, rotation=0, horizontalalignment="left")

ax.text(1850, -0.7 , 'r_value：%7.4f' %(r_value1), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.7 , 'r_value：%7.4f' %(r_value2), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.7 , 'r_value：%7.4f' %(r_value3), color='k', fontsize=10, rotation=0, horizontalalignment="left")

ax.text(1850, -0.8 , 'std_err：%7.4f' %(std_err1), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1920, -0.8 , 'std_err：%7.4f' %(std_err2), color='k', fontsize=10, rotation=0, horizontalalignment="left")
ax.text(1990, -0.8 , 'std_err：%7.4f' %(std_err3), color='k', fontsize=10, rotation=0, horizontalalignment="left")

# 设置坐标轴信息
# xlim，ylim：设置x轴y轴显示得最大最小值
# xticks，yticks： 设置坐标轴刻度，可控制隔几个标注，本图为两个一标
# xticklabels，yticklabels：设置坐标轴刻度名称，维度要和xticks，yticks一致奥
# xlabel，ylabel：坐标轴名称
# title：图名
xlabels = np.arange(1850, 2025, 20)
ylabels = np.arange(-12, 13, 4)/10
ax.set_xticks(xlabels)
ax.set_xticklabels(xlabels, rotation = 0, fontsize = 12)
ax.set_yticks(ylabels)
ax.set_yticklabels(ylabels, fontsize = 12)
ax.set_xlabel('Year',fontsize=12)
ax.set_ylabel('Temperature Anomaly(℃)',fontsize=12)
ax.set_title('Sea Surface Temperature (HadSST4) Anomalies Over the Last 100 Years', fontsize=14)


# 设置水平检验线
plt.axhline(0, color='grey', linestyle='--')
###############################################################################
# 保存图像
plt.savefig(figpath+'208 时间序列线性拟合.png', bbox_inches = 'tight', dpi=600, format='png')
plt.show()