在气象领域中NetCDF(nc)数据使用十分普遍。Python也有很多处理nc数据的库:netCDF4-python,Iris,xarray,UV-CDAT。在众多库中,xarray是我觉得最好用的一个,有大量的开发人员在不断的维护更新,还有很多人基于此之上做了很多的二次开发,功能完善,操作简单,nc数据处理必备。

安装

xarray的安装依旧推荐使用conda

conda install xarray

在终端里输入如上命令,之后输入y,等待安装结束就好了。

数据结构

xarray有两大数据类型:DataArray、Dataset。

DataArray

一个带有标签的多维数组,它有如下几个重要的属性 - values 获取数组的具体数值 - dims 获取维度的名字,如('x', 'y', 'z') - coords 获取一个类似于字典的结果,里面包含各个坐标 - attrs 获取原始数据的属性,比如变量的名字、单位等

Dataset

如果为了方便理解,Dataset可以简单的理解为由多个DataArray组成的集合,它比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.vards[var]来获取

实例

此处使用的是ERA-Interim中2018年的月数据,包含10米的径向风、纬向风和2米气温。在ECMWF注册过的都可以直接下载,此处为传送门

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
>>>import xarray as xr

# 由于数据包含了多个物理量(u10,v10,t2m),所以要用open_dataset来读取数据
>>>ds = xr.open_dataset('EC-Interim_monthly_2018.nc')

# ds的类型为Dataset,里面包含u10,v10,t2m三个物理量
# 每个物理量都有经度、纬度、时间三个坐标系
# 还可以从Attributes中看到数据的一些原始属性
>>>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)>
[1388160 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 temperature

通过xarray可以清晰的了解nc数据中的维度、坐标、物理量以及各种属性等信息。清晰的数据结构是准确、高效地分析数据的基础

简单的可视化

了解完数据结构,再来看一看数据可视化的结果。xarray封装了matplotlib的部分绘图函数,一行代码就可以将数据画出来,不过作为有良心的公众号,还是用cartopy加载了地图。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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()