datawhale2023第三期夏令营AI4S之大气科学笔记

本篇内容是datawhale第三期夏令营AI4S之大气科学笔记-第一篇

一 机器学习baseline流程

1 导入必要的库

import numpy as np  # 导入NumPy库,用于数值计算和数组操作
import pandas as pd  # 导入Pandas库,用于数据处理和分析
import xarray as xr  # 导入Xarray库,用于处理标签化的多维数据
from collections import defaultdict, Counter  # 导入defaultdict和Counter类,用于创建默认字典和计数器
import xgboost as xgb  # 导入XGBoost库,用于梯度提升树模型
import lightgbm as lgb  # 导入LightGBM库,也用于梯度提升树模型
from catboost import CatBoostRegressor  # 导入CatBoostRegressor类,用于CatBoost回归模型
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold  # 导入三种交叉验证类,分别是分层K折交叉验证、K折交叉验证和分组K折交叉验证
from sklearn.metrics import mean_squared_error  # 导入均方误差函数
from joblib import Parallel, delayed  # 导入Parallel和delayed类,用于并行化计算任务
from tqdm.auto import tqdm  # 导入tqdm类,用于在循环中显示进度条
import sys, os, gc, argparse, warnings, torch, glob  # 导入其他常用模块,如sys、os、gc、argparse、warnings、torch、glob

warnings.filterwarnings('ignore')  # 设置警告过滤器,将忽略警告信息

# pip install importlib-metadata==4.13.0
# pip install zarr lightgbm catboost
# pip install xarray[complete]

2 数据加载

path = './weather'  # 数据集文件路径

def chunk_time(ds):
    # 获取数据集维度的字典并赋值给dims变量
    dims = {k: v for k, v in ds.dims.items()}
    # 在'time'维度上将数据集进行分块处理,每块包含一个时间步
    dims['time'] = 1
    ds = ds.chunk(dims)
    return ds

def load_dataset():
    ds = []  
    for y in range(2007, 2012):  # 遍历年份2007到2011
        data_name = os.path.join(path, f'weather_round1_train_{y}')  # 构建数据文件名
        x = xr.open_zarr(data_name, consolidated=True)  # 打开Zarr格式的数据集文件
        print(f'{data_name}, {x.time.values[0]} ~ {x.time.values[-1]}')  # 输出数据文件名和时间范围
        ds.append(x)  # 将打开的数据集添加到ds列表中
        
    ds = xr.concat(ds, 'time')  # 在'time'维度上连接所有数据集
    ds = chunk_time(ds)  # 对数据集进行时间分块处理
    print(ds)  #查看数据信息
    return ds

ds = load_dataset().x  # 调用load_dataset函数加载数据集,并获取其中的'x'变量

打印信息如下:从2007-2011年5年的数据总共有7304条。

# 这里打印的信息的每次加载的文件名字与时间范围,共5条
./weather\weather_round1_train_2007, 2007-01-01T00:00:00.000000000 ~ 2007-12-31T18:00:00.000000000
./weather\weather_round1_train_2008, 2008-01-01T00:00:00.000000000 ~ 2008-12-31T18:00:00.000000000
./weather\weather_round1_train_2009, 2009-01-01T00:00:00.000000000 ~ 2009-12-31T18:00:00.000000000
./weather\weather_round1_train_2010, 2010-01-01T00:00:00.000000000 ~ 2010-12-31T18:00:00.000000000
./weather\weather_round1_train_2011, 2011-01-01T00:00:00.000000000 ~ 2011-12-31T18:00:00.000000000
# 这里打印的是加载后的数据形式,可以看见每一列的维度
<xarray.Dataset>
Dimensions:  (channel: 70, lat: 161, lon: 161, time: 7304)
Coordinates:

  * channel  (channel) <U5 'z50' 'z100' 'z150' 'z200' ... 'u10' 'v10' 'msl' 'tp'
  * lat      (lat) float32 50.0 49.75 49.5 49.25 49.0 ... 10.75 10.5 10.25 10.0
  * lon      (lon) float32 100.0 100.2 100.5 100.8 ... 139.2 139.5 139.8 140.0
  * time     (time) datetime64[ns] 2007-01-01 ... 2011-12-31T18:00:00
Data variables:
x        (time, channel, lat, lon) float32 dask.array<chunksize=(1, 70, 161, 161), meta=np.ndarray>

3 构建训练样本

3.1 取x份数据

