从xarray走向netCDF处理(四)：合并与计算
以下文章来源于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 删除
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
