空间管理您的位置: 51Testing软件测试网 » liqianqian1116的个人空间 » 日志

分析气象数据，向Python Cartopy地图中添加循环点

上一篇 / 下一篇 2023-05-04 10:48:45

查看( 121 ) / 评论( 0 ) / 评分( 0 / 0 )

　　在使用Python和Cartopy对气候数据进行可视化分析的过程中，有一个叫做循环点（cyclic point）的术语，它在地理空间栅格数据可视化领域中很重要。

　　什么是循环点，它有什么用？

　　有时，当我们试图绘制地理空间数据时，我们可能会在投影边缘遇到不连续（跳跃）。这通常是由于我们在地理空间数据集中存储经度数据的方式导致的。

　　假设我们有一组具有经度和纬度维度的数据。我们的经度从0°到359.9°（含），分辨率为0.1°（3600个值），纬度从-90°到90°（含），分辨率为0.1°（1801个值）。我们最终的数组形状为(1801, 3600)。

　　为什么我们的数据集不包含360度经度的值？答案很简单，因为360°经度与 0°经度相同。

　　现在，让我们尝试在cartopy地图上绘制数据。我们使用了可以从Climate Data Store下载的数据集，这是一个ERA5-Land从1950年到现在的月平均数据数据集。我们将使用带有cfgrib后端的xarray以grib格式加载此数据集，然后使用matplotlib可视化数据。

　　import xarray as xr

　　import cartopy.crs as ccrs

　　import matplotlib.pyplot as plt

　　import numpy as np

　　# 读取形状为（16, 1801, 3600）的数据集

　　# 选择第5个时间索引，以获得2D数组

　　# 将开尔文转换为摄氏度

　　ds = xr.open_dataset("data/era5_temperature.grib", engine="cfgrib")

　　ds = ds.isel(time=5) - 273.15

　　# 定义绘图函数

　　def plot_data(ds):

　　 # 指定投影

　　 proj = ccrs.PlateCarree()

　　 # 用投影对象创建一个图和轴

　　 fig, ax = plt.subplots(figsize=(16,9), subplot_kw={'projection': proj})

　　 # 为等高线图创建层次

　　 temperature_levels = np.arange(-40,41,2)

　　 # 增加更多的参数

　　 additional_kwargs = {

　　 "antialised": True,

　　 "transform": proj,

　　 "levels": temperature_levels

　　 }

　　 # 绘制数据

　　 ds["t2m"].plot.contourf(ax=ax, cmap='nipy_spectral', **additional_kwargs)

　　 # 设置范围

　　 ax.set_extent([-10, 20, 40, 60], crs=proj)

　　 plt.show()

　　plot_data(ds)

　　请注意从顶部跨越到底部的白线恰好在0°经度处，这就是我们要讨论的问题。

　　向数据添加循环点

　　幸运的是，Cartopy有一个功能可以让我们解决这个问题。该函数称为add_cyclic_point，可以从cartopy.util模块导入。

　　此函数采用数据数组和可选坐标，并返回循环数据和循环坐标。这意味着我们将增加一个等于360的经度坐标，数据等于经度359.9°。让我们用数据本身来证明这一点。

　　首先，我们将循环点添加到我们的数据中，然后比较经度359.9°和经度 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):

　　 # 加载数据

　　 data = ds["t2m"]

　　 # 用循环点生成数据，用循环点生成经度

　　 cyclic_data, cyclic_longitude = add_cyclic_point(data.values, coord=data['longitude'])

　　 # 创建新的坐标，将用于创建新的数据集

　　 # 复制现有数据集的坐标，用循环经度替换经度

　　 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()函数并将其命名为plot_data(new_ds)。当我们绘制使用循环点扩展的数据时，我们可以看到白线消失了，我们的绘图现在是正确的。

　　本文介绍了cartopy中不常见的函数add_cyclic_point。在绘制地理空间数据时，你可能会遇到这样的白线，希望本文能帮助你解决问题。

分析气象数据，向Python Cartopy地图中添加循环点

相关阅读:

用户菜单

标题搜索

日历

我的存档

数据统计

RSS订阅