Python处理NetCDF数据入门

发表于 分类于 阅读次数： Disqus：

前言

在前面做极值分析、气象数据可视化等的时候，我们主要用的是csv格式存储数据，而在气象科学等领域中，还有一种十分重要的数据格式——NetCDF。在本篇文章里，我们将简要介绍NetCDF以及如何用Python处理NetCDF文件。

NetCDF简介

NetCDF(Network Common Data Form)，即网络通用数据格式，是一种自描述、与机器无关、基于数组的科学数据格式，同时也是支持创建、访问和共享这一数据格式的函数库。NetCDF于1989年由美国大气科学研究大学联盟（UCAR）开发，最新版本为4.x（详见 项目主页）。NetCDF格式是一种开放标准，其经典格式和64位偏移量格式是开放地理空间协会采用的国际标准。

netCDF常用于气候学、气象学、海洋学等领域，如天气预报、气候变化等；也应用于地理信息系统，是许多GIS应用的输入输出格式。NCEP(美国国家环境预报中心)发布的再分析资料，NOAA的CDC(气候数据中心)发布的海洋与大气综合数据集(COADS)均采用NetCDF作为标准。NetCDF的广泛应用，在很大程度上得益于其灵活性，而它的灵活性则要归功于“自描述”的特点。

所谓“自描述”是指，数据本身包含了相关的数据描述信息，它有一个头部，描述文件余下部分的格局，特别是数组数据，连同以名称/值属性形式的任意文件元数据。程序自动判断维和变量等信息的前提条件是数据必须遵循某种约定（convensions）。气象上最常用的约定是气候和预报预定（Climate and Forecast (CF) conventions），对于维、变量、属性有详细的规定，这样以来软件才能通过约定对数据进行正确的判读。

更多关于NetCDF的介绍，参见这里。

Python处理NetCDF

由于NetCDF格式的流行，处理NetCDF数据的软件也层出不穷，既有开源、免费的，也有需要商业授权的，NetCDF项目网站上列出了许多可用的软件包：

NetCDF_Software_list

当然，我们关心的是Python下处理NetCDF的库，常用的有以下三个：

接下来，我们将用由NOAA提供的全球2018年日最高气温的NetCDF文件（https://www.esrl.noaa.gov/psd/repository/entry/show ）作为示例，介绍三个库的使用。

1
sample_file = 'tmax.2018.nc'

netCDF4

netCDF4，顾名思义可知，是一个用于处理NetCDF数据的专用库。netCDF4是一个在GIS中常用的库，我们之前在地理数据可视化的工作中曾经安装过，作为一个基本库，是其他一些库的依赖。netCDF4的功能比较简单，相对于其他的库（例如xarray）没有一些花里花哨的东西，反而提供了一些底层功能的易用性。

