AI夏令营天气预报大赛

AI夏令营第三期

比赛任务

提供10年的再分析数据，通过输入历史70个大气变量数据，预测华东区域未来1-5天的5个地面变量。

赛题数据集

==输入历史2个时刻的多个大气变量，输出未来1-5天每6小时的5个地面变量==

推理文件的输入Input: (2 x 70 x H x W)，输出Output: (20 X 5 X H X W)

==疑问：这个输入输出什么意思？矩阵的意思吗？==

其数据类型是float16（半精度的一个数据类型），存储为pt格式。==可以不使用70个变量作为输入，但是输出必须包含5个地面变量==，且顺序必须按照【T2M（2米温度），U10（10米纬向风），V10（10米径向风）， MSL（平均海平面气压），TP（6小时累计降水）】，必须包含20个step对应未来时刻：[6,12,18,24,30,36,42,48,54,60,66,72,78,84,90,96,102,108,114,120]（小时）

深度学习的宏观框架 训练和推理

训练

打个比方，想要训练一个能区分橙子和苹果的模型，你需要搜索一下苹果和橙子的图片，这些图片放在一起称为训练数据集（training dataset），训练数据集是有标签的，苹果的图片标签就是苹果，橙子依然。一个初始的神经网络通过不断地优化自身参数，来让自己变得准确，可能开始10个苹果图片，只有5张被认为是苹果，另外5张认错了，这个时候通过优化参数让另外5张也被认为是对的。这个过程称之为训练。

疑问：什么是神经网络？如何优化参数？（这些对于比赛没影响有影响的话我再补充）

推理

你训练好了一个模型，在训练数据集中表现良好，我们期望她可以对没看过的图片进行识别，你重新拍一张图片扔进网络让网络做判断，这种图片叫做现场数据，如果现场数据的区分准确率非常高，那么证明网络训练是非常好的，我们把训练好的模型拿出来遛一遛的过程称为推理

评价指标

指标还没搞懂，后面继续更新

解题思路

这个问题属于典型的回归问题

针对这类时间序列预测问题方法比较灵活，传统的时序模型、机器学习、深度学习方法均可以使用。

1、统计策略：使用最近时刻的结果进行均值、中位数、时间衰减等方式直接统计得到未来结果，这种方式比较简单，可以快速得到结果；

2、时序模型：比较常用的方法有指数平滑法、灰色预测模型、ARIMA预测、季节Sarima模型、VAR模型等，仅能刻画序列信息，无法加入其他信息进行训练，比如离散类特征；

3、机器学习模型：常见的为lightgbm、xgboost、catboost，需要构建大量时序相关特征；

4、深度学习模型：常见为rnn、lstm、cnn、transformer这类模型，可以直接输入序列信息，不需要构建大量的人工特征；

因为整体数据处理起来还是比较复杂的，为了快速的得到一个结果，所以我们选择使用统计策略作为baseline方案，后续进阶实践部分则使用机器学习模型。不过我们会在baseline实践部分花些时间进行数据解读，为进阶部分做铺垫。

Stept1 搞懂统计方法的BasLine

1、导入模块（导入工具包）

import os  #os库是Python的标准库
import numpy as np #是一个数学基础库，可以计算矩阵，导入numpy库将其命名为np
import pandas as pd #导入的是一个时间序列数据相关的库命名为pd
import xarray as xr #Python 模块，支持带有维度、坐标和属性标注的多维数组

import torch # 导入的是pytorch深度学习框架
torch.random.seed() #随机数
np.random.seed(0) #随机数

2、数据探索（将数据分块打印出来看看长啥样）

数据探索性分析，了解数据集，了解变量间的项目关系以及变量与预测值之间的关系，对已有数据在少量的先验假设下，通过作图、制表、方程拟合、计算特征量等手段，探索数据的结构合规律的一种数据分析方法、从而保住我们后期更好的进行特征工程和建立模型。

加载数据及基本的预处理，需要特别注意：

1. 数据分为两部分，训练数据和测试数据，并且读取方式是不一样的；

2. 训练数据非常大，不能直接展开，可以选择一部分数据进行展开，了解数据基本情况；

