使用Cartopy和Xarray可视化气候数据，使用Python进行气候数据可视化。
微信搜索关注《Python学研大本营》，加入读者群，分享更多精彩
尽管现在有很多BI工具，但Python仍然是非常优秀的数据可视化工具。它的简单性允许用户以有意义的方式快速阅读大量类型的数据。
Python在读取非结构化数据方面非常高效，例如有关气候的数据。大多数此类数据通常以grib/grib2格式（广泛用于气象学）或NC文件（用于存储多维数据的NetCDF格式）分发。
由于气候变化一直在敦促许多大公司和政府收集有关我们不断变化的星球的数据，现在有数PB的数据可以免费下载并用于良好目的和常识。
只需几行代码和很短的时间，你就可以创建有意义的图表。使用xarray和cartopy等库，我们可以快速选择我们想要的数据、裁剪数据、转换数据、在特定坐标系中绘制地图上的数据等等。
使用xarray读取数据
存在许多库来探索气候数据，如GDAL、rasterio、pygrib、cfgrib等等。这些库要么安装起来太复杂，要么包含导致意外错误的交叉依赖项。这并不是说它们不好。有些人正在使用这些库，它们可能适合特定用户。然而，对于气候和气象数据，xarray是最合适和最通用的库。同样重要的是要提到 cfgrib库用作后端引擎来读取xarray中的grib文件。
简单介绍一下xarray。我们将首先探索数据集，然后使用xarray及其绘图方法创建简单的绘图。
探索数据集
对于此分析，我们将使用我们在第一篇文章中探讨过的grib文件。它将存储在data/1month_anomaly_Global_ea_2t_201907_1991–2020_v02.grib下
该数据集包含全球温度异常，取2019年7月的月度平均值。该异常计算为温度与1991年至2020年长期平均值的偏差。
import xarray as xr
data_path = "data/1month_anomaly_Global_ea_2t_201907_1991-2020_v02.grib"
original_data = xr.open_dataset(data_path, engine="cfgrib")
original_data
从上面的代码返回的xarray的数据集类的结构
打开grib文件后，我们可以看到Dimensions、Coordinates、Data variables和Attributes。这些元素中的每一个都包含有关我们数据的有用信息。如我们所见，我们有两个维度——纬度和经度。这些维度也是我们数据的坐标——这意味着我们数据中的每个点都必须具有唯一的纬度和经度组合。坐标time、step、heightAboveGround和valid_time只有一个值——它是对我们底层DataArray中的所有元素都有效的单一信息。
对我们来说最重要的数据存储在数据变量中，在本例中为键t2m。这是我们的底层 DataArray，其中包含有关温度的所有信息。这个DataArray必须具有MxN的形状（在我们的例子中，M是经度，N是纬度）—— t2m DataArray的每个元素都包含纬度和经度的唯一组合。
创建简单的图
为了可视化DataArray，我们可以简单地键入data[“t2m”].plot()。请注意，数据以开尔文存储。然而，由于这个数据集包含异常，我们将其存储为开尔文或摄氏度并不重要——高于平均水平1摄氏度与高于平均水平1开尔文完全相同。
data["t2m"].plot()
正如我们所见，这是我们数据集的非常简单的图。虽然我们可以看到在经度和纬度轴上正确绘制的温度异常，但它缺乏更多的上下文。使用cartopy，我们可以添加网格线、海岸线、边界、河流等。此外，我们可以根据需要改变我们的坐标参考系统来转换数据。
使用Cartopy创建地图
使用cartopy，我们可以轻松地在不同的坐标参考系之间切换，改变我们地图的投影。我们可以添加网格线、海岸线、国家边界、河流、道路等等。首先，重要的是要了解坐标参考系统（CRS = 地理坐标系 — GCS）和投影参考系统 (PRS) 之间的区别。在我们的代码中，我们将使用projection来谈论PRS和crs来谈论CRS。
我们的数据以度数为单位存储纬度和经度。它也可以以不同的单位存储（例如，在 x 和y方向上距地图上特定点的距离（以米为单位），但这不是我们的情况。由于我们数据的crs以度为单位指定，因此我们将使用名为PlateCarree的cartopy的CRS。这就是我们确保在所需位置绘制数据的方式。
长话短说，我们将使用projection来指定地图的投影，并使用crs来告诉我们的代码应该在地图上的哪个位置绘制数据。让我们跳回代码。
创建底图
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
# 首先，我们指定地图投影的坐标参考系。
# 使用Mercator，它是一个圆柱形的保角投影。
# 高纬度地区有非常大的失真，不能完全到达极地地区。
projection = ccrs.Mercator()
# 指定CRS，它将被用来告诉代码，我们的数据应该被绘制在哪里。
crs = ccrs.PlateCarree()
# 现在我们将创建具有特定投影的轴对象。
plt.figure(dpi=150)
ax = plt.axes(projection=projection, frameon=True)
# 在Mercator地图上以度绘制网格线。
gl = ax.gridlines(crs=crs, draw_labels=True,
linewidth=.6, color='gray', alpha=0.5, linestyle='-.')
gl.xlabel_style = {"size" : 7}
gl.ylabel_style = {"size" : 7}
# 为了绘制边界和海岸线，我们可以使用cartopy功能。
import cartopy.feature as cf
ax.add_feature(cf.COASTLINE.with_scale("50m"), lw=0.5)
ax.add_feature(cf.BORDERS.with_scale("50m"), lw=0.3)
# 现在，我们将以最小/最大经度/纬度来指定我们地图的范围。注意，这些值是以经度和纬度来指定的。然而，我们可以用任何我们想要的经度来指定它们，但我们需要在ax.set_extent中提供适当的经度参数。
lon_min = -20
lon_max = 45
lat_min = 34
lat_max = 60
# crs是PlateCarree -- 我们明确地告诉轴，我们正在创建以度为单位的边界。
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=crs)
plt.title(f"Temperature anomaly over Europe in {original_data.valid_time.dt.strftime('%B %Y').values}")
plt.show()
这会生成我们的底图，其中包含墨卡托投影、表示纬度和经度的网格线、边界和海岸线，同时我们的范围设置为欧洲。您可以指定您想要的任何范围。
在底图上绘制数据
现在让我们在底图上绘制数据。有人会认为它需要一些额外的数据操作，但事实并非如此。要在我们的底图上绘制数据，我们只需添加 1 行代码。我将根据需要添加更多行来绘制颜色条。我们的地图代码如下所示：
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
def plot_dataset(dataset : xr.Dataset):
# 首先，我们指定地图投影的坐标参考系。
# 使用Mercator，它是一个圆柱形的保角投影。
projection = ccrs.Mercator()
# 指定CRS，它将被用来告诉代码，我们的数据应该被绘制在哪里。
crs = ccrs.PlateCarree()
# 现在我们将创建具有特定投影的轴对象。
plt.figure(figsize=(16,9), dpi=150)
ax = plt.axes(projection=projection, frameon=True)
# 在使用Mercator地图上以度绘制网格线。
gl = ax.gridlines(crs=crs, draw_labels=True,
linewidth=.6, color='gray', alpha=0.5, linestyle='-.')
gl.xlabel_style = {"size" : 7}
gl.ylabel_style = {"size" : 7}
# 为了绘制边界和海岸线，我们可以使用cartopy功能。
import cartopy.feature as cf
ax.add_feature(cf.COASTLINE.with_scale("50m"), lw=0.5)
ax.add_feature(cf.BORDERS.with_scale("50m"), lw=0.3)
# 现在，我们将以最小/最大经度/纬度指定我们地图的范围。
# 请注意，这些值是以经度和纬度为单位的。
# 然而，我们可以在任何我们想要的crs中指定它们，但我们需要提供适当的ax.set_extent中的crs参数。
# crs是PlateCarree -- 我们明确地告诉轴，我们正在创建以度为单位的边界。
lon_min = -20
lon_max = 45
lat_min = 34
lat_max = 60
##### WE ADDED THESE LINES #####
cbar_kwargs = {'orientation':'horizontal', 'shrink':0.6, "pad" : .05, 'aspect':40, 'label':'2 Metre Temperature Anomaly [K]'}
dataset["t2m"].plot.contourf(ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs=cbar_kwargs, levels=21)
################################
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=crs)
plt.title(f"Temperature anomaly over Europe in {dataset.valid_time.dt.strftime('%B %Y').values}")
plt.show()
plot_dataset(original_data)
请注意，绘制可能需要相当长的时间（此处用了49秒）。这是因为我们正在绘制整个数据集，无论我们指定的范围如何。在我们绘制数据后应用范围。这也可以通过自动选择的色标来证明。它的范围远远超出了我们在地图上看到的范围。
切片数据
为了加快这个过程，最佳做法是将数据切片到我们的特定范围。它将确保我们只绘制我们范围内的数据。这可以简单地使用 xarray 的sel和slice来完成：
from copy import deepcopy
# 创建数据的副本，以便我们保持原始数据。
cropped_dataset = deepcopy(original_data)
cropped_dataset = cropped_dataset.sel(
latitude=slice(lat_max, lat_min),
longitude=slice(lon_min, lon_max)
)
plot_dataset(cropped_dataset)
这样做，我们从49s变成了5s。这有很大的不同！
修复白线
如果仔细查看地图，您会发现问题所在——0 度经度上的白线。这是由我们存储经度的方式引起的。在此数据集中，我们使用0到360度约定。这意味着我们的数据集从经度0度开始，经度360度结束。我们有2个选项来改变这个：
-
添加循环点来环绕我们的数据集，
-
将经度转换为-180到180约定（请注意，这将在地球的另一侧创建白线 - 在180度经度上，但它不会在我们的地图上可见，因为它不在我们指定的范围内）。最好的选择是同时使用它们。为了本文的简单起见，我们使用了第二个选项，它不需要太多的编码和对循环点的深入理解。
请注意，我们首先创建名为_longitude_adjusted的虚拟维度。我们使用 xr.where() 方法——我们从每个大于 180 的经度中减去 360。然后我们将经度与_longitude_adjusted交换并删除经度。最后，我们将_longitude_adjusted重命名为经度。
def adjust_longitude(dataset: xr.Dataset) -> xr.Dataset:
"""Swaps longitude coordinates from range (0, 360) to (-180, 180)
Args:
dataset (xr.Dataset): xarray Dataset
Returns:
xr.Dataset: xarray Dataset with swapped longitude dimensions
"""
lon_name = "longitude"
# 调整lon值，以确保它们在（-180，180）范围内。
dataset["_longitude_adjusted"] = xr.where(
dataset[lon_name] > 180, dataset[lon_name] - 360, dataset[lon_name]
)
dataset = (
dataset.swap_dims({lon_name: "_longitude_adjusted"})
.sel(**{"_longitude_adjusted": sorted(dataset._longitude_adjusted)})
.drop(lon_name)
)
dataset = dataset.rename({"_longitude_adjusted": lon_name})
return dataset
整合起来
总结一下，这是我们的最终代码，用于分割数据并从我们的图中删除白线。
from copy import deepcopy
from cartopy.util import add_cyclic_point
# 创建数据的副本，以便我们保持原始数据。
cropped_dataset = deepcopy(original_data)
def adjust_longitude(dataset: xr.Dataset) -> xr.Dataset:
"""Swaps longitude coordinates from range (0, 360) to (-180, 180)
Args:
dataset (xr.Dataset): xarray Dataset
Returns:
xr.Dataset: xarray Dataset with swapped longitude dimensions
"""
lon_name = "longitude" # whatever name is in the data
# 调整lon值，以确保它们在（-180，180）范围内。
dataset["_longitude_adjusted"] = xr.where(
dataset[lon_name] > 180, dataset[lon_name] - 360, dataset[lon_name]
)
dataset = (
dataset.swap_dims({lon_name: "_longitude_adjusted"})
.sel(**{"_longitude_adjusted": sorted(dataset._longitude_adjusted)})
.drop(lon_name)
)
dataset = dataset.rename({"_longitude_adjusted": lon_name})
return dataset
cropped_dataset = adjust_longitude(cropped_dataset)
cropped_dataset = cropped_dataset.sel(
latitude=slice(lat_max, lat_min),
longitude=slice(lon_min, lon_max)
)
plot_dataset(cropped_dataset)
结论
这篇文章向我们展示了结合使用xarray和cartopy是多么容易。使用这两个库，数据分析师可以更快、更高效地可视化气候数据。有太多免费数据可供探索，尤其是在Climate Data Store（https://cds.climate.copernicus.eu/cdsapp）上。
推荐书单
《Python从入门到精通（第2版）》
《Python从入门到精通（第2版）》从初学者角度出发，通过通俗易懂的语言、丰富多彩的实例，详细介绍了使用Python进行程序开发应该掌握的各方面技术。全书共分23章，包括初识Python、Python语言基础、运算符与表达式、流程控制语句、列表和元组、字典和集合、字符串、Python中使用正则表达式、函数、面向对象程序设计、模块、异常处理及程序调试、文件及目录操作、操作数据库、GUI界面编程、Pygame游戏编程、网络爬虫开发、使用进程和线程、网络编程、Web编程、Flask框架、e起去旅行网站、AI图像识别工具等内容。所有知识都结合具体实例进行介绍，涉及的程序代码都给出了详细的注释，读者可轻松领会Python程序开发的精髓，快速提升开发技能。除此之外，该书还附配了243集高清教学微视频及PPT电子教案。
https://item.jd.com/13284890.html
精彩回顾
ChatGPT教你如何用Python和Matplotlib绘图（上）
微信搜索关注《Python学研大本营》
访问【IT今日热榜】，发现每日技术热点
1465
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
