2021 “AI Earth”人工智能创新挑战赛-AI助力精准气象和海洋预测BaseLine学习笔记

最新推荐文章于 2024-04-09 11:55:47 发布

JeffDingAI

最新推荐文章于 2024-04-09 11:55:47 发布

阅读量1.6k

点赞数 3

分类专栏： Datawhale学习笔记文章标签： python 机器学习

本文链接：https://blog.csdn.net/yichao_ding/article/details/114096881

版权

Datawhale学习笔记专栏收录该内容

80 篇文章 3 订阅

订阅专栏

一、比赛主页

链接：https://tianchi.aliyun.com/competition/entrance/531871/introduction

二、比赛介绍

本次赛题是一个时间序列预测问题。基于历史气候观测和模式模拟数据，利用T时刻过去12个月(包含T时刻)的时空序列（气象因子），构建预测ENSO的深度学习模型，预测未来1-24个月的Nino3.4指数，如下图所示：

三、背景数据描述

1. 数据简介

本次比赛使用的数据包括CMIP5/6模式的历史模拟数据和美国SODA模式重建的近100多年历史观测同化数据。每个样本包含以下气象及时空变量：海表温度异常(SST)，热含量异常(T300)，纬向风异常（Ua），经向风异常（Va），数据维度为（year,month,lat,lon）。对于训练数据提供对应月份的Nino3.4 index标签数据。

2. 训练数据标签说明

标签数据为Nino3.4 SST异常指数，数据维度为（year,month）。

CMIP(SODA)_train.nc对应的标签数据当前时刻Nino3.4 SST异常指数的三个月滑动平均值，因此数据维度与维度介绍同训练数据一致。

注：三个月滑动平均值为当前月与未来两个月的平均值。

3. 测试数据说明

测试用的初始场（输入）数据为国际多个海洋资料同化结果提供的随机抽取的n段12个时间序列，数据格式采用NPY格式保存，维度为（12，lat，lon, 4）,12为t时刻及过去11个时刻，4为预测因子，并按照SST,T300,Ua,Va的顺序存放。

测试集文件序列的命名规则：test_编号_起始月份_终止月份.npy，如test_00001_01_12_.npy。

四、BaseLine

（一）数据读取

1.引入包

