多維數(shù)組庫Xarray 在數(shù)據(jù)處理中的應(yīng)用

曹永正 王立鵬 高 山 溫連杰 王 翠

（1、山東省海洋生態(tài)環(huán)境與防災(zāi)減災(zāi)重點實驗室，山東 青島266033 2、國家海洋局北海預(yù)報中心，山東 青島 266033）

1 概述

地球科學(xué)領(lǐng)域一直與“大數(shù)據(jù)”密切相關(guān)，高分辨率的觀測或模擬資料能夠更準(zhǔn)確地刻畫不同尺度的物理過程，但也給數(shù)據(jù)分析處理帶來了困難。近年來隨著地球觀測手段和物理模型的進步，人類能夠獲取的氣象、海洋和陸地要素越來越多，其中再分析和模式預(yù)報結(jié)果多以正交網(wǎng)格形式呈現(xiàn)，而船舶、浮標(biāo)和測站觀測記錄則是散點形式。為了更方便處理這些不同要素、不同形式的數(shù)據(jù)，Xarray 提供了一套較為統(tǒng)一的格式和操作工具。

Xarray 是一個開源項目，面向需要N 維數(shù)組，尤其是需要數(shù)據(jù)分析的科研工作，旨在讓多維數(shù)據(jù)有關(guān)的工作更加簡單、高效。Xarray 在類似NumPy的原始數(shù)組之上引入了維度、坐標(biāo)和屬性形式的標(biāo)簽，氣象、海洋數(shù)據(jù)的地理空間和時間、單位、變量名稱、注釋等信息在處理過程中都得以保留。Xarray 通過數(shù)據(jù)引擎的方式，對netCDF、GRIB 和ZARR 等研究中常用格式都提供了相應(yīng)支持。

Xarray 最早開發(fā)與2014 年，期間經(jīng)過數(shù)十次更新迭代，功能已經(jīng)日趨成熟和強大，目前最新版本為v0.19.0，本文內(nèi)容基于該版本。

2 數(shù)據(jù)的查看與獲取

針對不同需求，Xarray 提供了兩種數(shù)據(jù)類型：

DataArray：帶有標(biāo)簽或名稱的多維數(shù)組。DataArray的基礎(chǔ)是一個多維數(shù)組矩陣，在此基礎(chǔ)上添加了用以描述該數(shù)組的元數(shù)據(jù)信息，既維度名稱、坐標(biāo)和屬性等。

Dataset：可以看做是多個DataArray的集合，這些DataArray具有相同的維度，是彼此“對齊”的，他們有各自的元數(shù)據(jù)信息，但也可以有公用的部分。大多數(shù)可以在單個DataArray 上執(zhí)行的操作都可以在Dataset 上執(zhí)行。

對于只包含單個要素的數(shù)據(jù)文件，構(gòu)建DataArray 類型對象即可，但實際工作中為了方便分發(fā)，一個數(shù)據(jù)文件中往往包含相同時空范圍的多個要素信息，因此Dataset 類型更常使用。

不管是DataArray 還是Dataset，其中的數(shù)組都以Variable 對象的形式呈現(xiàn)，各種運算可以通過維度名稱實現(xiàn)數(shù)組廣播。

我們以包含2020 年12 月一個月內(nèi)10m 風(fēng)速u、v 分量和海平面氣壓3 個要素的ERA5 再分析數(shù)據(jù)（文件名202012_atm.nc）為例，首先在Python 環(huán)境下導(dǎo)入Xarray 庫和Python 環(huán)境下的高效矩陣計算庫numpy、數(shù)據(jù)分析庫pandas：

可以看到，Xarray 以202012_atm.nc 數(shù)據(jù)文件為基礎(chǔ)，創(chuàng)建了一個Dataset 對象。對象有3 個維度（Dimensions），分別為緯度（latitude）、經(jīng)度（longitude）和時間（time），各維度方向的格點數(shù)為721、1440 和744。

所有數(shù)值數(shù)組均為Variable 對象格式，其中各維度數(shù)值歸類到“坐標(biāo)（Coordinates）”中，而變量10 米經(jīng)向風(fēng)矢量（u10）、10米緯向風(fēng)矢量（v10）和平均海平面氣壓（msl）則歸屬于“數(shù)據(jù)變量（Data variables）”。其中緯度、經(jīng)度和時間為1 維變量，即沿自身方向延伸的1 維數(shù)組，而u10、v10 和msl 為3 維變量，第0 維為時間，第1 維為緯度，第2 維為經(jīng)度。除時間外，各項數(shù)值以32bits 浮點數(shù)格式存儲，而時間則轉(zhuǎn)換為專有的datetime64 格式以方便計算。

本數(shù)據(jù)集中含有2 項共有屬性（Attributes），Conventions 和history，描述了數(shù)據(jù)來源和生成操作歷史。要想進一步查看每個變量的屬性信息，可以在環(huán)境中輸入“file1.msl”，顯示如下：