3. 与常规读取数据的方式不同，这里使用xarray进行数据加载；

   data_dir = 'weather'
   # 按时间进行数据块划分
   # 定义了一个`chunk_time` 函数（区分时间函数）
   def chunk_time(ds):
       # 接受一个xarray的数据集（`ds`）作为参数。
       #它使用`ds.dims`获取数据集的维度信息，并将其存储在`dims`字典中。
       dims = {k:v for k, v in ds.dims.items()}
       #将dims中的维度的长度修改为1，为什么要修改呢？
       dims['time'] = 1
       # 使用了xarray库中的chunk方法来对数据集ds进行分块操作。chunk方法需要传入一个参数dims，它表示要对哪些维度进行分块。
       ds = ds.chunk(dims)
       # ds.chunk会根据维度进行分块操作，之前把time设成了1所以每个时间维度上有一个元素。
       return ds
   
   # 加载训练数据
   def load_dataset():
       # 将ds定义为空数组
       ds = []
       # 从2007到2012年进行遍历，并赋值给y 
       for y in range(2007, 2012):
   
   # os.path.join()函数来创建一个路径字符串，里面传入数据目录路径变量，一开始定的weather，以及f'weather_round1_train_{y}'是一个表示文件名的格式化字符串，其中{y}是一个时间变量也就是从2007到2012。
   # data_name 是数据路径的字符串
   
           data_name = os.path.join(data_dir, f'weather_round1_train_{y}')
       
       # xarray库中的open_zarr函数打开Zarr数据集的代码。open_zarr函数用于读取Zarr格式的数据，并将其存储为xarray.Dataset对象。
   # 在这里，data_name是指Zarr数据集的文件路径或URL。consolidated=True参数表示如果数据集被分块存储，打开时将会加载所有块的元数据进行合并。
   # 通过执行这行代码，可以使用x变量来访问和操作Zarr数据集中的数据。
   
           x = xr.open_zarr(data_name, consolidated=True)
       
      # 这会将data_name打印为数据集的名称，并通过x.time.values[0]和x.time.values[-1]打印出时间范围的起始和结束时间。
   
           print(f'{data_name}, {x.time.values[0]} ~ {x.time.values[-1]}')
       
       # 通过执行ds.append(x)，可以将x数据集中的所有变量添加到ds数据集中，从而创建一个包含原始数据和新数据的完整数据集。
       
           ds.append(x)
           
          # xr.concat()函数是用于在指定维度上连接（拼接）多个xarray.Dataset对象的函数。
   #在代码中，xr.concat(ds, 'time')用于将 ds 数据集列表沿着 'time' 维度进行连接。
   #这个函数会将ds中的多个数据集按照时间维度进行拼接，生成一个新的xarray.Dataset对象。
   #注意，ds应该是一个包含多个xarray.Dataset对象的列表，而不是一个单个数据集。
   #通过执行ds = xr.concat(ds, 'time')，你可以将列表中的所有数据集按照时间维度进行连接，从而生成一个包含所有数据的新数据集，并将其赋值给变量ds
   
       ds = xr.concat(ds, 'time')
       # 将数据集按照时间分块
       ds = chunk_time(ds)
       return ds
   # 调用load_dataset()函数，加载训练数据   
   # .x表示从加载的数据集中提取一个名为x的变量
   ds = load_dataset().x
   
   num_step = 20 # 5天的数据，每6小时一个时刻，5*6/24=20，所以输出的时间维度一共有20个。5个地面变量，H和W都是161表示空间范围
   shape = ds.shape # 四维数据，time x channel x lat x lon
   
   times = ds.time.values # 提取对应的所有time值
   # times是一个可迭代对象（如列表或数组），而slice(1, -20)则是一个切片对象，用于选择times中的一部分元素。
   
   # 在这里，slice(1, -20)表示从索引1（第二个元素）开始，直到倒数第20个元素之前结束（不包括倒数第20个元素）。（为什么要这样切片呢？）这将创建一个切片对象，包含了times中的一部分元素。
   # 然后，init_times将通过这个切片对象来选择times中的相应元素，作为初始化时间（initial times）。
   
   init_times = times[slice(1, -num_step)] 
   
   num_data = len(init_times) # 计算init_times中元素的数量，并将结果保存在num_data变量中。
   
   names = list(ds.channel.values) # 提取所有变量名
   test_names = names[-5:] # 最后5个变量为目标变量，即我们需要预测的目标
   
   print(f'\n shape: {shape}')
   print('\n times: {} ~ {}'.format(times[0], times[-1]))
   print('\n init_times: {} ~ {}'.format(init_times[0], init_times[-1]))
   print(f'\n names: {names}')
   print(f'\n test_names: {test_names}\n')   

   • data_dir 是存储数据的目录，可以根据实际情况进行修改。

   • 定义一个chunk_time 函数（区分时间函数），接受一个xarray的数据集（ds）作为参数。

   • 它使用ds.dims获取数据集的维度信息，并将其存储在dims字典中。

   • ds.dims.items()返回的是一个迭代器，它遍历数据集的维度字典的每一个键值对，每次迭代返回一个键值对的元组。

   • 然后，使用ds.chunk方法根据dims字典对数据集进行划分，将时间维度划分为一个chunk。最后，返回划分后的数据集。

   • 分块操作是将数据集的维度划分为多个块，以提高计算性能和内存利用率。