安装工具
# !pip install netCDF4 
''' 
import pandas as pd
import numpy  as np
import tensorflow as tf
from tensorflow.keras.optimizers import Adam
import matplotlib.pyplot as plt
import scipy 
from netCDF4 import Dataset
import netCDF4 as nc
import gc
%matplotlib inline

2.SODA_label处理

标签含义：标签数据为Nino3.4 SST异常指数，数据维度为（year,month）。
CMIP(SODA)_train.nc对应的标签数据当前时刻Nino3.4 SST异常指数的三个月滑动平均值，因此数据维度与维度介绍同训练数据一致
注：三个月滑动平均值为当前月与未来两个月的平均值。

#将标签转化为我们熟悉的pandas形式
label_path       = './data/SODA_label.nc'
label_trans_path = './data/' 
nc_label         = Dataset(label_path,'r')
 
years            = np.array(nc_label['year'][:])
months           = np.array(nc_label['month'][:])

year_month_index = []
vs               = []
for i,year in enumerate(years):
    for j,month in enumerate(months):
        year_month_index.append('year_{}_month_{}'.format(year,month))
        vs.append(np.array(nc_label['nino'][i,j]))

df_SODA_label               = pd.DataFrame({'year_month':year_month_index}) 
df_SODA_label['year_month'] = year_month_index
df_SODA_label['label']      = vs

df_SODA_label.to_csv(label_trans_path + 'df_SODA_label.csv',index = None)
df_label.head()

3.SODA_train处理

SODA_train.nc中[0,0:36,:,:]为第1-第3年逐月的历史观测数据；

SODA_train.nc中[1,0:36,:,:]为第2-第4年逐月的历史观测数据；
…,
SODA_train.nc中[99,0:36,:,:]为第100-102年逐月的历史观测数据。
SODA_path        = './data/SODA_train.nc'
nc_SODA          = Dataset(SODA_path,'r')

4.自定义抽取对应数据&转化为df的形式；

def trans_df(df, vals, lats, lons, years, months):
    '''
        (100, 36, 24, 72) -- year, month,lat,lon 
    ''' 
    for j,lat_ in enumerate(lats):
        for i,lon_ in enumerate(lons):
            c = 'lat_lon_{}_{}'.format(int(lat_),int(lon_))  
            v = []
            for y in range(len(years)):
                for m in range(len(months)): 
                    v.append(vals[y,m,j,i])
            df[c] = v
    return df
year_month_index = []

years              = np.array(nc_SODA['year'][:])
months             = np.array(nc_SODA['month'][:])
lats             = np.array(nc_SODA['lat'][:])
lons             = np.array(nc_SODA['lon'][:])


for year in years:
    for month in months:
        year_month_index.append('year_{}_month_{}'.format(year,month))

df_sst  = pd.DataFrame({'year_month':year_month_index}) 
df_t300 = pd.DataFrame({'year_month':year_month_index}) 
df_ua   = pd.DataFrame({'year_month':year_month_index}) 
df_va   = pd.DataFrame({'year_month':year_month_index})
%%time
df_sst = trans_df(df = df_sst, vals = np.array(nc_SODA['sst'][:]), lats = lats, lons = lons, years = years, months = months)
df_t300 = trans_df(df = df_t300, vals = np.array(nc_SODA['t300'][:]), lats = lats, lons = lons, years = years, months = months)
df_ua   = trans_df(df = df_ua, vals = np.array(nc_SODA['ua'][:]), lats = lats, lons = lons, years = years, months = months)
df_va   = trans_df(df = df_va, vals = np.array(nc_SODA['va'][:]), lats = lats, lons = lons, years = years, months = months)
label_trans_path = './data/'
df_sst.to_csv(label_trans_path  + 'df_sst_SODA.csv',index = None)
df_t300.to_csv(label_trans_path + 'df_t300_SODA.csv',index = None)
df_ua.to_csv(label_trans_path   + 'df_ua_SODA.csv',index = None)
df_va.to_csv(label_trans_path   + 'df_va_SODA.csv',index = None)

5.CMIP_label处理

label_path       = './data/CMIP_label.nc'
label_trans_path = './data/'
nc_label         = Dataset(label_path,'r')
 
years            = np.array(nc_label['year'][:])
months           = np.array(nc_label['month'][:])

year_month_index = []
vs               = []
for i,year in enumerate(years):
    for j,month in enumerate(months):
        year_month_index.append('year_{}_month_{}'.format(year,month))
        vs.append(np.array(nc_label['nino'][i,j]))

df_CMIP_label               = pd.DataFrame({'year_month':year_month_index}) 
df_CMIP_label['year_month'] = year_month_index
df_CMIP_label['label']      = vs

df_CMIP_label.to_csv(label_trans_path + 'df_CMIP_label.csv',index = None)

6.CMIP_train处理

CMIP_train.nc中[0,0:36,:,:]为CMIP6第一个模式提供的第1-第3年逐月的历史模拟数据；
…,
CMIP_train.nc中[150,0:36,:,:]为CMIP6第一个模式提供的第151-第153年逐月的历史模拟数据；

CMIP_train.nc中[151,0:36,:,:]为CMIP6第二个模式提供的第1-第3年逐月的历史模拟数据；
…,
CMIP_train.nc中[2265,0:36,:,:]为CMIP5第一个模式提供的第1-第3年逐月的历史模拟数据；
…,
CMIP_train.nc中[2405,0:36,:,:]为CMIP5第二个模式提供的第1-第3年逐月的历史模拟数据；
…,
CMIP_train.nc中[4644,0:36,:,:]为CMIP5第17个模式提供的第140-第142年逐月的历史模拟数据。

其中每个样本第三、第四维度分别代表经纬度（南纬55度北纬60度，东经0360度），所有数据的经纬度范围相同。
CMIP_path       = './data/CMIP_train.nc'
CMIP_trans_path = './data'
nc_CMIP  = Dataset(CMIP_path,'r')
nc_CMIP.variables.keys()

# dict_keys(['sst', 't300', 'ua', 'va', 'year', 'month', 'lat', 'lon'])
nc_CMIP['t300'][:].shape

# (4645, 36, 24, 72)
year_month_index = []

years              = np.array(nc_CMIP['year'][:])
months             = np.array(nc_CMIP['month'][:])
lats               = np.array(nc_CMIP['lat'][:])
lons               = np.array(nc_CMIP['lon'][:])

last_thre_years = 1000
for year in years:
    '''
        因为内存限制,我们暂时取最后1000个year的数据，如果内存够强大可以注释掉if
    '''
    if year >= 4645 - last_thre_years:
        for month in months:
            year_month_index.append('year_{}_month_{}'.format(year,month))

df_CMIP_sst  = pd.DataFrame({'year_month':year_month_index}) 
df_CMIP_t300 = pd.DataFrame({'year_month':year_month_index}) 
df_CMIP_ua   = pd.DataFrame({'year_month':year_month_index}) 
df_CMIP_va   = pd.DataFrame({'year_month':year_month_index})

def trans_thre_df(df, vals, lats, lons, years, months, last_thre_years = 1000):
    '''
        (4645, 36, 24, 72) -- year, month,lat,lon 
    ''' 
    for j,lat_ in (enumerate(lats)):
#         print(j)
        for i,lon_ in enumerate(lons):
            c = 'lat_lon_{}_{}'.format(int(lat_),int(lon_))  
            v = []
            for y_,y in enumerate(years):
                '''
                    因为内存限制,我们暂时取最后1000个year的数据，如果内存够强大可以注释掉if
                '''
                if y >= 4645 - last_thre_years:
                    for m_,m in  enumerate(months): 
                        v.append(vals[y_,m_,j,i])
            df[c] = v
    return df
%%time
df_CMIP_sst  = trans_thre_df(df = df_CMIP_sst,  vals   = np.array(nc_CMIP['sst'][:]),  lats = lats, lons = lons, years = years, months = months)
df_CMIP_sst.to_csv(CMIP_trans_path + 'df_CMIP_sst.csv',index = None)
del df_CMIP_sst
gc.collect()

df_CMIP_t300 = trans_thre_df(df = df_CMIP_t300, vals   = np.array(nc_CMIP['t300'][:]), lats = lats, lons = lons, years = years, months = months)
df_CMIP_t300.to_csv(CMIP_trans_path + 'df_CMIP_t300.csv',index = None)
del df_CMIP_t300
gc.collect()

df_CMIP_ua   = trans_thre_df(df = df_CMIP_ua,   vals   = np.array(nc_CMIP['ua'][:]),   lats = lats, lons = lons, years = years, months = months)
df_CMIP_ua.to_csv(CMIP_trans_path + 'df_CMIP_ua.csv',index = None)
del df_CMIP_ua
gc.collect()

df_CMIP_va   = trans_thre_df(df = df_CMIP_va,   vals   = np.array(nc_CMIP['va'][:]),   lats = lats, lons = lons, years = years, months = months)
df_CMIP_va.to_csv(CMIP_trans_path + 'df_CMIP_va.csv',index = None)
del df_CMIP_va
gc.collect()

（二）数据建模

1.引入包

import pandas as pd
import numpy  as np
import tensorflow as tf
from tensorflow.keras.optimizers import Adam
import matplotlib.pyplot as plt
import scipy 
import joblib
from netCDF4 import Dataset
import netCDF4 as nc 
from tensorflow.keras.callbacks import LearningRateScheduler, Callback
import tensorflow.keras.backend as K
from tensorflow.keras.layers import *
from tensorflow.keras.models import *
from tensorflow.keras.optimizers import *
from tensorflow.keras.callbacks import *
from tensorflow.keras.layers import Input 
import gc
%matplotlib inline

2.SODA_label处理

label_path       = './data/SODA_label.nc' 
nc_label         = Dataset(label_path,'r')
tr_nc_labels     = nc_label['nino'][:]

3.原始特征数据读取

SODA_path        = './data/SODA_train.nc'
nc_SODA          = Dataset(SODA_path,'r') 

nc_sst           = np.array(nc_SODA['sst'][:])
nc_t300          = np.array(nc_SODA['t300'][:])
nc_ua            = np.array(nc_SODA['ua'][:])
nc_va            = np.array(nc_SODA['va'][:])

4.神经网络框架

def RMSE(y_true, y_pred):
    return tf.sqrt(tf.reduce_mean(tf.square(y_true - y_pred)))

def RMSE_fn(y_true, y_pred):
    return np.sqrt(np.mean(np.power(np.array(y_true, float).reshape(-1, 1) - np.array(y_pred, float).reshape(-1, 1), 2)))

def build_model():  
    inp    = Input(shape=(12,24,72,4))  
    
    x_4    = Dense(1, activation='relu')(inp)   
    x_3    = Dense(1, activation='relu')(tf.reshape(x_4,[-1,12,24,72]))
    x_2    = Dense(1, activation='relu')(tf.reshape(x_3,[-1,12,24]))
    x_1    = Dense(1, activation='relu')(tf.reshape(x_2,[-1,12]))
     
    x = Dense(64, activation='relu')(x_1)  
    x = Dropout(0.25)(x) 
    x = Dense(32, activation='relu')(x)   
    x = Dropout(0.25)(x)  
    output = Dense(24, activation='linear')(x)   
    model  = Model(inputs=inp, outputs=output)

    adam = tf.optimizers.Adam(lr=1e-3,beta_1=0.99,beta_2 = 0.99) 
    model.compile(optimizer=adam, loss=RMSE)

    return model

5.训练集验证集划分

### 训练特征，保证和训练集一致
tr_features = np.concatenate([nc_sst[:,:12,:,:].reshape(-1,12,24,72,1),nc_t300[:,:12,:,:].reshape(-1,12,24,72,1),\
                              nc_ua[:,:12,:,:].reshape(-1,12,24,72,1),nc_va[:,:12,:,:].reshape(-1,12,24,72,1)],axis=-1)

### 训练标签，取后24个
tr_labels = tr_nc_labels[:,12:] 

### 训练集验证集划分
tr_len     = int(tr_features.shape[0] * 0.8)
tr_fea     = tr_features[:tr_len,:].copy()
tr_label   = tr_labels[:tr_len,:].copy()
val_len     = tr_features.shape[0]-tr_len
val_fea     = tr_features[tr_len:,:].copy()
val_label   = tr_labels[tr_len:,:].copy()

6.模型训练

#### 构建模型
model_mlp     = build_model()
#### 模型存储的位置
model_weights = './model_baseline/model_mlp_baseline.h5'

checkpoint = ModelCheckpoint(model_weights, monitor='val_loss', verbose=0, save_best_only=True, mode='min',
                             save_weights_only=True)

plateau        = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=5, verbose=1, min_delta=1e-4, mode='min')
early_stopping = EarlyStopping(monitor="val_loss", patience=20)
history        = model_mlp.fit(tr_fea, tr_label,
                    validation_data=(val_fea, val_label),
                    batch_size=4096, epochs=200,
                    callbacks=[plateau, checkpoint, early_stopping],
                    verbose=2)

7.模型预测

prediction = model_mlp.predict(val_fea)

8. Metrics

from   sklearn.metrics import mean_squared_error
def rmse(y_true, y_preds):
    return np.sqrt(mean_squared_error(y_pred = y_preds, y_true = y_true))

def score(y_true, y_preds):
    accskill_score = 0
    rmse_scores    = 0
    a = [1.5] * 4 + [2] * 7 + [3] * 7 + [4] * 6
    y_true_mean = np.mean(y_true,axis=0) 
    y_pred_mean = np.mean(y_preds,axis=0) 
#     print(y_true_mean.shape, y_pred_mean.shape)

    for i in range(24): 
        fenzi = np.sum((y_true[:,i] -  y_true_mean[i]) *(y_preds[:,i] -  y_pred_mean[i]) ) 
        fenmu = np.sqrt(np.sum((y_true[:,i] -  y_true_mean[i])**2) * np.sum((y_preds[:,i] -  y_pred_mean[i])**2) ) 
        cor_i = fenzi / fenmu
    
        accskill_score += a[i] * np.log(i+1) * cor_i
        rmse_score   = rmse(y_true[:,i], y_preds[:,i])
#         print(cor_i,  2 / 3.0 * a[i] * np.log(i+1) * cor_i - rmse_score)
        rmse_scores += rmse_score 
    
    return  2 / 3.0 * accskill_score - rmse_scores
print('score', score(y_true = val_label, y_preds = prediction))

六、提交BaseLine

参见：https://blog.csdn.net/yichao_ding/article/details/113856633

七、BaseLine结果

JeffDingAI

关注

3
点赞
踩
14

收藏

觉得还不错? 一键收藏
1
评论
2021 “AI Earth”人工智能创新挑战赛-AI助力精准气象和海洋预测BaseLine学习笔记

一、比赛主页链接：https://tianchi.aliyun.com/competition/entrance/531871/introduction二、比赛介绍本次赛题是一个时间序列预测问题。基于历史气候观测和模式模拟数据，利用T时刻过去12个月(包含T时刻)的时空序列（气象因子），构建预测ENSO的深度学习模型，预测未来1-24个月的Nino3.4指数，如下图所示：三、背景数据描述1. 数据简介本次比赛使用的数据包括CMIP5/6模式的历史模拟数据和美国SODA模式重建的近
复制链接

扫一扫