多維數組庫Xarray 在數據處理中的應用

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

（1、山東省海洋生態環境與防災減災重點實驗室，山東 青島266033 2、國家海洋局北海預報中心，山東 青島 266033）

1 概述

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

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

Xarray 最早開發與2014 年，期間經過數十次更新迭代，功能已經日趨成熟和強大，目前最新版本為v0.19.0，本文內容基于該版本。

2 數據的查看與獲取

針對不同需求，Xarray 提供了兩種數據類型：

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

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

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

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

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

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

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

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

屬性中即顯示變量msl的單位為Pa，全名為“Mean sea level pressure”，標準名稱為“air_pressure_at_mean_sea_level”。

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

msl_data=file1.msl.data

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

3 子集選取與合并

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

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

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

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

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

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

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

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

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

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

4 數據插值與計算

為了方便數據分析，Xarray 還提供了一些較基礎的計算工具，包括時間、空間插值和簡單的統計計算。

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

圖1 氣壓變化折線圖（（a）為默認樣式，（b）為添加顏色和標記樣式參數）

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

5 可視化

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

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

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

圖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()

網格填色圖：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 分別對應著不同的圖表類型，即有的圖表類型以DataArray 為輸入，而另一些則以Dataset 為輸入，這是由圖表性質決定的。上述例子中，折線圖、直方圖、階梯圖、填色等值線圖和網格填色圖的輸入為DataArray，矢量圖和流線圖的輸入為Dataset。

6 數據的修改、創建與保存

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

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

用生成的新數據為目標對象賦值：

最后修改屬性中history 鍵對應值：

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

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

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

7 結論

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

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