屬性中即顯示變量msl的單位為Pa，全名為“Mean sea level pressure”，標(biāo)準(zhǔn)名稱為“air_pressure_at_mean_sea_level”。

為了進行復(fù)雜的分析計算，我們可以把數(shù)據(jù)從Dataset 或DataArray 中提取為numpy 支持的無屬性數(shù)組格式：

msl_data=file1.msl.data

提取出的數(shù)組不再含有坐標(biāo)和屬性信息，為單純的3 維格點氣壓值。

3 子集選取與合并

知道數(shù)據(jù)集基本結(jié)構(gòu)后，我們就可以根據(jù)需要提取相應(yīng)子集。對于DataArray 對象，Xarray 提供的loc()方法可以讓用戶像對待numpy 中Array 對象一樣進行數(shù)據(jù)切片操作。例如從全球范圍分辨率0.25°的原始氣壓數(shù)據(jù)中選取中國渤、黃海范圍1°分辨率氣壓場，通過以下命令就能得到：

target_field=file1.msl.loc[:,40:30:4,117:127:4]

但該方法無法用于Dataset 對象。對于Dataset，需要借助sel()或isel()方法。其中sel()方法的輸入?yún)?shù)為坐標(biāo)值或坐標(biāo)值列表，isel()方法的輸入?yún)?shù)為坐標(biāo)位置索引。兩種方法用法相似，但sel()對大多數(shù)使用者來說更直觀。例如選取u10、v10 和msl全部3 個變量的渤、黃海范圍1°分辨率氣壓場，可以：

除了提取子集，在實際工作中也經(jīng)常會遇到需要將兩個或多個較小數(shù)據(jù)集合并為一個大數(shù)據(jù)集，Xarray 中可以用concat()和merge()方法實現(xiàn)。這里增加2021 年1 月數(shù)據(jù)集：

file2=xr.open_dataset('202101_atm.nc')

concat()方法是維度上的合并，需要指定操作維度，需要合并的數(shù)據(jù)集以列表方式傳入：

target_file=xr.concat([file1,file2],dim="time")

而merge()方法是變量上的合并，需要合并的數(shù)據(jù)集同樣以列表方式傳入：

target_file=xr.merge([file1,file2])

merge()方法不需要指定操作維度，倘若file1與file2 坐標(biāo)上存在差異，其維度坐標(biāo)會自動補全，但變量場的缺失值會自動設(shè)置為缺省值Nan。

4 數(shù)據(jù)插值與計算

為了方便數(shù)據(jù)分析，Xarray 還提供了一些較基礎(chǔ)的計算工具，包括時間、空間插值和簡單的統(tǒng)計計算。

由于Xarray 將時間、空間維度都看做不同名稱命名的坐標(biāo)，因此各維度的插值均使用interp()方法實現(xiàn)。使用時，僅需將要調(diào)整的維度坐標(biāo)值傳入，未傳入的坐標(biāo)將保持原樣。例如要插值得到50-55°N、120-125°E 范圍內(nèi)分辨率0.1°，2020 年12 月15 日08-14 時期間間隔0.5 小時的氣壓場，可以通過以下命令：latitude=lat_new,time=time_new)

圖1 氣壓變化折線圖（（a）為默認樣式，（b）為添加顏色和標(biāo)記樣式參數(shù)）

Xarray 內(nèi)置了幾種基礎(chǔ)的統(tǒng)計方法，有最大值（max()）、最小值（min()）、平均值（mean()）、中值（median()）、總和（sum()）和標(biāo)準(zhǔn)差（std()）等。但來自外部的數(shù)學(xué)運算多數(shù)也可以直接進行，例如四則運算、三角函數(shù)計算等。這些運算會直接作用在“Data variables”下的每個Variable 對象中，相當(dāng)于對多維矩陣進行的計算。而最終結(jié)果會還原輸入形式，即DataArray 或Dataset。

5 可視化

Matplotlib 是Python 環(huán)境下功能最全面、應(yīng)用最廣泛的可視化工具庫，借助Matplotlib 強大的繪圖功能，Xarray 可以對DataArray 和Dataset 數(shù)據(jù)進行多種形式的可視化。由于集成了Matplotlib 應(yīng)用接口，使用Xarray 時不必先將數(shù)據(jù)導(dǎo)出到Matplotlib，而是直接調(diào)用plot()函數(shù)。這樣不僅避免了繁瑣的數(shù)據(jù)轉(zhuǎn)移流程，還可以方便地將DataArray 或Dataset 中的屬性和維度等信息應(yīng)用在繪圖中。

目前Xarray 支持的繪圖樣式涵蓋了一維、二維和三維形式下多種常見的圖表。plot()函數(shù)默認情況下為折線圖，例如，查看2020 年12 月期間（120°E，36°N）處氣壓變化情況，只需直接使用plot()命令即可得到相應(yīng)圖像，通過添加繪圖參數(shù)可進一步控制圖表樣式：

