本文介绍如何使用Python和netCDF库读取并处理哈德利中心的全球海冰和海面温度数据HadISST，包括数据预处理（异常值和缺失值处理），以及使用多层感知机（MLP）进行时间序列预测，展示了模型训练过程和结果可视化。
人工智能大作业，同学们可以拿去参考学习，或者作为课设。
问题：
简介
哈德利中心全球海冰和海面温度（HadISST）是 1871 年至今的每月全球完整 SST 和海冰浓度场的组合。
背景描述
哈德利中心全球海冰和海面温度（HadISST）是 1871 年至今的每月全球完整 SST 和海冰浓度场的组合。
时间范围
1870年1月——2020年12月
时间分辨率
Monthly 逐月
空间范围
Global 全球
原始HadISST数据的经度范围为 -180~180
而大部分气象要素数据的经度范围为 0~360
为便于进行数据处理与分析，该数据集已将经度范围调整为 0~360
空间分辨率
1° x 1°
数据格式
netcdf (.nc)
分析：
nc文件的读取
nc文件不同与普通的csv，对nc文件的分析我用了两种办法：
1. 简单分析
from netCDF4 import Dataset
sst_file = '../HadISST_sst.nc'
dataset = Dataset(sst_file)
# 查看文件中的变量名
print(dataset.variables)
这里用netCDF4读取，并用variables方法直接展示变量，得到下列结果：
OrderedDict([('time', <class 'netCDF4._netCDF4.Variable'>
float32 time(time)
standard_name: time
long_name: Time
units: days since 1870-1-1 0:0:0
calendar: gregorian
axis: T
unlimited dimensions: time
current shape = (1812,)
filling on, default _FillValue of 9.969209968386869e+36 used), ('lon', <class 'netCDF4._netCDF4.Variable'>
float64 lon(lon)
standard_name: longitude
long_name: longitude
units: degrees_east
axis: X
unlimited dimensions:
current shape = (360,)
filling on, default _FillValue of 9.969209968386869e+36 used), ('lat', <class 'netCDF4._netCDF4.Variable'>
float64 lat(lat)
standard_name: latitude
long_name: latitude
units: degrees_north
axis: Y
unlimited dimensions:
current shape = (180,)
filling on, default _FillValue of 9.969209968386869e+36 used), ('time_bnds', <class 'netCDF4._netCDF4.Variable'>
float32 time_bnds(time, nv)
unlimited dimensions: time
current shape = (1812, 2)
filling on, default _FillValue of 9.969209968386869e+36 used), ('sst', <class 'netCDF4._netCDF4.Variable'>
float32 sst(time, lat, lon)
standard_name: sea_surface_temperature
long_name: sst
units: C
_FillValue: -1e+30
missing_value: -1e+30
cell_methods: time: lat: lon: mean
unlimited dimensions: time
current shape = (1812, 180, 360)
filling on)])
Process finished with exit code 0
结果我们交直接给chatgpt解释，可以得到有效信息：
.nc文件包含以下变量：
- time: 时间变量，类型为float32，存储了从1870年1月1日起每天的时间（单位为“天”）。
- lon: 经度变量，类型为float64，存储了网格点的经度信息（单位为“度”）。
- lat: 纬度变量，类型为float64，存储了网格点的纬度信息（单位为“度”）。
- time_bnds: 时间边界变量，类型为float32，存储了每个时间步的开始和结束时间。
- sst: 海表温度变量，类型为float32，存储了每个时间步、每个网格点的海表温度数据（单位为“摄氏度”），缺失值为-1e+30。
也就是说，time、lon、lat、time_bnds是文件的一些介绍，真正的海温数据在sst中，通过以下代码读到sst的数据，之后部分暂时不作介绍了，因为这种方法较为原始。
sst_data = data['sst'][:]下面详细介绍用专业软件分析的方法：
2. 专业软件分析
软件：panoply
地址：没找到官网，可以自己百度搜索，我也直接从非官网处下载的
运行：直接运行这个exe就可以啦
打开数据集后是以下画面：
可以看到各维度的shape，sst是二维地图，双击进入sst：
选择用这种方式查看（注意x和y的选择不要反了）
打开后显示第一个时间跨度的全球海温：
左边一列可以选择年份月份，查看其他时间的海温。
可以看到数据集中有明显的异常值和缺失值，导致全球海温都是红色，北极部分应该有很多异常值，海温不可能达到-1000度左右，大陆部分是灰色，缺失值。待会儿数据预处理时候要注意。
对于异常值：
上图把Plot，切换成Array 1就能看到具体数据，大概浏览了一下数据集发现正常海温都有个区间，于是自己设置了一个阈值：-10到60度是正常温度，如果某处海温不在该范围，则替换为所在纬度正常值的平均值（海温在纬度上相差不多）。
# 数据预处理
for i in range(sst_data.shape[0]):
# 遍历每个纬度
for j in range(sst_data.shape[1]):
# 获取当前纬度的数据
lat_data = sst_data[i, j, :]
# 发现了异常值
if np.any((lat_data > 60) | (lat_data < -10)):
# 计算所在纬度正常值的平均值
mean_value = np.mean(lat_data[(lat_data <= 60) & (lat_data >= -10)])
# 替换异常值为平均值
lat_data[(lat_data > 60) | (lat_data < -10)] = mean_value
# 将处理后的数据赋值回原数组
sst_data[i, j, :] = lat_data
print("数据预处理完成")对于缺失值：
由于数据集里大陆全部当作缺失值了，这个值如果带进模型训练的计算，会成为一个无穷大值，导致数据溢出。
这里我选择把缺失值忽略掉再进行模型训练：
for epoch in range(total_epochs):
print(f"开始第 {epoch + 1}/{total_epochs} 次迭代")
for i in range(len(X_train)):
x = X_train[i]
y = y_train[i]
if np.isnan(x).any() or np.isnan(y).any():
continue # 如果数据中包含NaN值，则跳过这组数据
else:
model.partial_fit([x], [y]) # 针对每一对x, y进行部分训练
print("fit完成")
# 过滤NaN数据
X_test_filtered = []
y_test_filtered = []
for i in range(len(X_test)):
x = X_test[i]
y = y_test[i]
if np.isnan(x).any() or np.isnan(y).any():
continue
else:
X_test_filtered.append(x)
y_test_filtered.append(y)
y_pred = model.predict(X_test_filtered)
print("pre完成")
mse = mean_squared_error(y_test_filtered, y_pred)
mse_values.append(mse)
elapsed_time = time.time() - start_time
remaining_time = (total_epochs - epoch - 1) * elapsed_time / (epoch + 1)
print(f"Epoch {epoch + 1}/{total_epochs} - MSE: {mse:.6f} - Elapsed Time: {elapsed_time:.2f}s - Remaining Time: {remaining_time:.2f}s")模型设计：
划分输入和输出
offset=1，根据上一个月海温预测下一个月，比如根据1999年10月的海温数据预测1999年11月的海温。
你也可以尝试换成12，即用今年某月的海温，预测明年的这个月的海温，比如根据1999年10月的海温预测2000年10月的海温
# 划分x和y
offset = 1
X_train = train_set[:-offset]
y_train = train_set[offset:]模型参数
model = MLPRegressor(hidden_layer_sizes=(50,30), activation='relu', solver='adam', max_iter=500, verbose=True)-
hidden_layer_sizes：隐藏层的大小，指定一个元组来设置每个隐藏层中神经元的数量。例如，(50, 30) 表示有两个隐藏层，第一个隐藏层有50个神经元，第二个隐藏层有30个神经元。
-
activation：激活函数，用于激活神经元的输出。此处使用 'relu' 激活函数，它在输入为负数时返回0，输入为正数时返回原值。
-
solver：优化器，用于优化权重和偏置项的损失函数。此处使用 'adam' 优化器，它根据损失函数的梯度自适应地调整学习率。
-
max_iter：最大迭代次数，即训练过程中权重和偏置项更新的最大次数。此处设置为500次。
-
verbose：是否打印详细信息。此处设置为True，表示在训练过程中会打印出详细的迭代信息。
代码
最后献出完整训练代码：
from sklearn.neural_network import MLPRegressor
from netCDF4 import Dataset
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
import joblib
import time
import matplotlib.pyplot as plt
# 读取数据集
sst_file = '../HadISST_sst.nc'
dataset = Dataset(sst_file)
# 读取海表温度数据
sst_data = dataset.variables['sst'][:] # 读取有效数据
# 数据预处理
for i in range(sst_data.shape[0]):
# 遍历每个纬度
for j in range(sst_data.shape[1]):
# 获取当前纬度的数据
lat_data = sst_data[i, j, :]
# 发现了异常值
if np.any((lat_data > 60) | (lat_data < -10)):
# 计算所在纬度正常值的平均值
mean_value = np.mean(lat_data[(lat_data <= 60) & (lat_data >= -10)])
# 替换异常值为平均值
lat_data[(lat_data > 60) | (lat_data < -10)] = mean_value
# 将处理后的数据赋值回原数组
sst_data[i, j, :] = lat_data
print("数据预处理完成")
# 数据归一化
scaler = MinMaxScaler(feature_range=(0, 1))
sst_data_normalized = scaler.fit_transform(sst_data.reshape(-1, 1)).reshape(sst_data.shape)
print("归一化完成")
# 将数据转换为适合输入的形式
samples, features = sst_data_normalized.shape[0], sst_data_normalized.shape[1] * sst_data_normalized.shape[2]
sst_data_reshaped = sst_data_normalized.reshape(samples, features)
# 确定训练集和验证集数量
split_ratio = 0.2
train_test_samples = int((1 - split_ratio) * sst_data_reshaped.shape[0])# (int)
# 划分训练集、测试集和验证集
xunlian=sst_data_reshaped[:train_test_samples] # 训练+测试
# 划分训练集和测试集
train_set = xunlian[:int(xunlian.shape[0]*0.8)] # 训练
test_set = xunlian[int(xunlian.shape[0]*0.8):] # 测试
# 划分x和y
offset = 1
X_train = train_set[:-offset]
y_train = train_set[offset:]
X_test = test_set[:-offset]
y_test = test_set[offset:]
# 训练模型
model = MLPRegressor(hidden_layer_sizes=(50,30), activation='relu', solver='adam', max_iter=500, verbose=True)
print("模型定义完成")
start_time = time.time()
mse_values = []
total_epochs = 10
for epoch in range(total_epochs):
print(f"开始第 {epoch + 1}/{total_epochs} 次迭代")
for i in range(len(X_train)):
x = X_train[i]
y = y_train[i]
if np.isnan(x).any() or np.isnan(y).any():
continue # 如果数据中包含NaN值，则跳过这组数据
else:
model.partial_fit([x], [y]) # 针对每一对x, y进行部分训练
print("fit完成")
# 过滤NaN数据
X_test_filtered = []
y_test_filtered = []
for i in range(len(X_test)):
x = X_test[i]
y = y_test[i]
if np.isnan(x).any() or np.isnan(y).any():
continue
else:
X_test_filtered.append(x)
y_test_filtered.append(y)
y_pred = model.predict(X_test_filtered)
print("pre完成")
mse = mean_squared_error(y_test_filtered, y_pred)
mse_values.append(mse)
elapsed_time = time.time() - start_time
remaining_time = (total_epochs - epoch - 1) * elapsed_time / (epoch + 1)
print(f"Epoch {epoch + 1}/{total_epochs} - MSE: {mse:.6f} - Elapsed Time: {elapsed_time:.2f}s - Remaining Time: {remaining_time:.2f}s")
# 保存模型
joblib.dump(model, '../models/mlp_model.pkl')
# 绘制MSE曲线
plt.plot(range(1, total_epochs+1), mse_values, marker='o')
plt.xlabel('Epoch')
plt.ylabel('Mean Squared Error (MSE)')
plt.title('Model Evaluation')
plt.show()
# 保存图像到文件中
plt.savefig('MLP_mses.png')文件读取路径和模型保存路径，请根据自己需要修改哦。
数据集和验证的代码放在了资料里，进入我的个人空间下载。
2606
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
