从xarray走向netCDF处理(一)：数据结构及数据读取
想如今气象数据netCDF(.nc)为盛，用者甚多，初学者见之仰天长啸，倘若再由Python经手，netCDF4-python，Iris，xarray，UV-CDAT选择众多，劳心伤神事小，逼出选择困难症事大。多番比对，选用xarray，解查安抚，化繁为简，最为称心。
说人话就是，经学前班大队长亲测利用Python中的xarray库处理nc数据非常方便。
安装
xarray的安装依旧推荐使用conda，还不会的小伙伴移步：一文教你解决Python所有安装配置
conda install xarray
在终端里输入如上命令，之后输入y，等待安装结束就好了。
数据结构
xarray有两大数据类型：DataArray、Dataset。
DataArray
一个带有标签的多维数组，它有如下几个重要的属性
values获取数组的具体数值dims获取维度的名字，如('x', 'y', 'z')coords获取一个类似于字典的结果，里面包含各个坐标attrs获取原始数据的属性，比如变量的名字、单位等
Dataset
Dataset可以简单的理解为由多个DataArray组成的集合，它有如下几个重要的属性
dims获取维度的名字，结果类似于字典，如{'x': 6, 'y': 6, 'time': 8}data_vars获取物理量的名字coords获取一个类似于字典的结果，里面包含各个坐标attrs获取原始数据的属性，比如变量的名字、单位等
数据结构图示
数据类型的使用
读取数据:
- xarray.open_dataset()读取Dataset类型数据，即能读取多个物理量。
- xarray.open_dataarray()读取DataArray类型数据，即只能读取单个物理量。
如果nc文件中含有多个物理量，用open_dataarray()读取会报错，因此建议统一都用open_dataset()来读取文件。
提取物理量
从文件中读取数据ds = xarray.open_dataset()
假如数据中含有一个名为var的物理量可以通过ds.var或ds[var]来获取
实例
此处使用的是ERA-Interim中2018年的月数据，包含10米的径向风、纬向风和2米气温,在ECMWF注册过的都可以直接下载。下面代码部分建议大家横屏阅读。
>>>import xarray as xr
# 由于数据包含了多个物理量（u10，v10，t2m），所以要用open_dataset来读取数据
>>>ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
# ds的类型为Dataset，里面包含u10，v10，t2m三个物理量，每个物理量都有经度、纬度、时间三个坐标系
>>>ds
<xarray.Dataset>
Dimensions: (latitude: 241, longitude: 480, time: 12)
Coordinates:
* longitude (longitude) float32 0.0 0.75 1.5 2.25 .. 357.75 358.5 359.25
* latitude (latitude) float32 90.0 89.25 88.5 87.75 .. -88.5 -89.25 -90.0
* time (time) datetime64[ns] 2018-01-01 2018-02-01 .. 2018-12-01
Data variables:
u10 (time, latitude, longitude) float32 ..
v10 (time, latitude, longitude) float32 ..
t2m (time, latitude, longitude) float32 ..
Attributes:
Conventions: CF-1.6
history: 2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...
# 取出ds中名为t2m的物理量,可以看到它的维度，坐标系，以及t2m有单位和名字两个属性
>>>ds['t2m']
<xarray.DataArray 't2m' (time: 12, latitude: 241, longitude: 480)>
[13e88160 values with dtype=float32]
Coordinates:
* longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
* latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
* time (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
units: K
long_name: 2 metre temperatur通过xarray可以清晰的了解nc数据中的维度、坐标、物理量以及各种属性等信息。清晰的数据结构是准确、高效地分析数据的基础
简单的可视化
了解完数据结构，再来看一看数据可视化的结果。xarray封装了matplotlib的部分绘图函数，一行代码就可以将数据画出来，不过作为一个负责的公众号，还是用cartopy顺带加载了地图。
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
# 数据读取及时间平均处理
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
temp = (ds['t2m'] - 273.15).mean(dim='time') #把温度转换为℃并对其时间纬求平均
temp.attrs['units'] = 'deg C' #温度单位转换为℃
# 创建画图空间
proj = ccrs.PlateCarree() #创建投影
fig = plt.figure(figsize=(9,6)) #创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) #子图
# 设置地图属性:加载国界、海岸线、河流、湖泊
ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)
ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1)
ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1)
ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1)
# 设置网格点属性
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False #关闭顶端标签
gl.ylabels_right = False #关闭右侧标签
gl.xformatter = LONGITUDE_FORMATTER #x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER #y轴设为纬度格式
# 设置colorbar
cbar_kwargs = {
'orientation': 'horizontal',
'label': '2m temperature (℃)',
'shrink': 0.8,
'ticks': np.arange(-30,30+5,5)
}
# 画图
levels = np.arange(-30,30+1,1)
temp.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
fig.show()后面还有一些与nc数据相关的操作和画图，感兴趣的伙计们赶快关注我们。
本文数据
https://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/
腾讯云开发者