将数据打印出来看看

print(ds.time.values)
print(ds.channel.values)
print(ds.lat.values)
print(ds.lon.values)

ds[0,0,:,:].values
# 对应维度
# ds[0,0,:,:].shape
# (161, 161)

如果想要看具体的变量值，则可以通过如下代码查看，因为数据比较大，所以仅看第一个time值（2007-01-01T00:00:00.000000000），对应的z50，lat和lon下的所有z50值

可以看到一个时刻对应某个变量的值应该是161*161个，即不同经纬度下的值。所以需要预测的结果是20个时刻对应5个目标变量下不同经纬度下的值，即20 x 5 x 161 x 161。

以上都是训练数据。

接下来看看测试数据，同样先加载下来看看：

file_lists = os.listdir(f'{path}/weather_round1_test/input')
for file in file_lists:
    input_data = torch.load(f'{path}/weather_round1_test/input/{file}')
    break
print(file_lists) # 打印所有输入数据，每个pt文件为一个测试集输入数据，pt文件直接没有任何相关性
# file_lists是一个包含特定路径下所有文件名的列表。使用os.listdir()函数，从{path}/weather_round1_test/input路径中获取所有文件的名称。

# 在接下来的循环中，对file_lists中的每个文件进行迭代。对于每个文件，将使用torch.load()函数加载文件，加载的路径为{path}/weather_round1_test/input/{file}。这样可以得到一个input_data对象，它包含了文件的内容。

# 循环中的break语句用于在第一个文件迭代后退出循环，因此只会加载第一个文件的数据。

# 最后，使用print(file_lists)打印出所有的文件列表，每个.pt文件都是一个测试数据集的输入文件，这些.pt文件之间没有任何相关性。

input_data.shape, input_data
# shape为（2,70,161,161)，即输入数据为两个时刻的数据

一般时间序列预测问题都有较长的历史数据，而本次赛题却只有历史两个时刻单位的数据，这就导致常见的时间序列特征提取方式在这个赛题中无法使用，或者极具困难。

3、数据清洗

数据和特征决定了机器学习的上限，而模型和算法只是逼近这个上限而已。俗话说：garbage in, garbage out。分析完数据后，特征工程前，必不可少的步骤是对数据清洗。

数据清洗的作用是利用有关技术如数理统计、数据挖掘或预定义的清理规则将脏数据转化为满足数据质量要求的数据。主要包括缺失值处理、异常值处理、数据分桶、特征归一化/标准化等流程。

数据非常干净，不需要进行数据清洗工作。

4、统计策略进行

import torch,os
path = 'weather' 
file_lists = os.listdir(f'{path}/weather_round1_test/input')
# 遍历读取每个测试集输入数据文件
for file in file_lists:
    input_data = torch.load(f'{path}/weather_round1_test/input/{file}')
    # 计算两个时刻对应目标变量的均值，并复制20份，即20个时刻
    output_data = torch.tile(torch.mean(input_data, axis=0)[-5:,:,:], (20,1,1,1))
    # output_data = torch.tile(torch.mean(input_data, axis=0)[-5:,:,:], (20,1,1,1))通过计算两个时刻对应目标变量的均值，并将其复制20份，可以得到一个结果。

# 这个计算过程是基于某种特定的统计策略。具体的策略可能与数据集的特点和使用的模型有关。

