氣象海洋預測-Task2 數據分析

數據分析是我們解決一個數據挖掘任務的重要一環，通過數據分析，我們可以了解標簽的分布、數據中存在的缺失值和異常值、特征與標簽之間的相關性、特征之間的相關性等，并根據數據分析的結果，指導我們后續的特征工程以及模型的選擇和設計。

在本次任務中，我們將探索賽題中給出的兩份訓練數據，可視化分析四個氣象特征的分布情況，思考如何進行特征工程以及如何選擇或設計模型來實現我們的預測任務。

學習目標

學習如何探索并可視化分析氣象數據。

根據數據分析結果思考以下兩個問題：

• 能否構造新的特征？
• 選擇或設計什么樣的模型進行預測？

內容介紹

以SODA數據集為例，探索標簽的分布情況

• 可視化總體數據中Nino3.4指數隨時間的變化曲線
• 可視化樣本中Nino3.4指數隨時間的變化曲線

以SST特征為例，進行海陸掩膜和插值分析

• 進行陸地掩膜，并對海洋部分的SST進行插值
• 以樣本0為例，可視化每個月海洋的SST分布情況

以CMIP數據集為例，進行缺失值分析

• 統計四個氣象特征中缺失值的數量
• 以Ua特征為例，可視化分析缺失值的情況

溫度場和風場可視化

• 以樣本0的第0月為例，可視化溫度場和風場

總結

代碼示例

在本次任務中我們需要用到global_land_mask這個Python庫，這個庫可以根據輸入的緯度和經度判斷該點是否在陸地上。

安裝global_land_mask