得到的圖像中，Xarray 自動為橫軸采用了時間格式標(biāo)簽，為縱軸標(biāo)注了變量名稱和單位，并將經(jīng)緯度坐標(biāo)作為圖像標(biāo)題。

圖2 Xarray 支持的常用圖表樣式

除最基本的折線圖外，Xarray 還可以繪制直方圖、階梯圖、二維等值線圖、矢量圖、流線圖等，示例如下：

直方圖：file1.sel(longitude=120,latitude=36).msl.plot.hist()

階梯圖：file1.sel(longitude=120,latitude=36).msl[:10].plot.step()

填色等值線圖：file1.sel(longitude=np.linspace(120,130,11),latitude=np.linspace(30,40,11),time='2020-12-01T12').msl.plot.contour()

網(wǎng)格填色圖：file1.sel (longitude=np.linspace (120,130,11),latitude=np.linspace(30,40,11),time='2020-12-01T12').msl.plot.pcolormesh()

矢量圖：file1.sel (longitude=np.linspace (120,130,11),latitude=np.linspace (30,40,11),time='2020-12-01T12').plot.quiver(x='longitude',y='latitude',u='u10',v='v10')

流線圖：file1.sel (longitude=np.linspace (120,130,11),latitude=np.linspace (30,40,11),time='2020-12-01T12').plot.streamplot(x='longitude',y='latitude',u='u10',v='v10')

需要注意的是，DataArray 和Dataset 分別對應(yīng)著不同的圖表類型，即有的圖表類型以DataArray 為輸入，而另一些則以Dataset 為輸入，這是由圖表性質(zhì)決定的。上述例子中，折線圖、直方圖、階梯圖、填色等值線圖和網(wǎng)格填色圖的輸入為DataArray，矢量圖和流線圖的輸入為Dataset。

6 數(shù)據(jù)的修改、創(chuàng)建與保存

對于讀取為DataArray 或Dataset 格式的數(shù)據(jù)，可以使用Xarray 對其數(shù)值和維度、屬性等信息進行修改。由于這些信息在DataArray 或Dataset 中以對象的形式呈現(xiàn)，因此修改實質(zhì)上是對目標(biāo)對象的操作。Xarray 為不同對象提供了不同操作方式，其中變量數(shù)值可以直接用賦值的方式修改，維度需要通過專門函數(shù)assign_coords()，屬性則為“鍵-值對”的字典形式。例如，將file1 中的時間維由2020 年12 月改為1 月，變量msl的數(shù)值用隨機數(shù)替代，并將屬性中的“history”內(nèi)容改為其它文本，可以通過以下幾個步驟實現(xiàn)。首先借助pandas 生成新時間軸序列：

使用numpy 生成與變量msl 具有相同維度范圍的隨機數(shù)矩陣：

用生成的新數(shù)據(jù)為目標(biāo)對象賦值：

最后修改屬性中history 鍵對應(yīng)值：

file1.attrs['history']='Edited by XXX at 2021-9-8'

在已有其它來源數(shù)據(jù)的情況下，還可以使用Xarray 從頭創(chuàng)建一個DataArray 或Dataset，進而保存為netCDF 等文件格式。這里我們使用剛才生成的隨機數(shù)矩陣,可以通過如下方式構(gòu)建一個新DataArray 并以netCDF 格式保存：

其中da_new.nc 為此次保存的文件名。修改過的file1 也可以用同樣的方法保存在本地。

7 結(jié)論

本文探索和介紹了Xarray 在地球物理數(shù)據(jù)處理中的常用功能和使用方法，對于科研中常見的資料格式，例如netCDF、GRIB 和ZARR 等，Xarray 提供了一套強大的工具和命令，除了本文中介紹的數(shù)據(jù)讀取、查看、選取與合并、插值計算、可視化、修改保存等功能外，還可以借助Lamda 表達式實現(xiàn)自定義函數(shù)，以及與Numpy、Pandas、Dask、matplotlib 等工具結(jié)合使用，實現(xiàn)更復(fù)雜的數(shù)據(jù)分析、可視化，甚至并行計算。在數(shù)據(jù)量日益增大的今天，Python 環(huán)境下的數(shù)據(jù)科學(xué)生態(tài)系統(tǒng)成為了不可或缺的生產(chǎn)力來源，而Xarray 就為地球科學(xué)領(lǐng)域研究者提供了進入這套生態(tài)系統(tǒng)最方便的門戶和工具。

Xarray 具有使用簡單、集成度高、可拓展性強的優(yōu)點，但也存在著如對復(fù)雜格式文件支持不完善、部分函數(shù)的計算效率低、與其它工具同時使用偶現(xiàn)Bug 等，目前仍在開發(fā)中。當(dāng)前版本Xarray 功能已經(jīng)覆蓋了大部分數(shù)據(jù)操作需求，能夠為數(shù)據(jù)研究分析工作帶來很大方便。