1.“AI Earth”人工智能翻新挑战赛：助力精准气象和陆地预测 Baseline[1]、NetCDF4 应用教学、Xarray 应用教学，针对气象畛域.nc 文件读取解决
较量官网：https://tianchi.aliyun.com/specials/promotion/aiearth2021?spm…
1.1 背景形容
聚焦寰球大气陆地钻研前沿方向，将人工智能技术利用到天气气象预测畛域中，进步极其灾害性天气的预报程度，已成为整个行业钻研的热点方向。产生在寒带太平洋上的厄尔尼诺 - 北方涛动 (ENSO) 景象是地球上最强、最显著的年际气象信号，常常会引发洪涝、干旱、低温、雪灾等极其事件，2020 年底我国夏季极寒也与 ENSO 非亲非故。对于 ENSO 的预测，气象能源模式耗费计算资源大且存在秋季预测阻碍。基于历史气象观测和模仿数据，利用 T 时刻过来 12 个月 (蕴含 T 时刻) 的时空序列（气象因子），能够构建预测 ENSO 的深度学习模型，预测将来 1 -24 个月的 Nino3.4 指数，这对极其天气与气象事件的预测具备重要意义。
本数据集由气象与利用前沿研究院 ICAR 提供。数据包含 CMIP5/ 6 模式的历史模仿数据和美国 SODA 模式重建的近 100 多年历史观测异化数据。每个样本蕴含以下气象及时空变量：海表温度异样(SST)，热含量异样(T300)，纬向风异样（Ua），经向风异样（Va），数据维度为（year,month,lat,lon）。训练数据提供对应月份的 Nino3.4 index 标签数据。测试用的初始场数据为国内多个陆地材料异化后果提供的随机抽取的 n 段 12 个工夫序列，数据格式采纳 NPY 格局保留。
1.2 历史气象观测和模仿数据集数据阐明
我的项目链接以及码源见文末
- 训练数据
每个数据样本第一维度（year）表征数据所对应起始年份，对于 CMIP 数据共 4645 年，其中 1 -2265 为 CMIP6 中 15 个模式提供的 151 年的历史模仿数据（总共：151 年 15 个模式 =2265）；2266-4645 为 CMIP5 中 17 个模式提供的 140 年的历史模仿数据（总共：140 年 17 个模式 =2380）。对于历史观测异化数据为美国提供的 SODA 数据。
其中每个样本第二维度（month）表征数据对应的月份，对于训练数据均为 36，对应的从以后年份开始间断三年数据（从 1 月开始，共 36 月），比方：
SODA_train.nc 中 [0,0:36,:,:] 为第 1 - 第 3 年逐月的历史观测数据；
SODA_train.nc 中 [1,0:36,:,:] 为第 2 - 第 4 年逐月的历史观测数据；
…,
SODA_train.nc 中 [99,0:36,:,:] 为第 100-102 年逐月的历史观测数据。
和
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 度），所有数据的经纬度范畴雷同。
- 训练数据标签
标签数据为 Nino3.4 SST 异样指数，数据维度为（year,month）。
CMIP(SODA)_train.nc 对应的标签数据以后时刻 Nino3.4 SST 异样指数的三个月滑动平均值，因而数据维度与维度介绍同训练数据统一
注：三个月滑动平均值为以后月与将来两个月的平均值。
- 测试数据
测试用的初始场（输出）数据为国内多个陆地材料异化后果提供的随机抽取的 n 段 12 个工夫序列，数据格式采纳 NPY 格局保留，维度为（12，lat，lon, 4）,12 为 t 时刻及过来 11 个时刻，4 为预测因子，并依照 SST,T300,Ua,Va 的程序寄存。
测试集文件序列的命名规定：test_编号_起始月份_终止月份.npy，如 test_00001_01_12_.npy。
- 数据 (Netcdf 文件) 读取办法
（1）https://www.giss.nasa.gov/tools/panoply/ panoply 可视化文件
（2）Python 中 xarray/netCDF4 库
1.3 评估指标如下：
$Score = \frac{2}{3} \times accskill – RMSE$
其中 $accskill$ 为相关性技巧评分，计算形式如下：
$accskill = \sum_{i=1}^{24} a \times ln(i) \times cor_i \ (i \leq 4, a = 1.5; 5 \leq i \leq 11, a = 2; 12 \leq i \leq 18, a = 3; 19 \leq i, a = 4)$
能够看出，月份 $i$ 减少时系数 $a$ 也增大，也就是说，模型能精确预测的工夫越长，评分就越高。
是对于 $N$ 个测试集样本在时刻 $i$ 的预测值与理论值的相关系数，计算公式如下：
$cor_i = \frac{\sum_{j=1}^N(y_{truej}-\bar{y}{true})(y{predj}-\bar{y}{pred})}{\sqrt{\sum(y{truej}-\bar{y}{true})^2\sum(y{predj}-\bar{y}{pred})^2}}$
其中 $y{truej}$ 为时刻 $i$ 样本 $j$ 的理论 Nino3.4 指数，$\bar{y}{true}$ 为该时刻 $N$ 个测试集样本的 Nino3.4 指数的均值，$y{predj}$ 为时刻 $i$ 样本 $j$ 的预测 Nino3.4 指数，$\bar{y}_{pred}$ 为该时刻 $N$ 个测试集样本的预测 Nino3.4 指数的均值。
为 24 个月份的累计均方根误差，计算公式为：
$RMSE = \sum_{i=1}^{24}rmse_i \ rmse = \sqrt{\frac{1}{N}\sum_{j=1}^N(y_{truej}-y_{predj})^2}$
ENSO 景象会在世界大部分地区引起极其天气，对寰球的天气、气象以及粮食产量具备重要的影响，精确预测 ENSO，是进步东亚和寰球气象预测程度和防灾减灾的要害。Nino3.4 指数是 ENSO 景象监测的一个重要指标，它是指 Nino3.4 区 (170°W – 120°W，5°S – 5°N) 的均匀海温距平指数，用于反馈海表温度异样，若 Nino3.4 指数间断 5 个月超过 0.5℃就断定为一次 ENSO 事件。本赛题的指标，就是基于历史气象观测和模式模仿数据，利用 T 时刻过来 12 个月 (蕴含 T 时刻) 的时空序列，预测将来 1 – 24 个月的 Nino3.4 指数。
# !unzip -d input /home/aistudio/data/data221707/enso_final_test_data_B_labels.zip2. 气象陆地预测 - 气象数据分析常用工具
气象迷信中的数据通常蕴含多个维度，本赛题中给出的数据就蕴含年、月、经度、纬度四个维度，为了便于数据的读取和操作，气象数据通常采纳 netCDF 文件来存储，文件后缀为.nc。
对于以 netCDF 文件存储的气象数据，有两个罕用的数据分析库，即 NetCDF4 和 Xarray。首先将学习这两个库的根本对象和基本操作，把握用这两个库读取和解决气象数据的根本办法。
2.1 NetCDF4 应用教学
官网文档
NetCDF4 是 NetCDF C 库的 Python 模块，反对 Groups、Dimensions、Variables 和 Attributes 等对象类型及其相干操作。
- 装置 NetCDF4
!pip install netCDF42.1.1 创立、关上和敞开 netCDF 文件
NetCDF4 能够通过调用 Dataset 创立 netCDF 文件或关上已存在的文件，并通过查看 data_model 属性确定文件的格局。须要留神创立或关上文件后要先敞开文件能力再次调用 Dataset 关上文件。
- 创立 netCDF 文件
import netCDF4 as nc
from netCDF4 import Dataset
# Dataset 蕴含三个输出参数：文件名，模式（其中 'w', 'r+', 'a' 为可写入模式），文件格式
test = Dataset('test.nc', 'r+', 'NETCDF4')
# test = Dataset('test.nc', 'w', 'NETCDF4')# - 关上已存在的 netCDF 文件
# 关上训练样本中的 SODA 数据
soda= Dataset("SODA_train.nc", "w", format="NETCDF4")
# soda = Dataset('SODA_train.nc')# - 查看文件格式
print(soda.data_model)NETCDF4
# - 敞开 netCDF 文件
soda.close()2.1.2 Groups
NetCDF4 反对按层级的组（Groups）来组织数据，相似于文件系统中的目录，Groups 中能够蕴含 Variables、Dimenions、Attributes 对象以及其余 Groups 对象，Dataset 会创立一个非凡的 Groups，称为根组（Root Group），相似于根目录，应用 Dataset.createGroup 办法创立的组都蕴含在根组中。
- 创立 Groups
# 承受一个字符串参数作为 Group 名称
group1 = test.createGroup('group1')
group2 = test.createGroup('group2')# - 查看文件中的所有 Groups
# 返回一个 Group 字典
test.groups{'group1': <class 'netCDF4._netCDF4.Group'>
group /group1:
dimensions(sizes):
variables(dimensions):
groups: group11,
'group2': <class 'netCDF4._netCDF4.Group'>
group /group2:
dimensions(sizes):
variables(dimensions):
groups: group21}
# - Groups 嵌套
# 在 group1 和 group2 下别离再创立一个 Group
group1_1 = test.createGroup('group1/group11')
group2_1 = test.createGroup('group2/group21')
test.groups
test.groups.values()dict_values([<class 'netCDF4._netCDF4.Group'>
group /group1:
dimensions(sizes):
variables(dimensions):
groups: group11, <class 'netCDF4._netCDF4.Group'>
group /group2:
dimensions(sizes):
variables(dimensions):
groups: group21])
# - 遍历查看所有 Groups
# 定义一个生成器函数用来遍历所有目录树
def walktree(top):
values = top.groups.values()
yield values
for value in top.groups.values():
for children in walktree(value):
yield children
for groups in walktree(test):
for group in groups:
print(group)<class 'netCDF4._netCDF4.Group'>
group /group1:
dimensions(sizes):
variables(dimensions):
groups: group11
<class 'netCDF4._netCDF4.Group'>
group /group2:
dimensions(sizes):
variables(dimensions):
groups: group21
<class 'netCDF4._netCDF4.Group'>
group /group1/group11:
dimensions(sizes):
variables(dimensions): float64 ftemp(time, level, lat, lon)
groups:
<class 'netCDF4._netCDF4.Group'>
group /group2/group21:
dimensions(sizes):
variables(dimensions):
groups:
2.1.3 Dimensions
NetCDF4 用维度来定义各个变量的大小，例如本赛题中训练样本的第二维度 month 就是一个维度对象，每个样本蕴含 36 个月的数据，因而 month 维度内的变量的大小就是 36。变量是蕴含在维度中的，因而在创立每个变量时要先创立其所在的维度。
- 创立 Dimensions
Dataset.createDimension 办法承受两个参数：维度名称，维度大小。维度大小设置为 None 或 0 时示意无穷维度。
# 创立无穷维度
level = test.createDimension('level', None)
time = test.createDimension('time', None)
# 创立无限维度
lat = test.createDimension('lat', 180)
lon = test.createDimension('lon', 360)---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_2348/1064279341.py in <module>
1 # 创立无穷维度
----> 2 level = test.createDimension('level', None)
3 time = test.createDimension('time', None)
4 # 创立无限维度
5 lat = test.createDimension('lat', 180)
src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dataset.createDimension()
src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dimension.__init__()
src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4._ensure_nc_success()
RuntimeError: NetCDF: String match to name in use
# - 查看 Dimensions
test.dimensions# 查看维度大小
print(len(lon))
# Dimension 对象存储在字典中
print(level)
# 判断维度是否是无穷
print(time.isunlimited())
print(lat.isunlimited())2.1.4 Variables
NetCDF4 的 Variables 对象相似于 Numpy 中的多维数组，不同的是，NetCDF4 的 Variables 变量能够存储在无穷维度中。
- 创立 Variables
Dataset.createVariable 办法承受的参数为：变量名，变量的数据类型，变量所在的维度。
变量的无效数据类型包含：’f4′(32 位浮点数)、’f8′(64 位浮点数)、’i1′(8 位有符号整型)、’i2′(16 位有符号整型)、’i4′(32 位有符号整型)、’i8′(64 位有符号整型)、’u1′(8 位无符号整型)、’u2′(16 位无符号整型)、’u4′(32 位无符号整型)、’u8′(64 位无符号整型)、’s1′(单个字符)。
# 创立单个维度上的变量
times = test.createVariable('time', 'f8', ('time',))
levels = test.createVariable('level', 'i4', ('level',))
lats = test.createVariable('lat', 'f4', ('lat',))
lons = test.createVariable('lon', 'f4', ('lon',))
# 创立多个维度上的变量
temp = test.createVariable('temp', 'f4', ('time', 'level', 'lat', 'lon'))# - 查看 Variables
print(temp)# 通过门路的形式在 Group 中创立变量
ftemp = test.createVariable('/group1/group11/ftemp', 'f8', ('time', 'level', 'lat', 'lon'))
# 能够通过门路查看变量
print(test['/group1/group11/ftemp'])
print(test['/group1/group11'])
# 查看文件中的所有变量
print(test.variables)2.1.5 Attributes
Attributes 对象用于存储对文件或维变量的形容信息，netcdf 文件中蕴含两种属性：全局属性和变量属性。全局属性提供 Groups 或整个文件对象的信息，变量属性提供 Variables 对象的信息，属性的名称能够本人设置，上面例子中的 description 和 history 等都是自定义的属性名称。
import time
# 设置对文件的形容
test.description = 'bogus example script'
# 设置文件的历史信息
test.history = 'Created' + time.ctime(time.time())
# 设置文件的起源信息
test.source = 'netCDF4 python module tutorial'
# 设置变量属性
lats.units = 'degrees north'
lons.units = 'degrees east'
levels.units = 'hPa'
temp.units = 'K'
times.units = 'hours since 0001-01-01 00:00:00.0'
times.calendar = 'gregorian'# 查看文件属性名称
print(test.ncattrs())
# 查看变量属性名称
print(test['lat'].ncattrs())
print(test['time'].ncattrs())
# 查看文件属性
for name in test.ncattrs():
print('Global attr {} = {}'.format(name, getattr(test, name)))2.1.6 写入或读取变量数据
相似于数组，能够通过切片的形式向变量中写入或读取数据。
- 向变量中写入数据
from numpy.random import uniform
nlats = len(test.dimensions['lat'])
nlons = len(test.dimensions['lon'])
print('temp shape before adding data = {}'.format(temp.shape))# 无穷维度的大小会随着写入的数据的大小主动扩大
temp[0:5, 0:10, :, :] = uniform(size=(5, 10, nlats, nlons))
print('temp shape after adding data = {}'.format(temp.shape))
print('levels shape after adding pressure data = {}'.format(levels.shape))- 读取变量中的数据
print(temp[1, 5, 100, 300])
print(temp[1, 5, 10:20, 100:110].shape)
# 能够用 start:stop:step 的模式进行切片
print(temp[1, 5, 10, 100:110:2])2.1.7 利用
咱们尝试用 NetCDF4 来操作一下训练样本中的 SODA 数据。
# 关上 SODA 文件
soda = Dataset('SODA_train.nc')
# 查看文件格式
print('SODA 文件格式：', soda.data_model)
# 查看文件中蕴含的对象
print(soda)# 查看维度和变量
print(soda.dimensions)
print(soda.variables)能够看到，SODA 文件中蕴含 year、month、lat、lon 四个维度，维度大小别离是 100、36、24 和 72，蕴含 sst、t300、ua、va 四个变量，每个变量都定义在 (year, month, lat, lon) 维度上。
# # 读取每个变量中的数据
# soda_sst = soda['sst'][:]
# print(soda_sst[1, 1, 1, 1])
# soda_t300 = soda['t300'][:]
# print(soda_t300[1, 2, 12:24, 36])
# soda_ua = soda['ua'][:]
# print(soda_ua[1, 2, 12:24:2, 36:38])
# soda_va = soda['va'][:]
# print(soda_va[5:10, 0:12, 12, 36])
# # 敞开文件
# soda.close()2.2 Xarray 应用教学
官网文档
Xarray 是一个开源的 Python 库，反对在相似 Numpy 的多维数组上引入维度、坐标和属性标记并能够间接应用标记的名称进行相干操作，可能读写 netcdf 文件并进行进一步的数据分析和可视化。
Xarray 有两个根本的数据结构：DataArray 和 Dataset，这两个数据结构都是在多维数组上建设的，其中 DataArray 用于标记的实现，Dataset 则是一个相似于字典的 DataArray 容器。
装置 Xarray 要求满足以下依赖包：
- Python(3.7+)
- setuptools(40.4+)
- Numpy(1.17+)
- Pandas(1.0+)
!pip install xarray2.2.1 创立 DataArray
xr.DataArray 承受三个输出参数：数组，维度，坐标。其中维度为数组的维度名称，坐标以字典的模式给维度赋予坐标标签。
import numpy as np
import pandas as pd
import xarray as xr
# 创立一个 2x3 的数组，将维度命名为 'x' 和 'y'，并赋予 'x' 维度 10 和 20 两个坐标标签
data = xr.DataArray(np.random.randn(2, 3), dims=('x', 'y'), coords={'x': [10, 20]})
data
# 也能够用 Pandas 的 Series 或 DataFrame 数据创立 DataArray。# index 名称会主动转换成坐标标签
xr.DataArray(pd.Series(range(3), index=list('abc'), name='foo'))# 查看数据
print(data.values)
# 查看维度
print(data.dims)
# 查看坐标
print(data.coords)
# 能够用 data.attrs 字典来存储任意元数据
print(data.attrs)- 索引
Xarray 反对四种索引形式。
# 通过地位索引，相似于 numpy
print(data[0, :], '\n')
# 通过坐标标签索引
print(data.loc[10], '\n')
# 通过维度名称和地位索引，isel 示意 "integer select"
print(data.isel(x=0), '\n')
# 通过维度名称和坐标标签索引，sel 示意 "select"
print(data.sel(x=10), '\n')- 属性
和 NetCDF4 一样，Xarray 也反对自定义 DataArray 或标记的属性形容。
# 设置 DataArray 的属性
data.attrs['long_name'] = 'random velocity'
data.attrs['units'] = 'metres/sec'
data.attrs['description'] ='A random variable created as an example'
data.attrs['ramdom_attribute'] = 123
# 查看属性
print(data.attrs)# 设置维度标记的属性形容
data.x.attrs['units'] ='x units'
print('Attributes of x dimension:', data.x.attrs, '\n')- 计算
DataArray 的计算形式相似于 numpy ndarray。
data + 10np.sin(data)data.Tdata.sum()能够间接应用维度名称进行聚合操作。
data.mean(dim='x')DataArray 之间的计算操作能够依据维度名称进行播送。
a = xr.DataArray(np.random.randn(3), [data.coords['y']])
b = xr.DataArray(np.random.randn(4), dims='z')
print(a, '\n')
print(b, '\n')
print(a+b, '\n')data - data.Tdata[:-1] - data[:1]- GroupBy
Xarray 反对应用相似于 Pandas 的 API 进行分组操作。
labels = xr.DataArray(['E', 'F', 'E'], [data.coords['y']], name='labels')
labels# 将 data 的 y 坐标对齐 labels 后按 labels 的值分组求均值
data.groupby(labels).mean('y')# 将 data 的 y 坐标按 labels 分组后减去组内的最小值
data.groupby(labels).map(lambda x: x - x.min())2.2.2 绘图
Xarray 反对简略不便的可视化操作，这里只做简略的介绍，更多的绘图办法感兴趣的同学们能够自行去摸索。
%matplotlib inline
data.plot()<matplotlib.collections.QuadMesh at 0x7f5bf77c4490>
2.2.3 与 Pandas 对象相互转换
Xarray 能够不便地转换成 Pandas 的 Series 或 DataFrame，也能够由 Pandas 对象转换回 Xarray。
# 转换成 Pandas 的 Series
series = data.to_series()
seriesx y
10 0 -0.188398
1 0.612211
2 0.741517
20 0 -0.547693
1 -0.021661
2 0.183370
dtype: float64
# Series 转换成 Xarray
series.to_xarray()# 转换成 Pandas 的 DataFrame
df = data.to_dataframe(name='colname')
df2.2.4 Dataset
Dataset 是一个相似于字典的 DataArray 的容器，能够看作是一个具备多为构造的 DataFrame。比照 NetCDF4 库中的 Dataset，咱们能够发现两者的作用是类似的，都是作为容器用来存储其余的对象。
# 创立一个 Dataset，其中蕴含三个 DataArray
ds = xr.Dataset({'foo': data, 'bar': ('x', [1, 2]), 'baz': np.pi})
ds能够通过字典的形式或者点索引的形式来查看 DataArray，然而只有采纳字典形式时才能够进行赋值。
# 通过字典形式查看 DataArray
print(ds['foo'], '\n')
# 通过点索引的形式查看 DataArray
print(ds.foo)<xarray.DataArray 'foo' (x: 2, y: 3)>
array([[-0.18839763, 0.6122107 , 0.74151684],
[-0.547693 , -0.02166147, 0.18336995]])
Coordinates:
* x (x) int64 10 20
Dimensions without coordinates: y
Attributes:
long_name: random velocity
units: metres/sec
description: A random variable created as an example
ramdom_attribute: 123
<xarray.DataArray 'foo' (x: 2, y: 3)>
array([[-0.18839763, 0.6122107 , 0.74151684],
[-0.547693 , -0.02166147, 0.18336995]])
Coordinates:
* x (x) int64 10 20
Dimensions without coordinates: y
Attributes:
long_name: random velocity
units: metres/sec
description: A random variable created as an example
ramdom_attribute: 123
同样能够通过坐标标记来索引。
ds.bar.sel(x=10)2.2.5 读 / 写 netCDF 文件
# 写入到 netcdf 文件
ds.to_netcdf('xarray_test.nc')
# 读取已存在的 netcdf 文件
xr.open_dataset('xarray_test.nc')2.2.6 利用
尝试用 Xarray 来操作一下训练样本中的 SODA 数据。
# 关上 SODA 文件
soda = xr.open_dataset('SODA_train.nc')
# 查看文件属性
print(soda.attrs)
# 查看文件中蕴含的对象
print(soda){}
<xarray.Dataset>
Dimensions: ()
Data variables:
*empty*
# 查看维度和坐标
print(soda.dims)
print(soda.coords)Frozen(SortedKeysDict({}))
Coordinates:
*empty*
# # 读取数据
# soda_sst = soda['sst']
# print(soda_sst[1, 1, 1, 1], '\n')
# soda_t300 = soda['t300']
# print(soda_t300[1, 2, 12:24, 36], '\n')
# soda_ua = soda['ua']
# print(soda_ua[1, 2, 12:24:2, 36:38], '\n')
# soda_va = soda['va']
# print(soda_va[5:10, 0:12, 12, 36])我的项目链接以及码源见文末
链接跳转到文末，欢送订阅
更多文章请关注公重号：汀丶人工智能