由于评测数据只给了两条历史数据来预测未来20条数据,因此这里构建训练样本时,前两条数据构建特征,后20条数据来构建目标值,组成多步预测。
且预测的要求是:必须包含20个step对应未来时刻如下:[6,12,18,24,30,36,42,48,54,60,66,72,78,84,90,96,102,108,114,120](小时)
每天有4条数据:[6, 12, 18, 24]。因此我们构建特征时,采用每天的后两条数据来预测未来的20个step刚好数据连续(这就是为什么idx= i * 24的原因)。

lats = ds.lat.values.tolist()  # 将数据集中的纬度值转换为Python列表,并赋值给lats变量
lons = ds.lon.values.tolist()  # 将数据集中的经度值转换为Python列表,并赋值给lons变量

# 对齐测试数据,每份数据应该是22个时刻
# 获取多份训练数据,滑动5次,每次24个时刻(6天),每次滑动提取最后22个时刻作为一份训练数据
train_data = []
for i in range(5):
    if i == 0:  # 如果是第一份训练数据
        print(-22, 0)  # 输出起始和结束索引
        train_data.append(ds[-22:, :, :, :].values)  # 提取最后22个时刻的数据,添加到训练数据列表中
    else:
        idx = i * 24  # 计算索引偏移量
        train_data.append(ds[-22-idx:-idx, :, :, :].values)  # 提取最后22个时刻的数据,添加到训练数据列表中
        print(-22-idx, -idx)  # 输出起始和结束索引

# 经纬度标识字段
latlon_vals = []
for i in range(161):  # 遍历纬度索引
    for j in range(161):  # 遍历经度索引
        latlon_vals.append([lats[i], lons[j]])  # 将每个纬度和经度对应的值添加到列表中
latlon_vals = np.array(latlon_vals)  # 将列表转换为NumPy数组,并赋值给latlon_vals变量

# 确定好x份训练数据后,接下来需要提取特征
for flag, dat in tqdm(enumerate(train_data)):  # 遍历每份训练数据,并显示进度条
    # 第一时刻特征
    first_feat = dat[0,:,:,:].reshape(70,161*161).transpose()  # 将第一时刻的数据转换为特征向量,reshape后的形状为(70, 161*161),然后进行转置得到(161*161, 70)
    # 第二时刻特征
    second_feat = dat[1,:,:,:].reshape(70,161*161).transpose()  # 将第二时刻的数据转换为特征向量,reshape后的形状为(70, 161*161),然后进行转置得到(161*161, 70)
    # 两个时刻差分特征
    diff_feat = (dat[1,:,:,:] - dat[0,:,:,:]).reshape(70,161*161).transpose()  # 将两个时刻的数据差分并转换为特征向量,reshape后的形状为(70, 161*161),然后进行转置得到(161*161, 70)
    
    # 构建训练样本
    tmp_dat = dat[2:,-5:,:,:]  # 从第3个时刻开始获取最后5个时刻的数据,形状为(20, 5, 161, 161)
    for i in range(20):  # 遍历20个时刻
        # 时间特征、经纬特征
        time_vals = np.array([i]*161*161).reshape(161*161,1)  # 创建一个维度为(161*161,1)的数组,每个元素为当前时刻的时间特征i
        sub_vals = np.concatenate((time_vals, latlon_vals), axis=1)  # 将时间特征和经纬度特征拼接在一起,形状为(161*161,3)
        
        # 添加历史特征,第一时刻特征、第二时刻特征、两个时刻差分特征
        sub_vals = np.concatenate((sub_vals, first_feat), axis=1)  # 将第一时刻特征拼接到子特征向量中,形状变为(161*161,73)
        sub_vals = np.concatenate((sub_vals, second_feat), axis=1)  # 将第二时刻特征拼接到子特征向量中,形状变为(161*161,143)
        sub_vals = np.concatenate((sub_vals, diff_feat), axis=1)  # 将两个时刻差分特征拼接到子特征向量中,形状变为(161*161,213)
        
        # 添加5个目标变量,形状变为(161*161,218)
        for j in range(5):
            var_vals = tmp_dat[i,j,:,:].reshape(161*161, 1)  # 将目标变量的数据转换为特征向量,形状为(161*161,1)
            sub_vals = np.concatenate((sub_vals, var_vals), axis=1)  # 将目标变量拼接到子特征向量中
        
        if i == 0 :
            all_vals = sub_vals  # 如果是第一个时刻,则将子特征向量赋值给all_vals
        else:
            all_vals = np.concatenate((all_vals, sub_vals), axis=0)  # 如果不是第一个时刻,则将子特征向量与all_vals进行垂直拼接
        
    if flag == 0:
        final_vals = all_vals  # 如果是第一份训练数据,则将all_vals赋值给final_vals
    else:
        final_vals = np.concatenate((final_vals, all_vals), axis=0)  # 如果不是第一份训练数据,则将all_vals与final_vals进行垂直拼接