# 数据集的特点：目标变量是一个时间序列，其中每个时刻都有多个观测值（多个特征）。每个时刻的观测值是一个向量，其维度为 (T, C, H, W)，其中 T 是时间，C 是通道数，H 和 W 是高度和宽度。

# 统计策略：要获取一段时间内的整体趋势，可以计算这段时间内观测值的平均值。这样可以将多个时刻的观测值合并为一个均值，从而捕捉它们的整体特征。

# 在给定的代码中，通过 torch.mean(input_data, axis=0) 对 input_data 沿着第0维（通常是批次维度）进行平均操作，得到了一个表示整体趋势的张量。

# 然后，通过 [-5:,:,:] 选择最后5个时刻的观测值。这可能是基于一个假设，即最近的观测值对于预测下一个时刻的目标变量更具有重要性。

# 最后，通过 torch.tile() 函数将选择的最后5个时刻的观测值复制20次，得到形状为 (20, 5, :, :) 的张量。这样做可能是为了将这段时间内的整体趋势扩展到更长的时间范围，以便与模型的输入形状匹配。


    # 计算结果保存到output文件作为输出结果
    torch.save(output_data.half(), f'output/{file}', ) # 精度转为半精度浮点类型

# 现在就可以得到答案的压缩包啦~~~
_ = !zip -r output.zip output/

#如果上面代码无响应，运行下面这段代码
import zipfile,os
path='./output/' #压缩当前路径所有文件，输出zip文件
zipName = 'output.zip' #压缩后文件的位置及名称
f = zipfile.ZipFile( zipName, 'w', zipfile.ZIP_DEFLATED )
for dirpath, dirnames, filenames in os.walk(path):
    for filename in filenames:
        print(filename)
        f.write(os.path.join(dirpath,filename))
f.close()

Step2 搞懂机器学习BaseLine

我们将使用机器学习的方法搭建方案，模型选择catboost，在baseline部分已经进行了数据探索和数据清洗两个部分，进阶部分直接跳过。主要围绕数据预处理、特征工程、模型训练、模型验证进展开实践。

实践难点：

1. 数据格式与以往的不同：在数据挖掘比赛中，数据是四维的形式的呈现（time,channl,lat,lon），与以往直接转为dataframe格式完全不同，需要我们自行转换，并且自行构建训练样本和测试样本。
2. 可用的历史数据非常少：每份训练数据只提供历史2个时刻的输入文件，这也表明仅能使用2个时刻的数据提取特征，如何使用2个时刻数据刻画未来20个时刻的目标变量信息能为难点。
3. 数据量非常大：每份测试数据需要预测的结果为（20，5，161，161），对应205161*161条样本，有300份测试数据，那就是777630000条样本，就单单进行样本构建、特征转换和结果预测就非常耗时了。训练数据是远大于测试数据，不过可以提取部分数据输入到模型进行训练，并不需要全部使用。

1、导入模块

# 是一个数学基础库，可以计算矩阵，导入numpy库将其命名为np
import numpy as np

# 导入的是一个时间序列数据相关的库命名为pd
import pandas as pd

# Python 模块，支持带有维度、坐标和属性标注的多维数组
import xarray as xr
# 从collections模块导入了两个类：defaultdict和Counter。

# defaultdict是一个用于创建字典的类，它的特点是在访问不存在的键时不会引发KeyError异常，而是返回一个默认值。这个默认值可以在创建defaultdict对象时指定。使用defaultdict可以简化字典操作，特别是在处理不存在的键时非常有用。

# Counter是一个用于计数的类，它接受一个可迭代对象（如列表或字符串）作为输入，并返回一个字典，其中包含了每个元素出现的频率。可以使用Counter方便地进行单词统计、频率计数等操作。
from collections import defaultdict, Counter

# 从catboost模块导入了CatBoostRegressor类。CatBoost是一个梯度提升框架，专门用于处理分类和回归问题。CatBoostRegressor是CatBoost库中用于回归问题的模型类。

# CatBoost框架在梯度提升的基础上引入了一些优化策略和创新功能，以提高模型的性能和准确性。它具有处理大量类别特征的能力，并且可以自动处理缺失值。此外，CatBoost还支持并行训练，具有梯度剪切和自动特征缩放等功能。
from catboost import CatBoostRegressor

