分析氣象數(shù)據(jù)，向Python Cartopy地圖中添加循環(huán)點
在使用Python和Cartopy對氣候數(shù)據(jù)進(jìn)行可視化分析的過程中，有一個叫做循環(huán)點（cyclic point）的術(shù)語，它在地理空間柵格數(shù)據(jù)可視化領(lǐng)域中很重要。
什么是循環(huán)點，它有什么用？
有時，當(dāng)我們試圖繪制地理空間數(shù)據(jù)時，我們可能會在投影邊緣遇到不連續(xù)（跳躍）。這通常是由于我們在地理空間數(shù)據(jù)集中存儲經(jīng)度數(shù)據(jù)的方式導(dǎo)致的。
假設(shè)我們有一組具有經(jīng)度和緯度維度的數(shù)據(jù)。我們的經(jīng)度從0°到359.9°（含），分辨率為0.1°（3600個值），緯度從-90°到90°（含），分辨率為0.1°（1801個值）。我們最終的數(shù)組形狀為(1801, 3600)。
為什么我們的數(shù)據(jù)集不包含360度經(jīng)度的值？答案很簡單，因為360°經(jīng)度與 0°經(jīng)度相同。
現(xiàn)在，讓我們嘗試在cartopy地圖上繪制數(shù)據(jù)。我們使用了可以從Climate Data Store下載的數(shù)據(jù)集，這是一個ERA5-Land從1950年到現(xiàn)在的月平均數(shù)據(jù)數(shù)據(jù)集。我們將使用帶有cfgrib后端的xarray以grib格式加載此數(shù)據(jù)集，然后使用matplotlib可視化數(shù)據(jù)。
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
# 讀取形狀為（16, 1801, 3600）的數(shù)據(jù)集
# 選擇第5個時間索引，以獲得2D數(shù)組
# 將開爾文轉(zhuǎn)換為攝氏度
ds = xr.open_dataset("data/era5_temperature.grib", engine="cfgrib")
ds = ds.isel(time=5) - 273.15
# 定義繪圖函數(shù)
def plot_data(ds):
# 指定投影
proj = ccrs.PlateCarree()
# 用投影對象創(chuàng)建一個圖和軸
fig, ax = plt.subplots(figsize=(16,9), subplot_kw={'projection': proj})
# 為等高線圖創(chuàng)建層次
temperature_levels = np.arange(-40,41,2)
# 增加更多的參數(shù)
additional_kwargs = {
"antialised": True,
"transform": proj,
"levels": temperature_levels
}
# 繪制數(shù)據(jù)
ds["t2m"].plot.contourf(ax=ax, cmap='nipy_spectral', **additional_kwargs)
# 設(shè)置范圍
ax.set_extent([-10, 20, 40, 60], crs=proj)
plt.show()
plot_data(ds)
請注意從頂部跨越到底部的白線恰好在0°經(jīng)度處，這就是我們要討論的問題。
向數(shù)據(jù)添加循環(huán)點
幸運的是，Cartopy有一個功能可以讓我們解決這個問題。該函數(shù)稱為add_cyclic_point，可以從cartopy.util模塊導(dǎo)入。
此函數(shù)采用數(shù)據(jù)數(shù)組和可選坐標(biāo)，并返回循環(huán)數(shù)據(jù)和循環(huán)坐標(biāo)。這意味著我們將增加一個等于360的經(jīng)度坐標(biāo)，數(shù)據(jù)等于經(jīng)度359.9°。讓我們用數(shù)據(jù)本身來證明這一點。
首先，我們將循環(huán)點添加到我們的數(shù)據(jù)中，然后比較經(jīng)度359.9°和經(jīng)度 360°處的值。
from cartopy.util import add_cyclic_point
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
def add_cyclic_point_to_dataset(ds):
# 加載數(shù)據(jù)
data = ds["t2m"]
# 用循環(huán)點生成數(shù)據(jù)，用循環(huán)點生成經(jīng)度
cyclic_data, cyclic_longitude = add_cyclic_point(data.values, coord=data['longitude'])
# 創(chuàng)建新的坐標(biāo)，將用于創(chuàng)建新的數(shù)據(jù)集
# 復(fù)制現(xiàn)有數(shù)據(jù)集的坐標(biāo)，用循環(huán)經(jīng)度替換經(jīng)度
coords = {dim: data.coords[dim] for dim in data.dims}
coords["longitude"] = cyclic_longitude
new_ds = xr.Dataset(
data_vars={
"t2m": xr.DataArray(cyclic_data, dims=data.dims, coords=coords)
})
return new_ds
new_ds = add_cyclic_point_to_dataset(ds)
new_ds["t2m"].shape
# (1801, 3601)
difference = ds["t2m"].sel(lnotallow=0, method="nearest") - new_ds["t2m"].sel(lnotallow=360, method="nearest")
difference = difference.fillna(0) # fill nan values with 0
any(difference != 0)
# False讓我們重用plot_data()函數(shù)并將其命名為plot_data(new_ds)。當(dāng)我們繪制使用循環(huán)點擴(kuò)展的數(shù)據(jù)時，我們可以看到白線消失了，我們的繪圖現(xiàn)在是正確的。
本文介紹了cartopy中不常見的函數(shù)add_cyclic_point。在繪制地理空間數(shù)據(jù)時，你可能會遇到這樣的白線，希望本文能幫助你解決問題。
