从xarray走向netCDF处理(四)：合并与计算

自学气象人

发布于 2022-11-02 10:22:17

1.9K0

文章被收录于专栏：自学气象人自学气象人

以下文章来源于MeteoAI ，作者学前班大队长

前面有关xarray已经讲了3期了，介绍了数据索引，数据结构还有插值和掩膜。今天这是最后一期介绍用xarray处理nc数据了，打算聊一下如何做数据合并与计算。

数据合并

数据合并主要是两种形式

维度的拼接：如将日数据合成为年数据，就属于在时间维度上的合并。
变量的合并：如将多个物理量合到同一个Dataset中。

xarray围绕着这两种合并方式介绍了concatenate, merge, combine, update四种方法。

我在这里就挑最常用的跟大家聊聊。

维度拼接

使用 concat() 方法可以实现维度的拼接。

下面是演示数据，来源于2018年和2019年前三个月的ERA-Interim月平均数据。

>>> ds2018
<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

>>> ds2019
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 3)
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] 2019-01-01 2019-02-01 2019-03-01
Data variables:
    u10        (time, latitude, longitude) float32 ...
    v10        (time, latitude, longitude) float32 ...
    t2m        (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6

ds2018时间维度为12，ds2019时间维度为3，下面使用 concat() 合并后时间维度为15

>>> xr.concat([ds2018, ds2019], dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 15)
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 ... 2019-03-01
Data variables:
    u10        (time, latitude, longitude) float32 -0.9599868 ... 4.5229325
    v10        (time, latitude, longitude) float32 3.1737509 ... -2.289166
    t2m        (time, latitude, longitude) float32 248.46857 ... 225.19632
Attributes:
    Conventions:  CF-1.6

变量合并

使用 merge() 方法，可以将ds2018中的u10和ds2019中的t2m合并到一起，而且在时间维上缺失会自动设置为nan。

>>> xr.merge([ds2018.u10, ds2019.t2m])
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 15)
Coordinates:
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2019-03-01
  * 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
Data variables:
    u10        (time, latitude, longitude) float32 -0.9599868 -0.9599868 ... nan
    t2m        (time, latitude, longitude) float32 nan nan ... 225.19632

数据计算

最基本的计算就是进行加减乘除，任意一个DataArray或者Dataset都可以直接进行四则运算。

除此以外，xarray还可以帮你快速地求出平均值，方差，最小值，最大值等。你可以指定具体对那个维度进行计算，如果不指定维度默认会对所有维度进行计算。

比如要对经、纬两个维度进行平均，最后的结果只有时间维的12个值。

>>> ds2018.mean(dim=['latitude', 'longitude'])
<xarray.Dataset>
Dimensions:  (time: 12)
Coordinates:
  * time     (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Data variables:
    u10      (time) float32 -0.120867714 0.13738841 ... -0.016953295 0.08418254
    v10      (time) float32 -0.21417202 -0.12471106 ... -0.082666814
    t2m      (time) float32 277.57446 277.32916 276.72095 ... 278.1613 278.01758

如果对时间维进行求方差，则结果会保留空间场上的信息。

>>> ds2018.std(dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480)
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
Data variables:
    u10        (latitude, longitude) float32 1.9192954 1.9192954 ... 1.2133
    v10        (latitude, longitude) float32 1.3066719 1.3066719 ... 1.577495
    t2m        (latitude, longitude) float32 9.5681305 9.5681305 ... 11.313364

而且xarray在时间维上的计算还有很多贴心的用法，比如月数据转年数据，月数据转季节数据。

>>> ds2018.groupby('time.season').min(dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, season: 4)
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
  * season     (season) object 'DJF' 'JJA' 'MAM' 'SON'
Data variables:
    u10        (season, latitude, longitude) float32 -3.2505732 ... 3.4194348
    v10        (season, latitude, longitude) float32 0.9903783 ... -3.318876
    t2m        (season, latitude, longitude) float32 248.46857 ... 214.61685

一个栗子

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

def create_map(ax):
    # 设置地图属性
    ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8)  # 加载分辨率为50的国界线
    ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6)  # 加载分辨率为50的海岸线
    ax.add_feature(cfeat.RIVERS.with_scale('50m'))  # 加载分辨率为50的河流
    ax.add_feature(cfeat.LAKES.with_scale('50m'))  # 加载分辨率为50的湖泊
    # 设置网格点属性
    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轴设为纬度的格式
    return ax

if __name__ == '__main__':
    # 创建画图空间
    proj = ccrs.PlateCarree()
    fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(7, 14), subplot_kw={'projection': proj})
    ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
    temp = (ds['t2m'] - 273.15).groupby('time.season').mean('time')
    # --画图
    cbar_kwargs = {'label': '2m temperature (℃)',
       'ticks': np.arange(-30, 30+5, 5)}
    levels = np.arange(-30, 30+1, 1)
    for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):
        ax = create_map(axes[i])
        temp.sel(season=season).plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', extend='both',
            cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
    fig.show()

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-12-17，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自自学气象人微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度