!pip install global_land_mask Collecting global_land_maskDownloading global_land_mask-1.0.0-py3-none-any.whl (1.8 MB) [K |████████████████████████████████| 1.8 MB 601 kB/s eta 0:00:01 [?25hInstalling collected packages: global-land-mask Successfully installed global-land-mask-1.0.0 [33mWARNING: Running pip as root will break packages and permissions. You should install packages reliably by using venv: https://pip.pypa.io/warnings/venv[0m import numpy as np import matplotlib.pyplot as plt import random import netCDF4 import seaborn as sns from global_land_mask import globe from scipy import interpolate plt.rcParams['font.sans-serif'] = ['SimHei'] #中文支持 %matplotlib inline

以SODA數據集為例，探索標簽的分布情況

# 讀取SODA數據# 存放數據的路徑 path = '/kaggle/input/ninoprediction/' data = netCDF4.Dataset(path + 'SODA_train.nc') label = netCDF4.Dataset(path + 'SODA_label.nc') label = np.array(label.variables['nino']) print(data.variables['sst'].shape) (100, 36, 24, 72) label.shape (100, 36)

可以看到，數據集中的每個樣本是以某一年為起始的接下來36個月的氣象數據，同樣的，標簽也是以這一年為起始的接下來36個月中每個月的Nino3.4指數。然而一年只有12個月，怎么會有36個月的數據呢？

我們不妨來查看一下，將每個樣本中的12個月進行拼接時，Nino3.4指數的變化曲線。

# 分別將樣本中的0-12月、12-24月、24-36月進行拼接，繪制Nino3.4指數的變化曲線 plt.figure(figsize=(8, 10)) plt.subplots_adjust(hspace=0.4, wspace=0.4) for i in range(3):label_all = [label[j, 12*i:12*(i+1)] for j in range(label.shape[0])]label_all = np.concatenate(label_all, axis=0)plt.subplot(3, 1, i+1)plt.plot(label_all, 'k', linewidth=1)plt.xlabel('Time / month')plt.ylabel('Nino')plt.title('nino: month{}-{}'.format(12*i, 12*(i+1))) plt.show()

從圖中可以看出，三條曲線的變化趨勢是完全相同的，只是在時間上有12個月的位移。這說明，重疊部分的標簽是相同的，給出這樣的樣本的目的是可以以前12個月作為模型的輸入X、后24個月為預測目標Y構建訓練樣本。

進一步思考，將每個樣本構建一個訓練樣本，那么我們所能得到的訓練數據量就只有4645（CMIP數據）+100（SODA數據）=4745條，這樣小的數據量用于模型訓練必然是不夠的，那么如何增加數據量呢？這個問題留待大家思考。

同樣的，我們可以隨機抽取三個樣本，查看每個樣本中的Nino3.4指數變化曲線。

# 隨機抽取三個樣本，繪制Nino3.4指數變化曲線 plt.figure(figsize=(8, 10)) plt.subplots_adjust(hspace=0.4, wspace=0.4) for i in range(3):y = random.randint(0, label.shape[0])plt.subplot(3, 1, i+1)plt.plot(label[y], 'k', linewidth=1)plt.xlabel('Time / month')plt.ylabel('Nino')plt.title('nino: year{}'.format(y)) plt.show()

以SST特征為例，進行海陸掩膜和插值分析

在給定數據中，經度和緯度坐標都是離散的，每隔5度有一個坐標點，在這樣的經緯度坐標下的SST值也是離散的，因此我們以樣本0第0月的SST數據為例，用插值函數來擬合經緯度坐標與SST值之間的函數關系，得到平滑的SST分布。

lon = np.array(data.variables['lon']) lat = np.array(data.variables['lat']) x = lon y = lat # 以緯度和經度生成網格點坐標矩陣 xx, yy = np.meshgrid(x, y) # 取樣本0第0月的SST值 z = data.variables['sst'][0, 0] # 采用三次多項式插值，得到z = f(x, y)的函數f f = interpolate.interp2d(x, y, z, kind = 'cubic') x array([ 0., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50.,55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 105.,110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,165., 170., 175., 180., 185., 190., 195., 200., 205., 210., 215.,220., 225., 230., 235., 240., 245., 250., 255., 260., 265., 270.,275., 280., 285., 290., 295., 300., 305., 310., 315., 320., 325.,330., 335., 340., 345., 350., 355.])

經度的實際取值是從-180°到0°到180°，而給定的數據中經度的取值是0°到360度每間隔5°取一個坐標值，因此在之后判斷是否為陸地時需要轉換為實際的經度。

y array([-55., -50., -45., -40., -35., -30., -25., -20., -15., -10., -5.,0., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50.,55., 60.])

數據中緯度的取值也是每間隔5°取一個坐標值。

判斷每個經緯度坐標點是否在陸地上，用空值遮蓋陸地部分，便于觀察海洋上SST的分布。

# 將經度x轉換為實際經度重新生成網格點坐標矩陣 lon_grid, lat_grid = np.meshgrid(x-180, y) # 判斷坐標矩陣上的網格點是否為陸地 is_on_land = globe.is_land(lat_grid, lon_grid) is_on_land = np.concatenate([is_on_land[:, x >= 180], is_on_land[:, x < 180]], axis=1) # 進行陸地掩膜，將陸地的SST設為空值 z[is_on_land] = np.nan # 可視化掩膜結果，黃色為陸地，紫色為海洋 plt.imshow(is_on_land[::-1, :]) <matplotlib.image.AxesImage at 0x7f0044d0edd0>

z.shape (24, 72) # 可視化海洋上的SST分布 plt.imshow(z[::-1, :], cmap=plt.cm.RdBu_r) <matplotlib.image.AxesImage at 0x7f0044c82710>

由上圖可以看到，SST的分布是離散的，我們用之前得到的插值函數來平滑SST值，可視化平滑后的SST分布。

# 設置間隔為1°的經緯度坐標網格，用插值函數得到該坐標網格點的SST值 xnew = np.arange(0, 356, 1) ynew = np.arange(-65, 66, 1) znew = f(xnew, ynew) lon_grid, lat_grid = np.meshgrid(xnew-180, ynew) is_on_land = globe.is_land(lat_grid, lon_grid) is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1) # 同樣進行陸地掩膜 znew[is_on_land] = np.nan # 繪制平滑后的SST分布圖 plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r) <matplotlib.image.AxesImage at 0x7f0044d61c50>

我們用同樣的方法繪制樣本0中每個月的SST分布圖，觀察SST分布的變化。

# 繪制0年36個月的海陸掩膜 for i in range(1):plt.figure(figsize=(15, 18))for j in range(36):x = lony = latxx, yy = np.meshgrid(x, y)z = data.variables['sst'][i, j]f = interpolate.interp2d(x, y, z, kind='cubic')xnew = np.arange(0, 356, 1)ynew = np.arange(-65, 66, 1)znew = f(xnew, ynew)lon_grid, lat_grid = np.meshgrid(xnew-180, ynew)is_on_land = globe.is_land(lat_grid, lon_grid)is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1)znew[is_on_land] = np.nanplt.subplot(9, 4, j+1)plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r)plt.title('sst - year:{}, month:{}'.format(i+1, j+1))

可以看到，SST在每12個月中的前4個月和后4個月都較低，在中間4個月時較高，這說明，海表溫度在春季和冬季較低，在夏季和秋季呈現逐漸升高到最高點然后逐漸降低的變化趨勢，這與我們的認知常識是相符的。

大家也可以用同樣的方法觀察分析其它三個氣象特征的變化趨勢。

以CMIP數據集為例，進行缺失值分析

# 讀取CMIP數據# 存放數據的路徑 path = '/kaggle/input/ninoprediction/' data = netCDF4.Dataset(path + 'CMIP_train.nc') label = netCDF4.Dataset(path + 'CMIP_label.nc') label = np.array(label.variables['nino']) # 獲得陸地的掩膜 lon_grid, lat_grid = np.meshgrid(x-180, y) is_on_land = globe.is_land(lat_grid, lon_grid) is_on_land = np.concatenate([is_on_land[:, x >= 180], is_on_land[:, x < 180]], axis=1) mask = np.zeros(data.variables['sst'].shape, dtype=int) mask[:, :, :, :] = is_on_land[np.newaxis, np.newaxis, :, :] # 查看SST特征的缺失值數量 name = 'sst' data_ = np.array(data.variables[name]) before_nan = np.sum(np.isnan(data_)) print('before:', before_nan) before: 0 # 查看T300特征的缺失值數量 name = 't300' data_ = np.array(data.variables[name]) before_nan = np.sum(np.isnan(data_)) print('before:', before_nan) before: 3055032 # 查看Va特征的缺失值數量 name = 'va' data_ = np.array(data.variables[name]) before_nan = np.sum(np.isnan(data_)) print('before:', before_nan) before: 13921123 # 查看Ua特征的缺失值數量 name = 'ua' data_ = np.array(data.variables[name]) before_nan = np.sum(np.isnan(data_)) print('before:', before_nan) before: 13921123

四個氣象特征中，SST特征不存在缺失值，Va和Ua特征中的缺失值數量最多。

接下來以Ua特征為例，可視化分析缺失值的情況。

# 統計每年每月中Ua特征的缺失值數量 m = np.zeros(data_.shape[0:2]) for i in range(data_.shape[0]):for j in range(data_.shape[1]):if np.sum(np.isnan(data_[i][j])) != 0:m[i, j] = np.sum(np.isnan(data_[i, j])) # 計算每一年的缺失值 before = np.sum(m, axis=1) # 可視化每一年的缺失值數量 plt.plot(before, 'k') plt.ylabel('nan count') plt.xlabel('year') plt.show()

可以看到在某些年份中存在較多缺失值。

# 查看Ua特征中存在缺失值的年數 len(np.where(before!=0)[0]) 755

我們取樣本1900來觀察Ua特征的分布。

# 可視化樣本1900中Ua特征的分布 plt.imshow(data_[1900, 0][::-1, :]) <matplotlib.image.AxesImage at 0x7f003ff92a50>

上圖中白色部分即為缺失值，可以看到，缺失值多數分布在陸地上，我們將陸地部分進行填充，觀察填充后Ua的分布。

# 將陸地位置填0 data_[mask==1] = 0 # 可視化填充后樣本1900中Ua特征的分布 plt.imshow(data_[1900, 0][::-1, :]) <matplotlib.image.AxesImage at 0x7f003ff1bc90>

對陸地部分進行填充后缺失值數量大大減少。

# 統計填充后缺失值的數量 after_nan = np.sum(np.isnan(data_)) print('before: %d \nafter: %d \npercentage: %f'%(before_nan, after_nan, 1 - float(after_nan) / before_nan)) before: 13921123 after: 2440742 percentage: 0.824673

陸地部分填充處理了82%的缺失值。

# 可視化填充后每一年的缺失值數量 m = np.zeros(data_.shape[0: 2]) for i in range(data_.shape[0]):for j in range(data_.shape[1]):if np.sum(np.isnan(data_[i, j])) != 0:m[i, j] = np.sum(np.isnan(data_[i, j])) after = np.sum(m, axis=1) plt.plot(after, 'k') plt.ylabel('nan count') plt.xlabel('year') plt.show()

# 對比填充前后每一年缺失值的數量 plt.plot(before, 'k') plt.plot(after, 'r') plt.legend(['before', 'after']) plt.title(name) plt.ylabel('nan count') plt.xlabel('year') plt.show()

溫度場和風場可視化

在氣候問題中，溫度與風向往往是密切相關的。當溫度越高時，氣壓越低，空氣向上流動，溫度越低時，氣壓越高，空氣向下流動，于是溫度高的地方上方的空氣就會向溫度低的地方流動，形成風。因此在分析氣候問題時，我們往往會把溫度和風向放在一起進行可視化。

如何把風向可視化呢？這里我們要用到plt.quiver()這個函數。

plt.quiver()用于繪制二維的向量場，主要輸入參數有：

• X：向量起始點的X軸坐標
• Y：向量起始點的Y軸坐標
• U：向量的X軸分量
• V：向量的Y軸分量

詳細的用法可以參考官網：https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html

# 對溫度場SST進行插值，得到插值函數 x = lon y = lat xx, yy = np.meshgrid(x, y) z = data.variables['sst'][0, 0] f = interpolate.interp2d(x, y, z, kind='cubic') # 獲得陸地掩膜 lon_grid, lat_grid = np.meshgrid(x-180, y) is_on_land = globe.is_land(lat_grid, lon_grid) is_on_land = np.concatenate([is_on_land[:,x>=180], is_on_land[:,x<180]], axis=1)# 對Ua和Va進行陸地掩膜 ua = data.variables['ua'][0, 0] ua[is_on_land] = np.nan va = data.variables['va'][0, 0] va[is_on_land] = np.nan # 插值后生成平滑的SST分布 xnew = np.arange(0, 356, 1) ynew = np.arange(-65, 66, 1) znew = f(xnew, ynew)# 對平滑后的SST進行陸地掩膜 lon_grid, lat_grid = np.meshgrid(xnew-180, ynew) is_on_land = globe.is_land(lat_grid, lon_grid) is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1) znew[is_on_land] = np.nan# 繪制溫度場 plt.figure(figsize=(15, 10)) plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r) plt.colorbar(orientation='horizontal') # 顯示水平顏色條 # 繪制風向場，其實這里準確來說繪制的是風向異常的向量，而非實際的風向 plt.quiver(lon, lat+65, ua[::-1, :], va[::-1, :], alpha=0.8) # 在坐標(lon, lat)處繪制與sqrt(ua^2, va^2)成正比長度的箭頭 plt.title('year0-month0: SST/UA-VA') plt.show()

從上圖中可以看出，溫度異常SST在0值附近時沒有明顯的風向異常，而在其他區域風向異常通常由SST值大的地方指向SST值小的地方。

ENSO現象是指在溫度場上赤道東太平洋溫度持續異常增暖，在風向場上熱帶東太平洋與熱帶西太平洋氣壓變化（表現為風向）相反的現象。在上圖這個樣本中沒有出現ENSO現象，大家可以用同樣的方法繪制并觀察存在ENSO現象（Nino3.4指數連續5個月超過0.5℃）的樣本的溫度和風場。

總結

在以上的數據分析中，不難看出我們在分析氣象問題時，采用的仍然是通用的數據分析思路：分析標簽 -> 分析特征（包括特征分布、特征與特征的關系、特征與標簽的關系） -> 分析數據的基本情況（包括缺失值、異常值、重復值等）。這個思路大家可以靈活應用到各種問題的分析中去，不至于拿到數據后無從下手。

通過以上的數據分析，我們可以得到以下結論：

重疊部分的標簽是相同的，為了增加數據量，我們可以從每條數據中取固定的12個月拼接起來，用滑窗構建訓練數據集。

SST特征中沒有缺失值，在其他特征中，缺失值基本上在陸地部分，將陸地部分用0填充可以解決絕大部分的缺失值。

現在我們回到開篇的學習目標中提出的第一個問題：能否構造新的特征？目前的答案是不能。從各位TOP選手的方案以及相關的ENSO預測的論文來看，大家的目光都聚焦在如何構建模型上，而鮮少有人會去構造新的特征。這與其說是不能，更不如說是不必要，因為一般我們構造新的特征是為了從給出的特征中得到與預測目標更相關的特征，由此提高模型的學習效果，但是就本賽題而言，構造統計特征或者其他新的特征收效不高，我們更希望通過模型來挖掘給定數據之間在時間和空間上的依賴關系。

于是，解題的重點就落到了第二個問題上：選擇或設計什么樣的模型進行預測？在接下來的三個任務中，我們就要來學習TOP選手們的建模方案。

作業

在學習TOP選手的方案之前，請嘗試思考你會如何構造模型。

參考文獻

優秀選手的EDA分享：https://gitee.com/Little_Six/aiweather-ocean-forecasts/blob/master/code/EDA.ipynb?spm=5176.21852664.0.0.4aab7f46PqZUih&file=EDA.ipynb

Climate Data Guide：https://climatedataguide.ucar.edu/climate-data/nino-sst-indices-nino-12-3-34-4-oni-and-tni

中國氣象局國家氣候中心：https://cmdp.ncc-cma.net/pred/cn_enso.php?product=cn_enso_nino_indices

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的【算法竞赛学习】气象海洋预测-Task2 数据分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。