print(final_vals.shape)

打印信息如下

# 这里取了5份数据,这是每一份数据的起始索引与结束索引
-22 0
-46 -24
-70 -48
-94 -72
-118 -96
# 这里是最后5份数据的总共数据量。由(161*161*20*5,218) 20表示未来的20个step、5表示取5份数据
(2592100, 218)

3.2 将取出的数据转化为DataFrame格式

train_df = pd.DataFrame(final_vals)  # 将提取的特征final_vals转换为DataFrame格式,并赋值给train_df
train_df.columns = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)] + ['t2m','u10','v10','msl','tp']
# 添加列名,分别为时间特征、经度特征、纬度特征、210个特征向量、目标变量t2m、u10、v10、msl和tp

cols = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)]
# 创建列名列表,包括时间特征、经度特征、纬度特征和210个特征向量列名

数据太大,只展示前5行看看:

timelatlonfeat_0feat_1feat_2feat_3feat_4feat_5feat_6feat_205feat_206feat_207feat_208feat_209t2mu10v10msltp
00.050.0100.00-4.514620-4.506903-4.288303-3.989979-3.628431-3.401888-3.131553-0.674277-0.088560-0.0237650.0733700.0-6.2113150.609452-0.3061122.303155-0.604822
10.050.0100.25-4.530884-4.521255-4.301620-4.008497-3.653694-3.426081-3.152735-0.418350-0.1109390.0414250.0565710.0-5.6083990.484209-0.1889612.248176-0.604822
20.050.0100.50-4.547148-4.533813-4.316594-4.028692-3.680758-3.450274-3.173910-0.384310-0.0873700.0756550.0693160.0-5.6054040.4860610.0806052.248359-0.604822
30.050.0100.75-4.563421-4.546361-4.328246-4.045526-3.704220-3.472449-3.192436-0.443774-0.0530890.0124270.0285760.0-5.9144480.9789300.0917162.149359-0.604822
40.050.0101.00-4.577876-4.558918-4.338231-4.062355-3.729483-3.496637-3.213619-0.446717-0.0283300.005887-0.0214240.0-6.0663761.0168970.3351952.054393-0.604822

5 rows × 218 columns

4 训练模型

def train_model(train_x, train_y, label, seed=2023):
    # 初始化交叉验证的变量和参数
    oof = np.zeros(train_x.shape[0])  # 存储交叉验证的预测结果
    kf = KFold(n_splits=5, shuffle=True, random_state=seed)  # 5折交叉验证
    cv_scores = []  # 存储每一折的验证集分数

    for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
        print('************************************ {} {}************************************'.format(str(i+1), str(seed)))
        trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
        #trn_x, trn_y, val_x, val_y:根据索引划分出当前折的训练集和验证集。

        # 设置CatBoost模型的参数
        params = {'learning_rate': 0.5,
                  'depth': 5,
                  'bootstrap_type':'Bernoulli',
                  'random_seed':2023,
                  'od_type': 'Iter',
                  'od_wait': 100,
                  'random_seed': 11,
                  'allow_writing_files': False,
                  'task_type':"GPU",
                  'devices':'0:1'}

        model = CatBoostRegressor(iterations=100, **params)  # 创建CatBoost回归模型
        model.fit(trn_x, trn_y, eval_set=(val_x, val_y),
                  cat_features=[],
                  use_best_model=True,
                  verbose=1)  # 在训练集上拟合模型并在验证集上进行评估
        # model.save_model(f'catboost_{label}_fold{i}.model')  # 保存模型
        model.save_model(f'model/catboost_{label}_fold{i}.model')  # 保存模型到指定路径
        val_pred  = model.predict(val_x)  # 在验证集上进行预测
        oof[valid_index] = val_pred  # 将验证集的预测结果填充到oof数组中

        cv_scores.append(np.sqrt(mean_squared_error(val_y, val_pred)))  # 计算并存储验证集的均方根误差
        
        if i == 0: #如果是第一折(i==0),计算特征重要性并保存到CSV文件。
            fea_ = model.feature_importances_
            fea_name = model.feature_names_
            fea_score = pd.DataFrame({'fea_name':fea_name, 'score':fea_})
            fea_score = fea_score.sort_values('score', ascending=False)
            fea_score.to_csv('feature_importances.csv', index=False)

        print(cv_scores) #打印当前折的验证集分数。
    return oof
    