# 从sklearn.model_selection模块分别导入了三个类：StratifiedKFold、KFold和GroupKFold。这些类都是用于划分数据集的交叉验证策略。

# StratifiedKFold是用于分层抽样的交叉验证类。它保持了每个折叠中类别的分布与原始数据集中的类别分布相似。这对于不平衡数据集或需要保持样本类别比例的任务非常有用。

# KFold是标准的K折交叉验证类。它将数据集分成K个子集，其中每个子集都作为验证集一次，其余的K-1个子集组成训练集。这种方法适用于一般的机器学习任务。

# GroupKFold是一种针对分组数据的交叉验证类。它确保同一组的样本要么完全出现在训练集中，要么完全出现在验证集中。这对于处理具有相关性或依赖性的样本（例如时间序列数据或多个观察来自同一实体）非常有用。
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold

#从sklearn.metrics模块导入了mean_squared_error函数。mean_squared_error是一个用于计算均方误差（Mean Squared Error, MSE）的函数。

# 均方误差是回归问题中常用的评估指标之一，用于度量预测值与实际值之间的差异。它计算预测值与真实值之差的平方，并对所有样本进行求和后取平均。MSE越接近0，表示回归模型的预测值与实际值越接近。

# 使用mean_squared_error函数，您可以计算模型的预测值与真实值之间的均方误差。
from sklearn.metrics import mean_squared_error

# 从joblib模块导入了Parallel和delayed两个函数。这两个函数在并行计算和批量处理任务时非常有用。

# Parallel是一个用于并行计算的类，它可以方便地将函数应用于多个输入，并行地执行这些函数。通过指定n_jobs参数，您可以指定同时执行的并行任务数量。Parallel类可以显著提高计算密集型任务的执行效率。

# delayed是一个装饰器函数，用于将函数延迟执行。它用于在Parallel类的上下文中对函数进行修饰，以指定要并行执行的函数。通过使用delayed装饰器，可以将需要并行化的函数包装起来，以便在并行执行时被调用。

# 这两个函数通常一起使用，可以简化并行化任务的实现。可以创建一个Parallel对象，并将需要并行处理的函数使用delayed装饰器进行修饰，然后通过调用Parallel对象来并行执行这些函数。
from joblib import Parallel, delayed

# 从tqdm.auto模块导入了tqdm函数。tqdm是一个用于在循环中显示进度条的 Python 包。它可以在命令行界面中实时显示循环迭代的进度，并根据迭代的速度动态调整进度条的更新频率。

# 使用tqdm函数，可以将迭代对象作为参数传递给它，并在循环中使用它来实时显示进度条。根据迭代对象的大小和循环的进度，tqdm会自动更新并显示当前的进度条和预计剩余时间。这对于长时间运行的循环或处理大量数据的任务非常有用，可以提供可视化的进度反馈。
from tqdm.auto import tqdm
# 导入了一系列Python模块和库，它们具有以下功能：

# sys模块提供了对Python运行时环境变量和命令行参数的访问。
# os模块提供了与操作系统交互的函数，例如文件和目录操作。
# gc模块提供了对Python垃圾回收机制的控制，可以手动管理内存。
# argparse模块用于解析命令行参数和选项，可以帮助您编写具有可选参数和帮助信息的命令行工具。
# warnings模块用于控制警告消息的显示和处理。
# torch是一个用于深度学习的Python库，提供了张量操作、神经网络构建和训练等功能。
# glob模块用于在文件系统中查找符合指定模式的文件路径。
import sys,os,gc,argparse,warnings,torch,glob

# warnings.filterwarnings('ignore')是一个用于控制警告消息显示的函数调用。通过调用filterwarnings函数并传递参数'ignore'，您可以设置Python解释器忽略特定类型的警告消息。

# 警告消息是提示用户潜在问题或错误的信息，但并不会导致程序终止。有时，警告消息可能会在特定情况下频繁出现，但您已经确认这些警告对程序的正常执行没有影响，或者您没有兴趣在控制台中看到这些警告。

# 通过调用warnings.filterwarnings('ignore')，您可以将所有类型的警告都忽略，从而阻止它们显示在控制台上。这样，您就可以抑制不必要的警告消息，并使程序的输出更加干净。
warnings.filterwarnings('ignore')

