写点什么

“AI Earth”人工智能创新挑战赛:助力精准气象和海洋预测 Baseline[1]、NetCDF4 使用教学、Xarray 使用教学,针对气象领域.nc 文件读取处理

  • 2023-06-06
    浙江
  • 本文字数:8825 字

    阅读完需:约 29 分钟

“AI Earth”人工智能创新挑战赛:助力精准气象和海洋预测Baseline[1]、NetCDF4使用教学、Xarray 使用教学,针对气象领域.nc文件读取处理

1.“AI Earth”人工智能创新挑战赛:助力精准气象和海洋预测 Baseline[1]、NetCDF4 使用教学、Xarray 使用教学,针对气象领域.nc 文件读取处理

比赛官网:https://tianchi.aliyun.com/specials/promotion/aiearth2021?spm=a2c22.12281976.0.0.4d0d19efK2FngK


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 评估指标如下:


其中为相关性技巧评分,计算方式如下:


可以看出,月份增加时系数也增大,也就是说,模型能准确预测的时间越长,评分就越高。


是对于个测试集样本在时刻的预测值与实际值的相关系数,计算公式如下:


其中为时刻样本的实际 Nino3.4 指数,为该时刻个测试集样本的 Nino3.4 指数的均值,为时刻样本的预测 Nino3.4 指数,为该时刻个测试集样本的预测 Nino3.4 指数的均值。


为 24 个月份的累计均方根误差,计算公式为:


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.zip
复制代码

2.气象海洋预测-气象数据分析常用工具

气象科学中的数据通常包含多个维度,本赛题中给出的数据就包含年、月、经度、纬度四个维度,为了便于数据的读取和操作,气象数据通常采用 netCDF 文件来存储,文件后缀为.nc。


对于以 netCDF 文件存储的气象数据,有两个常用的数据分析库,即 NetCDF4 和 Xarray。首先将学习这两个库的基本对象和基本操作,掌握用这两个库读取和处理气象数据的基本方法。

2.1 NetCDF4 使用教学

官方文档



NetCDF4 是 NetCDF C 库的 Python 模块,支持 Groups、Dimensions、Variables 和 Attributes 等对象类型及其相关操作。


  • 安装 NetCDF4


!pip install netCDF4
复制代码

2.1.1 创建、打开和关闭 netCDF 文件

NetCDF4 可以通过调用 Dataset 创建 netCDF 文件或打开已存在的文件,并通过查看 data_model 属性确定文件的格式。需要注意创建或打开文件后要先关闭文件才能再次调用 Dataset 打开文件。


  • 创建 netCDF 文件


import netCDF4 as ncfrom 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下分别再创建一个Groupgroup1_1 = test.createGroup('group1/group11')group2_1 = test.createGroup('group2/group21')test.groupstest.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
复制代码


# - 查看Dimensionstest.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'))
复制代码


# - 查看Variablesprint(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 xarray
复制代码

2.2.1 创建 DataArray

xr.DataArray 接受三个输入参数:数组,维度,坐标。其中维度为数组的维度名称,坐标以字典的形式给维度赋予坐标标签。


import numpy as npimport pandas as pdimport 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 支持四种索引方式。


# 通过位置索引,类似于numpyprint(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 + 10
复制代码


np.sin(data)
复制代码


data.T
复制代码


data.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.T
复制代码


data[:-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 inlinedata.plot()
复制代码


<matplotlib.collections.QuadMesh at 0x7f5bf77c4490>
复制代码


2.2.3 与 Pandas 对象互相转换

Xarray 可以方便地转换成 Pandas 的 Series 或 DataFrame,也可以由 Pandas 对象转换回 Xarray。


# 转换成Pandas的Seriesseries = data.to_series()series
复制代码


x   y10  0   -0.188398    1    0.612211    2    0.74151720  0   -0.547693    1   -0.021661    2    0.183370dtype: float64
复制代码


# Series转换成Xarrayseries.to_xarray()
复制代码


# 转换成Pandas的DataFramedf = data.to_dataframe(name='colname')df
复制代码

2.2.4 Dataset

Dataset 是一个类似于字典的 DataArray 的容器,可以看作是一个具有多为结构的 DataFrame。对比 NetCDF4 库中的 Dataset,我们可以发现两者的作用是相似的,都是作为容器用来存储其他的对象。


# 创建一个Dataset,其中包含三个DataArrayds = xr.Dataset({'foo': data, 'bar': ('x', [1, 2]), 'baz': np.pi})ds
复制代码


可以通过字典的方式或者点索引的方式来查看 DataArray,但是只有采用字典方式时才可以进行赋值。


# 通过字典方式查看DataArrayprint(ds['foo'], '\n')
# 通过点索引的方式查看DataArrayprint(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 20Dimensions without coordinates: yAttributes:    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 20Dimensions without coordinates: yAttributes: 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])
复制代码

项目链接以及码源见文末

链接跳转到文末,欢迎订阅


更多文章请关注公重号:汀丶人工智能



发布于: 2023-06-06阅读数: 17
用户头像

本博客将不定期更新关于NLP等领域相关知识 2022-01-06 加入

本博客将不定期更新关于机器学习、强化学习、数据挖掘以及NLP等领域相关知识,以及分享自己学习到的知识技能,感谢大家关注!

评论

发布
暂无评论
“AI Earth”人工智能创新挑战赛:助力精准气象和海洋预测Baseline[1]、NetCDF4使用教学、Xarray 使用教学,针对气象领域.nc文件读取处理_人工智能_汀丶人工智能_InfoQ写作社区