oofs = []  # 存储所有特征的交叉验证预测结果

# 对每个特征进行训练和预测
for label in ['t2m', 'u10', 'v10', 'msl', 'tp']:
    print(f'==================== {label} ====================')
    
    # 调用train_model函数进行训练,并将返回的交叉验证预测结果存储到cat_oof变量中
    cat_oof = train_model(train_df[cols], train_df[label], label)
    
    # 将交叉验证预测结果添加到oofs列表中
    oofs.append(cat_oof)

训练使用的学习率比较大、周期比较短。为了快速得到结果,后续可以改进。每一次训练都保存了训练好的模型,用于后续的预测。
在这里插入图片描述

根据保存的特征重要性文件,可以看出time占比非常大,其次是feat_139。

5 测试集预测

5.1 构建预测函数

def inter_model(test_x, label, seed=2023):
    kf = KFold(n_splits=5, shuffle=True, random_state=seed)  # 5折交叉验证
    test = np.zeros(test_x.shape[0])  # 存储测试集的预测结果

    params = {'learning_rate': 0.5, 'depth': 6, 'bootstrap_type':'Bernoulli','random_seed':seed, 'metric_period': 300,
              'od_type': 'Iter', 'od_wait': 100, 'random_seed': 11, 'allow_writing_files': False, 'task_type': 'CPU',
              'task_type':"GPU", 'devices':'0:1'}

    for i in range(0, 5):
        # print('************************************ {} {}************************************'.format(str(i+1), str(seed)))
        model = CatBoostRegressor(**params)
        # model.load_model(f'catboost_{label}_fold{i}.model')
        model.load_model(f'model/catboost_{label}_fold{i}.model')  # 加载训练好的模型
        test_pred = model.predict(test_x)  # 在测试集上进行预测
        test += test_pred / kf.n_splits  # 将每一折的预测结果加权平均得到最终的测试集预测结果

    return test

5.2 构建测试样本

def get_test_data(dat):
    # 第一时刻特征
    first_feat = dat[0,:,:,:].reshape(70,161*161).transpose() # 将第一时刻的特征数据从(70, 161, 161)转换为(161*161, 70)

    # 第二时刻特征
    second_feat = dat[1,:,:,:].reshape(70,161*161).transpose() # 将第二时刻的特征数据从(70, 161, 161)转换为(161*161, 70)

    # 两个时刻差分特征
    diff_feat = (dat[1,:,:,:] - dat[0,:,:,:]).reshape(70,161*161).transpose() # 计算两个时刻特征之间的差分,并将结果从(70, 161, 161)转换为(161*161, 70)
    
    # 构建测试样本
    for i in range(20): # 20个时刻
        # 时间特征、经纬特征
        time_vals = np.array([i]*161*161).reshape(161*161,1) # 创建一个形状为(161*161, 1)的数组,每个元素都是i
        sub_vals = np.concatenate((time_vals, latlon_vals), axis=1) # 将时间特征和经纬度特征按列进行拼接,形状为(161*161, 3)
        
        # 添加历史特征,第一时刻特征、第二时刻特征、两个时刻差分特征
        sub_vals = np.concatenate((sub_vals, first_feat), axis=1) # 将第一时刻特征拼接到sub_vals中,形状为(161*161, 73)
        sub_vals = np.concatenate((sub_vals, second_feat), axis=1) # 将第二时刻特征拼接到sub_vals中,形状为(161*161, 143)
        sub_vals = np.concatenate((sub_vals, diff_feat), axis=1) # 将两个时刻差分特征拼接到sub_vals中,形状为(161*161, 213)
        
        if i == 0 :
            all_vals = sub_vals
        else:
            all_vals = np.concatenate((all_vals, sub_vals), axis=0) # 将sub_vals沿着行方向进行拼接
            
    df = pd.DataFrame(all_vals) # 将结果转换为DataFrame格式
    df.columns = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)] # 给DataFrame的列命名
    return df

5.3 生成结果文件

