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])我的项目链接以及码源见文末
链接跳转到文末，欢送订阅
更多文章请关注公重号：汀丶人工智能
发表回复