2、加载数据（与前面基本一致）

path = 'weather'

def chunk_time(ds):
    dims = {k:v for k, v in ds.dims.items()}
    dims['time'] = 1
    ds = ds.chunk(dims)
    return ds

def load_dataset():
    ds = []
    for y in range(2007, 2012):
        data_name = os.path.join(path, f'weather_round1_train_{y}')
        x = xr.open_zarr(data_name, consolidated=True)
        print(f'{data_name}, {x.time.values[0]} ~ {x.time.values[-1]}')
        ds.append(x)
    ds = xr.concat(ds, 'time')
    ds = chunk_time(ds)
    return ds

ds = load_dataset().x

3、构建训练样本（将特征量进行维度转换和合并，从三维变成一维）

因为数据量比较大，所以仅提取10份训练数据，每份数据大小应该为(22,70,161,161)，共22个时刻的数据，前面2个时刻用来提取特征，后面20个时刻用来提取目标变量。另外需要特别注意的是：

1. 整个流程涉及很多数组的操作，如提取数据、维度转换、数组合并等，需要仔细操作；

2. 分别提取时间、经纬度、第一个时刻的70个变量、第二个时刻的70个变量、两个时刻70个变量的差值。需要特别注意，对于未来20个时刻，除了时间信息是不一样的外，其余特征变量完全一致，主要刻画相同特征不同时间的情况下，目标变量的变化；

   # 将ds.lat和ds.lon的值分别转换为列表并赋值给变量lats和lons。
   lats = ds.lat.values.tolist()
   lons = ds.lon.values.tolist()
   
   # 对齐测试数据，每份数据应该是22个时刻
   # 获取多份训练数据，滑动10次，每次24个时刻（6天），每次滑动提取最后22个时刻作为一份训练数据
   
   #首先训练数据是空数组
   train_data = []
   for i in range(10):
       if i == 0:
           # 选用最后22个索引，
           print(-22, 0)
           train_data.append(ds[-22:, :, :, :].values)
       else:
           # 选择倒数22个时刻的数据作为训练数据
           idx = i * 24
           train_data.append(ds[-22-idx:-idx, :, :, :].values)
           print(-22-idx,-idx)
           
   # 经纬度标识字段
   # 每个循环迭代，i 和 j 的范围是从0到160，这将生成 161*161 个坐标对。
   
   # 最后，通过使用 numpy 的 array 函数，将 latlon_vals 转换为一个 numpy 数组，维度为 (161*161, 2)。这意味着该数组的行数是 161*161，每行有两个元素代表一个经纬度坐标。
   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) #（161*161，2）
   
   # 确定好20份训练数据后，接下来需要提取特征
   for flag, dat in tqdm(enumerate(train_data)):
       # 第一时刻特征
       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)
       
       # 在这段代码中，使用了 tqdm 函数来迭代 train_data 列表，并给每个迭代元素增加一个索引标志 flag。
   
   # 在每个迭代中，对数据 dat 进行了一系列操作，如下所示：
   
   # 第一时刻特征（first_feat）的计算：通过使用索引 [0,:,:,:]，从 dat 中选择第一个时间步的数据，然后通过重塑和转置操作将形状从 (70, 161, 161) 变为 (161*161, 70)。这样做，以便将原始的三维特征数据转换为一个二维特征矩阵。
   
   # 第二时刻特征（second_feat）的计算：通过使用索引 [1,:,:,:]，从 dat 中选择第二个时间步的数据，然后进行与上述相似的重塑和转置操作，将形状变为 (161*161, 70)。
   
   # 两个时刻差分特征（diff_feat）的计算：通过将第二个时间步的数据减去第一个时间步的数据，然后进行与上述相似的重塑和转置操作，将形状变为 (161*161, 70)。
   # 通过将这个三维矩阵重塑为 (161*161, 70) 的形状，可以将整个矩阵展平为一个长向量，其中每个元素就代表了矩阵中的一个特征值。这样，整个二维矩阵可以被表示为一个特征向量的集合。
   # 这些操作的目的是将原始的三维特征数据转换为二维特征矩阵，以便在后续的分析或处理中使用。
   
       # 构建训练样本
       tmp_dat = dat[2:,-5:,:,:] # (20, 5, 161, 161)
       # 通过切片 [2:,-5:,:,:]，你选取了 dat 数组的第3到最后一个时间步、倒数第5到最后一个通道、以及所有行和列。这样，你得到了一个新的数组 tmp_dat，其形状为 (20, 5, 161, 161)。
       
      # 使用了一个循环，迭代了20个时间步。
   
      # 在每个时间步中，进行了以下操作：
   
       for i in range(20): # 20个时刻
           # 时间特征、经纬特征
           # 时间特征的计算：通过使用 np.array([i]*161*161).reshape(161*161, 1)，创建了一个形状为 (161*161, 1) 的数组 time_vals，其中每个元素都是当前时间步 i。
           time_vals = np.array([i]*161*161).reshape(161*161,1) #（161*161，1）
           
           # 经纬特征的计算：将之前定义的 latlon_vals 数组与时间特征 time_vals 拼接起来，使用 np.concatenate 函数，按照列的方向（axis=1）进行拼接，得到形状为 (161*161, 3) 的数组 sub_vals。
           sub_vals = np.concatenate((time_vals, latlon_vals), axis=1) #（161*161，3）
           
           # 添加历史特征，第一时刻特征、第二时刻特征、两个时刻差分特征
           # 添加历史特征：将之前计算的 第一时刻特征 first_feat、第二时刻特征 second_feat 和两个时刻差分特征 diff_feat，依次与 sub_vals 进行拼接。通过多次使用 np.concatenate 函数，将这些特征按照列的方向拼接到 sub_vals 的右侧。最终，sub_vals 的形状将变为 (161*161, 213)，其中包含时间特征、经纬特征和历史特征。
           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)
               sub_vals = np.concatenate((sub_vals, var_vals), axis=1) 
           # 将当前时间步的特征数据 sub_vals 添加到 all_vals 中。如果 i 是第一个时间步（i == 0），则直接将 sub_vals 赋值给 all_vals；否则，通过使用 np.concatenate 函数将 sub_vals 拼接到 all_vals 的下方，按行的方向进行拼接。
           if i == 0 :
               all_vals = sub_vals 
           else:
               all_vals = np.concatenate((all_vals, sub_vals), axis=0)
       # 最后，通过外层循环的迭代，将每个时间步的特征数据添加到 final_vals 中。如果 flag 是0，则直接将 all_vals 赋值给 final_vals；否则，通过使用 np.concatenate 函数将 all_vals 拼接到 final_vals 的下方，按行的方向进行拼接。        
       if flag == 0:
           final_vals = all_vals
       else:
           final_vals = np.concatenate((final_vals, all_vals), axis=0)
           # 最终，final_vals 将包含所有时间步的特征数据，形状为 (161*161*20, 218)。