def get_parallel_result(file):
    # 加载输入数据
    input_data = torch.load(file)  # 加载保存在文件中的输入数据

    # 生成测试样本,并构建特征
    test_df = get_test_data(np.array(input_data))  # 根据输入数据生成测试样本并构建特征

    # 保存结果
    res = test_df[['time']]  # 创建一个只包含时间列的DataFrame用于保存结果

    # 模型推理,预测测试集结果
    for label in ['t2m', 'u10', 'v10', 'msl', 'tp']:
        # 通过模型推理预测测试集结果,inter_model是一个未给出的模型推理函数
        cat_pred = inter_model(test_df[cols], label)
        res[label] = cat_pred  # 将预测结果保存到res中

    # 转为提交格式
    for label in ['t2m', 'u10', 'v10', 'msl', 'tp']:
        for i in range(20):
            sub_vals = res[res['time'] == i][label].values.reshape(1, 161, 161)
            # 将时间为i的某个标签的预测结果转换为形状为(1, 161, 161)的数组

            if i == 0:
                all_vals = sub_vals
            else:
                all_vals = np.concatenate((all_vals, sub_vals), axis=0)
                # 将所有时间步的预测结果按时间序列连接起来,形状为(20, 161, 161)

        all_vals = all_vals.reshape(20, 1, 161, 161)

        if label == 't2m':
            final_vals = all_vals
        else:
            final_vals = np.concatenate((final_vals, all_vals), axis=1)
            # 将不同标签的预测结果按照通道连接起来,形状为(20, 5, 161, 161)

    final_vals = torch.tensor(final_vals)  # 将最终的结果转换为PyTorch张量
    save_file = file.split('/')[-1]  # 获取保存结果的文件名

    # torch.save(final_vals.half(), f'./{save_file}')  # 将结果以半精度浮点类型保存到指定路径中
    torch.save(final_vals.half(), f'output/{save_file}')  # 将结果以半精度浮点类型保存到指定路径中

注意:在保存结果文件时,由于使用环境是linux,所以使用/进行文件名分割,但是如果在windows下,最后读取的文件路径会改变,如下所示:

./weather/weather_round1_test/input\000.pt
./weather/weather_round1_test/input\001.pt
./weather/weather_round1_test/input\002.pt

最后一个变为了反斜杠\,如果仍然采用save_file = file.split('/')[-1] ,你读取到的是input\000.pt,如果在最后保存时,你的output目录下没有input目录,你将会报错,如果在windows下,可以改为save_file = file.split('\\')[-1],加两个反斜杠的目的是由于它本身是转义字符。

5.4 并行运行

# 并行处理指定目录下的文件
# 设置为32的前提是电脑能跑动,否则会卡死,跑不动可以设置为2、4、8、16等等
Parallel(n_jobs=32)( #创建一个并行处理的上下文,指定使用 32 个工作线程。

    # 对每个文件路径调用 get_parallel_result 函数进行处理
    delayed(get_parallel_result)(file_path)
    #使用 delayed 函数将 get_parallel_result 函数及其参数 file_path 封装起来,表示要对每个文件路径调用 get_parallel_result(file_path) 进行处理
    
    # 遍历指定目录下以 .pt 为后缀的文件路径
    for file_path in tqdm(glob.glob(f'{path}/weather_round1_test/input/*.pt'))
    #使用glob.glob函数获取指定目录下所有以 .pt 结尾的文件路径,并使用 tqdm 函数在进度条中显示遍历的进度。
)

6 结果输出

#将output文件夹压缩,如果这段代码无法生成output.zip,运行下面的代码
_ = !zip -r output.zip output/

#如果上面代码无响应,运行下面这段代码
import zipfile, os
path = './output/'  # 要压缩的文件路径
zipName = 'output.zip'  # 压缩后的zip文件路径及名称

# 创建一个新的zip文件
f = zipfile.ZipFile(zipName, 'w', zipfile.ZIP_DEFLATED)
#使用zipfile模块创建一个新的zip文件对象,指定为写模式('w')并采用ZIP_DEFLATED压缩算法。

# 遍历指定路径下的所有文件和文件夹
for dirpath, dirnames, filenames in os.walk(path): #使用os.walk函数遍历指定路径下的所有文件和文件夹,包括子文件夹
    for filename in filenames: #遍历每个文件夹中的文件名。
        print(filename)
        # 将文件添加到zip文件中
        f.write(os.path.join(dirpath, filename))

# 关闭zip文件
f.close()

7 总结

以上是采用机器学习方法跑完整个流程,整个过程大概2小时左右,真的跑的久,特别是刚开始运行在本地,环境是windows,导致出现了上面写的注意的错误,找错误又找了很久。运行在阿里云上也试过了,舍不得在阿里云上边运行,边看代码,总怕免费包用完,才用了本地一步步的看代码,看数据。跑完整个流程,基础了解了数据情况与运行流程,接下来朝着优化方向前进。

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值