1
2
3
import netCDF4
data = netCDF4.Dataset(sample_file)
data
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4_CLASSIC data model, file format HDF5):
    Conventions: CF-1.0
    Source: ftp://ftp.cpc.ncep.noaa.gov/precip/wd52ws/global_temp/
    version: V1.0
    title: CPC GLOBAL TEMP V1.0
    dataset_title: CPC GLOBAL TEMP
    References: https://www.esrl.noaa.gov/psd/data/gridded/data.cpc.globaltemp.html
    history: Updated 2019-01-28 17:21:38
    dimensions(sizes): lat(360), lon(720), time(365)
    variables(dimensions): float32 [4mlat[0m(lat), float32 [4mlon[0m(lon), float64 [4mtime[0m(time), float32 [4mtmax[0m(time,lat,lon)
    groups:

我们可以像访问词典的方式访问其中的数据，并且查看“变量”、“属性”等信息：

1
2
print(data.dimensions)
print(data.variables)
OrderedDict([('lat', <class 'netCDF4._netCDF4.Dimension'>: name = 'lat', size = 360
), ('lon', <class 'netCDF4._netCDF4.Dimension'>: name = 'lon', size = 720
), ('time', <class 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'time', size = 365
)])
OrderedDict([('lat', <class 'netCDF4._netCDF4.Variable'>
float32 lat(lat)
    actual_range: [ 89.75 -89.75]
    long_name: Latitude
    units: degrees_north
    axis: Y
    standard_name: latitude
    coordinate_defines: center
unlimited dimensions:
current shape = (360,)
filling on, default _FillValue of 9.969209968386869e+36 used
), ('lon', <class 'netCDF4._netCDF4.Variable'>
float32 lon(lon)
    long_name: Longitude
    units: degrees_east
    axis: X
    standard_name: longitude
    actual_range: [2.5000e-01 3.5975e+02]
    coordinate_defines: center
unlimited dimensions:
current shape = (720,)
filling on, default _FillValue of 9.969209968386869e+36 used
), ('time', <class 'netCDF4._netCDF4.Variable'>
float64 time(time)
    long_name: Time
    units: hours since 1900-01-01 00:00:00
    axis: T
    standard_name: time
    coordinate_defines: start
    delta_t: 0000-00-01 00:00:00
    avg_period: 0000-00-01 00:00:00
    actual_range: [1034376. 1043112.]
unlimited dimensions: time
current shape = (365,)
filling on, default _FillValue of 9.969209968386869e+36 used
), ('tmax', <class 'netCDF4._netCDF4.Variable'>
float32 tmax(time, lat, lon)
    missing_value: -9.96921e+36
    units: degC
    level_desc: Surface
    statistic: Mean
    parent_stat: Other
    long_name: Daily Maximum Temperature
    cell_methods: time: mean
    valid_range: [-90.  50.]
    avg_period: 0000-00-01 00:00:00
    dataset: CPC Global Temperature
    comments: GTS data and is gridded using the Shepard Algorithm
    max_period: 6z to 6z
    var_desc: Maximum Temperature
    actual_range: [-58.620934  53.65277 ]
unlimited dimensions: time
current shape = (365, 360, 720)
filling on, default _FillValue of 9.969209968386869e+36 used
)])

1
2
3
4
5
6
7
8
9
tmax= data.variables['tmax']
print('Variable:\n', tmax)

# Attributes can be accessed as properties of a variable
units = tmax.units
print("Attribute:\n", units)

# Variables can indexed numpy-style
tmax[0, 1].shape
Variable:
 <class 'netCDF4._netCDF4.Variable'>
float32 tmax(time, lat, lon)
    missing_value: -9.96921e+36
    units: degC
    level_desc: Surface
    statistic: Mean
    parent_stat: Other
    long_name: Daily Maximum Temperature
    cell_methods: time: mean
    valid_range: [-90.  50.]
    avg_period: 0000-00-01 00:00:00
    dataset: CPC Global Temperature
    comments: GTS data and is gridded using the Shepard Algorithm
    max_period: 6z to 6z
    var_desc: Maximum Temperature
    actual_range: [-58.620934  53.65277 ]
unlimited dimensions: time
current shape = (365, 360, 720)
filling on, default _FillValue of 9.969209968386869e+36 used

Attribute:
 degC

(720,)

Iris

与netCDF4采用类似字典的方式访问数据不同，Iris则是基于CF约定，采用类似list的方式实现了数据的访问。Iris可以识别多种文件格式，包括符合CF约定的NetCDF、GRIB和PP等，并且可以通过插件扩展添加其他格式。

由于Iris要求matplotlib<3，与我其他的包不兼容，也不想再弄虚拟环境，所以就不再写代码了。

xarray

xarray是我们要介绍的主角。个人感觉xarray库应该是目前用于处理NetCDF数据最为方便的库了：xarry有两个核心数据结构：DataArray和Dataset，xarray可以看作是对Numpy和pandas的扩展，同时也借鉴了许多pandas的优点，具有良好的易用性。xarray与netCDF4在使用上存在许多相似之处，但是从输出的结果我们可以发现，同样是Dataset，xarray的结果显然有更好的可读性。

xarray_datastructure

1
2
3
4
5
import xarray

# Open a file
data = xarray.open_dataset(sample_file)
data
<xarray.Dataset>
Dimensions:  (lat: 360, lon: 720, time: 365)
Coordinates:
  * lat      (lat) float32 89.75 89.25 88.75 88.25 87.75 87.25 86.75 86.25 ...
  * lon      (lon) float32 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ...
  * time     (time) datetime64[ns] 2018-01-01 2018-01-02 2018-01-03 ...
Data variables:
    tmax     (time, lat, lon) float32 ...
Attributes:
    Conventions:    CF-1.0
    Source:         ftp://ftp.cpc.ncep.noaa.gov/precip/wd52ws/global_temp/
    version:        V1.0
    title:          CPC GLOBAL TEMP V1.0
    dataset_title:  CPC GLOBAL TEMP
    References:     https://www.esrl.noaa.gov/psd/data/gridded/data.cpc.globa...
    history:        Updated 2019-01-28 17:21:38

1
2
3
4
5
6
7
8
9
10
11
# Variables can be assessed either as properties or as a dict
tmax = data['tmax']
tmax = data.tmax
print('Variable:\n', tmax)

# Same for attributes
units = tmax.units
units = tmax.attrs['units']

print()
print("Attribute:\n", units)
Variable:
 <xarray.DataArray 'tmax' (time: 365, lat: 360, lon: 720)>
[94608000 values with dtype=float32]
Coordinates:
  * lat      (lat) float32 89.75 89.25 88.75 88.25 87.75 87.25 86.75 86.25 ...
  * lon      (lon) float32 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ...
  * time     (time) datetime64[ns] 2018-01-01 2018-01-02 2018-01-03 ...
Attributes:
    units:         degC
    level_desc:    Surface
    statistic:     Mean
    parent_stat:   Other
    long_name:     Daily Maximum Temperature
    cell_methods:  time: mean
    valid_range:   [-90.  50.]
    avg_period:    0000-00-01 00:00:00
    dataset:       CPC Global Temperature
    comments:      GTS data and is gridded using the Shepard Algorithm
    max_period:    6z to 6z
    var_desc:      Maximum Temperature
    actual_range:  [-58.620934  53.65277 ]

Attribute:
 degC

1
2
3
4
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['figure.figsize']=(20,9)
mpl.style.use('ggplot')

在xarray中，调用DataArray的plot()方法可以作图，而且根据数据维数的不同会绘制不同类型的图：

1
2
tmax.plot(bins=50)
plt.show()

output_21_0

我们可以看一下北京（39.9042° N, 116.4074° E）附近网格的2018年的日最高气温情况，可以看到该网格2018年的日最高气温的最大值也就是年最高气温出现在6月29日，温度为39.73摄氏度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
import pandas as pd
lon = int(116.4//0.5)
lat = int((90-39.9)//0.5)
bj_tmax = tmax.isel(lon=lon,lat=lat)

year_max_time = bj_tmax.time[ np.argmax(bj_tmax)].values
ts = pd.to_datetime(year_max_time)
year_max_time_str = ts.strftime('%Y-%m-%d')
year_max_tmax = np.round(np.max(bj_tmax).values.tolist(),2)

bj_tmax.plot()
plt.plot(year_max_time,year_max_tmax, color='blue', marker='o')
plt.text(year_max_time,year_max_tmax,(year_max_time_str,year_max_tmax))
plt.show()

output_23_0

1
2
tmax.isel(time=0).plot()
plt.show()

output_24_0

1
2
3
4
5
6
7
8
9
ax={}
for i in range(4):
    ax[i] = plt.subplot(2,2,i+1)

tmax.isel(time=0).plot.pcolormesh(ax=ax[0])
tmax.isel(time=0).plot.imshow(ax=ax[1])
tmax.isel(time=0).plot.contourf(ax=ax[2])
tmax.isel(time=0).plot.contour(ax=ax[3])
plt.show()

output_25_0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import os
import os.path
import imageio

path = 'img/'

img_files = [f for f in os.listdir('img/')
             if os.path.splitext(f)[1] == ".png"]
img_files.sort(key=lambda f: int(f[:-4]))
img_files = [os.path.join(path, f) for f in img_files]
frames=[]
for image in img_files:
    frames.append(imageio.imread(image))

# save as gif
imageio.mimsave('tmax_2018.gif', frames,
                'GIF', duration=0.5, loop=3)

tmax_2018

此外，我们也可以用cartopy作出更加精美的图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import cartopy.crs as ccrs

# Define the map projection
ax = [
    plt.subplot(121, projection=ccrs.Orthographic(central_longitude=90, central_latitude = 20)),
    plt.subplot(122, projection=ccrs.LambertConformal(central_longitude=115, central_latitude=35)),
]

# ccrs.PlateCarree() is the projection of data
ax[0].set_title('Orthographic')
tmax.isel(time=179).plot.pcolormesh(ax=ax[0], transform=ccrs.PlateCarree(), add_colorbar=False, add_labels=False)
ax[0].coastlines()

ax[1].set_title('LambertConformal')
ax[1].set_extent([60,140, 0, 70])
tmax.isel(time=360).plot.pcolormesh(ax=ax[1], transform=ccrs.PlateCarree(), add_colorbar=False, add_labels=False)
ax[1].coastlines()
plt.show()

output_29_0

小结

本文简单介绍了地理、气象科学中常用的NetCDF数据格式，以及如何使用Python尤其是xarray库处理NetCDF数据。当然，我们目前做的只是非常初步的工作，接下来或许可以尝试以其他气象指标（例如厄尔尼诺？）等例进行分析。