将二维数组转为dataframe格式，并设置新的列名，然后选择入模特征，代码如下：

# 首先使用 pd.DataFrame 函数将 final_vals 数组转换为一个 Pandas DataFrame 对象，命名为 train_df。
train_df = pd.DataFrame(final_vals)
# 通过为 train_df 的列指定名称，对其列进行命名。列的名称包括：‘time’、‘lat’、‘lon’，以及以’feat_‘开头的210列特征名称，最后加上’t2m’、‘u10’、‘v10’、'msl’和’tp’这5个目标变量的名称。
train_df.columns = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)] + ['t2m','u10','v10','msl','tp']

cols = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)]
# 将 final_vals 数组转换为了一个带有适当的列名称的 Pandas DataFrame 对象，并定义了一个存储选取列名称的列表 cols

4、训练模型

​ 模型选择使用catboost模型，并且将模型训练和模型推理单独构建，主要是方便对300份测试数据进行单独推理。另外需要注意，模型学习率learning_rate设置的比较大。

def train_model(train_x, train_y, label, seed=2023):
    # train_x：训练数据的特征集合，是一个 Pandas DataFrame。
    # 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)))
         # 设置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)
#这段代码定义了一个 train_model 函数，用于训练一个 CatBoostRegressor 模型。

5、测试集预测

在训练模型阶段提前将训练好的模型进行保存，在模型推理函数部分直接加载对应目标变量的模型即可，代码如下：

def inter_model(test_x, label, seed=2023):
    # test_x：测试数据的特征集合，是一个 Pandas DataFrame。
	# label：模型训练的目标变量名称。
	# seed：随机种子，用来控制训练过程的随机性，默认值为 2023。
    
    # n_splits：表示将数据划分为几个折。在这个例子中，将数据划分为 5 个折。
	# shuffle：表示是否在进行划分前对数据进行洗牌。设置为 True 表示进行洗牌，即在划分前对数据随机打乱顺序。
	# random_state：随机种子，用于控制洗牌时的随机性。设置为 seed，该值作为参数传入。
    kf = KFold(n_splits=5, shuffle=True, random_state=seed)
    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'}

    for i in range(0, 5):
        # print('************************************ {} {}************************************'.format(str(i+1), str(seed)))
        model = CatBoostRegressor(**params) 
        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个目标变量外。代码如下：

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）
        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）
        
        if i == 0 :
            all_vals = sub_vals
        else:
            all_vals = np.concatenate((all_vals, sub_vals), axis=0)
            
    df = pd.DataFrame(all_vals)
    df.columns = ['time','lat','lon'] + [f'feat_{i}' for i in range(210)]
    return df

结果生成函数，代码流程：

1. 生成测试样本，并构建特征；
2. 调用推理模型函数，分别预测5个目标变量的结果；
3. 将结果转为可提交格式(20,5,161,161)；
4. Array转为Tensor；
5. 保存为pt文件，精度为半精度浮点类型。

代码如下：

def get_parallel_result(file):
    input_data = torch.load(file)
    
    test_df = get_test_data(np.array(input_data)) # 生成测试样本，并构建特征
    
    res = test_df[['time']] # 保存结果
    for label in ['t2m','u10','v10','msl','tp']:
        cat_pred = inter_model(test_df[cols], label) # 模型推理，预测测试集结果
        res[label] = cat_pred  # 将预测结果存储在字典中
        
    # 转为提交格式
    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)
            
            if i == 0:
                all_vals = sub_vals
            else:
                all_vals = np.concatenate((all_vals, sub_vals), axis=0)
        
        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)
    
    final_vals = torch.tensor(final_vals)
    save_file = file.split('/')[-1]
    torch.save(final_vals.half(), f'output/{save_file}', ) # 精度转为半精度浮点类型

为了加快运行速度，这里使用Joblib中的Parallel和delayed实现并行化处理，执行get_parallel_result函数，代码如下：

注：n_jobs为开启的进程数，即是用n_jobs个CPU来计算。

Parallel(n_jobs=32)(
    delayed(get_parallel_result)(file_path)
      for file_path in tqdm(glob.glob(f'{path}/weather_round1_test/input/*.pt'))
)

6、结果输出（统一的没啥区别几乎）

# 现在就可以得到答案的压缩包啦~~~
_ = !zip -r output.zip output/
#如果上面代码无响应，运行下面这段代码
import zipfile,os
path='./output/' #压缩当前路径所有文件，输出zip文件
zipName = 'output.zip' #压缩后文件的位置及名称
f = zipfile.ZipFile( zipName, 'w', zipfile.ZIP_DEFLATED )
for dirpath, dirnames, filenames in os.walk(path):
    for filename in filenames:
        print(filename)
        f.write(os.path.join(dirpath,filename))
f.close()

优化1、扩大训练数据为15份（后续想想怎么能更多）

优化 2、学习率和迭代次数105（测试中）

优化3、特征优化

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 = np.var([dat[1,:,:,:], dat[0,:,:,:]], axis=0).reshape(70, 161*161).transpose() # (70, 161, 161) --> (161*161, 70)

优化4、评价指标

回归问题与梯度下降

梯度下降：

首先梯度类似于导数，

首先随便找一个点X0，然后算出Y0,然后计算机算出来Y0’,然后给计算机一个迭代的过程。

这个过程大概就是梯度下降的过程。

详细回归问题与梯度下降视频：

https://www.bilibili.com/video/BV1Wd4y1X7Yi/?spm_id_from=333.337.search-card.all.click&vd_source=cf0b651a9f3e6c7822ca8f914aeb9a28

那么回归又是什么意思？

由于我们现实中的数据存在一些误差，我们允许一些误差的存在，然后得到一组整体更好的参数，回归