写在前面
最近看论文和写论文，感觉脑袋都快成浆糊了，简直需要给自己换换脑子。一直有个想法，就是给新进组的师弟师妹们整一个基于Python的海洋气象数据处理和绘图教程。这样他们可以更快上手编程，少走弯路，避免“重复造轮子”。虽然网上也有很多教程，但是我觉得自己写应该看着更顺手，也算是回顾一些hh。顺便，也能让自己在写教程的过程中找点乐子，毕竟“教是最好的学”嘛！
虽然我会分享很多学习中收藏的网页，但作为初学者或者说跨专业的朋友（包括我自己最开始也很抗拒编程）来说，初学者总是有点“懒癌”发作。同时，尽管 python-xarray 官网就存在相关的教程且非常全面且高阶，但是我觉得对于初学者可能那并不友好，过于丰富了。
因此，我决定直接写一份手把手的实操手册，尽可能简单的语言去描述每一个步骤的目的以及获取到的信息,以及可能遇到的问题。让大家可以在短时间内“无脑”上手编程。这样一来，就能轻松跨过那些“起步难”的坑，直接奔向编程的快车道！
海洋气象编程工具-Python
在海洋与气象科学领域，编程工具的掌握对于研究生来说至关重要。在这个领域中，NCL、Matlab和Python是最为广泛使用的工具。
NCL早在气象领域应用广泛，提供了大量便捷的函数，但由于其已经停止更新，并且需要基于Linux的平台进行安装和使用，应用受限。Matlab则因其高效的矩阵计算能力而被广泛使用，但其使用受限于版权问题。
近年来，随着神经网络和深度学习模型在海洋气象研究中的广泛应用，Python逐渐成为科研高校和研究机构的首选工具。
相较之下，Python作为一种开源免费的工具，具有广泛的吸引力。通过提供各种不同的包，Python能够方便地读取海洋和气象数据。其中，xarray 是我特别推荐的。
目前，xarray 支持多种文件格式，并且能够与其他包无缝集成，使其成为海洋与气象数据分析中的一种多功能选择。接下来的部分将进一步探讨 xarray 的功能及其在海洋和气象科学数据处理与分析中的重要作用。
xarray官网：
https://docs.xarray.dev/en/stable/index.html
Xarray 介绍
本文介绍 python 软件包 xarray：
★Xarray 在类似 NumPy 的原始数组之上引入了
维度、坐标和属性形式的标签，从而为开发人员带来更直观、更简洁、更不易出错的体验。该软件包包含一个庞大且不断增长的领域无关函数库，用于使用这些数据结构进行高级分析和可视化。
★Xarray 的灵感来源于并大量借鉴了
pandas，pandas 是一款流行的数据分析软件包，主要用于标记表格数据。它特别适合处理netCDF文件（xarray 数据模型的来源），并与用于并行计算的 dask 紧密集成。
本系列尤其侧重于使用 xarray 分析存储在 netCDF 文件中的气候数据。
该系列包括以下几个部分：
NetCDF 文件读入 xarray 数据集 按时间和空间对数据集进行切片 绘图 计算一些指标，如平均值、最大值 分组和及时重采样 掩膜 将多个文件作为一个数据集打开 将数据集保存为 NetCDF
1.1 打开一个nc数据
第一步，导入 xarray 软件包
import xarray
# 或许我们可以使用简化的名称例如
# import xarray as xr
打开 netCDF 文件时，文件元数据会被读取并存储为一个 xarray.DataSet。
★这里需要注意的是，当你的文件与你的脚本在同一个目录下，你需要导入文件的相对路径，即文件名及其后缀即可；但是如果不在同一个路径下，则需要导入文件的绝对路径。
# 以绝对路径打开
ds = xarray.open_dataset(r'I:/tas_Amon_CESM2-WACCM_historical_r1i1p1f1_gn_185001-201412.nc')
# 以相对路径打开
# ds = xarray.open_dataset(r'tas_Amon_CESM2-WACCM_historical_r1i1p1f1_gn_185001-201412.nc')
# 如果你在import的部分对于xarray进行了简化的命名，则打开文件的方式为：
# ds = xr.open_dataset(r'tas_Amon_CESM2-WACCM_historical_r1i1p1f1_gn_185001-201412.nc')
该数据集的元数据现在存储在变量 ds 中。这个变量名称是自由修改命名的，但 ds 这个名称很方便！要查看内容，只需将变量名放在单元格中并运行即可，这相当于 python 程序中的 print(ds)
ds
<xarray.Dataset> Size: 438MB
Dimensions: (time: 1980, lat: 192, lon: 288, nbnd: 2)
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 355.0 356.2 357.5 358.8 time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Dimensions without coordinates: nbnd
Data variables:
tas (time, lat, lon) float32 438MB ...
time_bnds (time, nbnd) object 32kB ...
lat_bnds (lat, nbnd) float32 2kB ...
lon_bnds (lon, nbnd) float32 2kB ...
Attributes: (12/45)
Conventions: CF-1.7 CMIP-6.2
activity_id: CMIP
case_id: 4
cesm_casename: b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.001
contact: cesm_cmip6@ucar.edu
creation_date: 2019-01-30T22:37:15Z
... ...
variable_id: tas
variant_info: CMIP6 CESM2 hindcast (1850-2014) with high-top at...
variant_label: r1i1p1f1
branch_time_in_parent: 20075.0
branch_time_in_child: 674885.0
branch_method: standard
有四个部分需要注意：维度(Dimensions)、坐标(Coordinates)、数据变量(Data variables)和 属性(Attributes)。
维度 给出了每个命名维度的大小。这是一个符合 CF 标准（http://cfconventions.org）的数据集，这意味着任何与维度同名的变量都会被假定为坐标。还有其他元数据可以用来表示坐标。在这种情况下，xarray 表示与带有*的变量相关联的坐标。
数据变量列出了所有不是坐标的变量。在本例中，其中三个是边界变量(*_bnds)，定义了三个坐标的起始值和终止值。唯一真正的数据变量是 tas。
1.2 读取数据
打开数据集的命令只读取 netCDF 文件中的元数据。在进行其他操作之前，它不会尝试读取任何数据。
xarray.DataSet 对象有许多方法用于访问坐标、属性和数据。数据变量保存在一个类似于 dict(字典) 的结构中，即 ds.data_vars：
ds.data_vars
Data variables:
tas (time, lat, lon) float32 438MB ...
time_bnds (time, nbnd) object 32kB ...
lat_bnds (lat, nbnd) float32 2kB ...
lon_bnds (lon, nbnd) float32 2kB ...
只需循环数据集，就可以循环查看数据变量，数据集会依次返回每个变量的名称：
for varname in ds:
print( "读取的变量为：" + varname)
读取的变量为：tas
读取的变量为：time_bnds
读取的变量为：lat_bnds
读取的变量为：lon_bnds
可以使用单个变量的名称访问该变量，或者将其作为一个类似于 dict 的键
ds['tas']
<xarray.DataArray 'tas' (time: 1980, lat: 192, lon: 288)> Size: 438MB
[109486080 values with dtype=float32]
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8 time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
为方便使用，xarray 还提供了作为 python 属性的数据变量访问权限
ds.tas
<xarray.DataArray 'tas' (time: 1980, lat: 192, lon: 288)> Size: 438MB
[109486080 values with dtype=float32]
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8 time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
因此，ds.tas 是一个xarray.DataArray，它有自己的元数据，可提供有关变量本身的更多信息。在本例中，它是以开尔文为单位的近地面空气温度。
总结
使用xarray打开nc文件的方法为:
import xarray as xr
ds = xr.open_dataset(r'file_path')
其中，file_path是你的文件所放置的位置，请注意你的文件是否为绝对路径。
注意点： 请养成你的**文件路径中不存在中文**的好习惯，同时请将路径放置在' '中。
读取变量为：
import xarray as xr
ds = xr.open_dataset(r'file_path')
tas = ds['tas']
或者
tas = ds.tas
★https://docs.xarray.dev/en/stable/user-guide/data-structures.html
2 按时间和空间对数据集进行提取（切片和切割）
Xarray 总共支持四种不同类型的索引，如下所示，这里只需要有个概念即可：
但是我这里只介绍其中两种个人常用的：sel 和 isel
xarray 可为符合 CF 标准的数据创建带标签的坐标索引。这使得在时间和空间上选择数据子集变得简单易行。
下面介绍了数据切片的一些常见使用方法，重新加载库和数据集,这里将导入的 xarray 库重新简化命名为xr
import xarray as xr
ds = xr.open_dataset(r'I:/tas_Amon_CESM2-WACCM_historical_r1i1p1f1_gn_185001-201412.nc')
d:\anaconda\ANACONDA\envs\pji\Lib\site-packages\xarray\conventions.py:286: SerializationWarning: variable 'tas' has multiple fill values {1e+20, 1e+20} defined, decoding all values to NaN.
var = coder.decode(var, name=name)
选择数据变量 tas 并将其引用保存在另一个 python 变量中,命名为：tas
tas = ds.tas
tas
<xarray.DataArray 'tas' (time: 1980, lat: 192, lon: 288)> Size: 438MB
[109486080 values with dtype=float32]
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8 time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
xarray 构建在 numpy 的基础之上，内部将数据存储为 numpy 数组。它支持许多 numpy 操作，因此可以找出底层数据的形状，并使用 numpy 风格的索引
观测上面的变量，可以发现：
tas 是由 time x lat x lon 组成的三维数据，这里的顺序不是确定好的的，不是 latxlonxtime ，也不是 timexlonxlat ；使用 tas.shape 可以查看三个坐标对应的数量
tas.shape
(1980, 192, 288)
1980 表示时间time有1980个数值，
192 表示纬度lat有192个间隔，
288 表示经度lon有288个间隔
那么如何切片呢，比如说我想提取出第一个时刻对应的经纬度组成的网格数据
tas[0,:]
<xarray.DataArray 'tas' (lat: 192, lon: 288)> Size: 221kB
[55296 values with dtype=float32]
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
time object 8B 1850-01-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
通过只选择第一个时间索引，它创建了一个没有时间维度的 DataArray，但时间仍然是一个与任何变量都不相关联的坐标，因为它旁边不再有 *。现在它只有一个值：第一个时间索引的值。上面的索引选择等同于使用 isel，如下所示
tas.isel(time=0)
<xarray.DataArray 'tas' (lat: 192, lon: 288)> Size: 221kB
[55296 values with dtype=float32]
Coordinates:
lat (lat) float64 2kB -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0 lon (lon) float64 2kB 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
time object 8B 1850-01-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
使用方法 tas[0,:] 更简洁，使用tas.isel(time=0)方法更具描述性，但结果是一样的。
xarray 的强大之处在于数据与坐标的紧密联系。因此，可以使用等效的 .sel 操作符，但要使用坐标值。例如，要选择包括印度洋和澳大利亚在内的区域，可以使用 slice 来指明所需的经纬度值范围，并将其作为键/值对传递给 sel。slice将包括小于或等于上限值的坐标值，而不像基本 python 中的range那样排除上限值。
tas.sel(lon=slice(20,160),lat=slice(-80,25))
<xarray.DataArray 'tas' (time: 1980, lat: 112, lon: 113)> Size: 100MB
[25058880 values with dtype=float32]
Coordinates:
lat (lat) float64 896B -79.63 -78.69 -77.75 ... 23.09 24.03 24.97 lon (lon) float64 904B 20.0 21.25 22.5 23.75 ... 157.5 158.8 160.0 time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
以上命令可以串联，因此可以依次执行多个操作。例如，选择上述区域和第一次索引
tas.isel(time=0).sel(lon=slice(20,160), lat=slice(-80,25))
<xarray.DataArray 'tas' (lat: 112, lon: 113)> Size: 51kB
[12656 values with dtype=float32]
Coordinates:
lat (lat) float64 896B -79.63 -78.69 -77.75 ... 23.09 24.03 24.97 lon (lon) float64 904B 20.0 21.25 22.5 23.75 ... 157.5 158.8 160.0
time object 8B 1850-01-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
在这种情况下，使用 isel 选择时间比指定日期更方便，但也可以使用 sel 明确指定日期
tas.sel(time='1850-01-15T12:00:00', lon=slice(20,160), lat=slice(-80,25))
<xarray.DataArray 'tas' (time: 1, lat: 112, lon: 113)> Size: 51kB
[12656 values with dtype=float32]
Coordinates:
lat (lat) float64 896B -79.63 -78.69 -77.75 ... 23.09 24.03 24.97 lon (lon) float64 904B 20.0 21.25 22.5 23.75 ... 157.5 158.8 160.0 time (time) object 8B 1850-01-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
也可以在 时间 维度上使用 slice ,选择某个时间范围内的数据，比如说选择 1871 年 3 月至 11 月：
tas.sel(time=slice('1871-03','1871-11'), lon=slice(20,160), lat=slice(-80,25))
<xarray.DataArray 'tas' (time: 9, lat: 112, lon: 113)> Size: 456kB
[113904 values with dtype=float32]
Coordinates:
lat (lat) float64 896B -79.63 -78.69 -77.75 ... 23.09 24.03 24.97 lon (lon) float64 904B 20.0 21.25 22.5 23.75 ... 157.5 158.8 160.0 time (time) object 72B 1871-03-15 12:00:00 ... 1871-11-15 00:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
slice - 切片 运算符选择上界和下界之间的值。如果使用 sel 时需要单个坐标值，它必须与坐标数组中的精确值相对应，或者使用指定的 method 参数告诉 xarray 如何选择值。例如，只选择最靠近布里斯班的单元格中的值
tas.sel(lat=-27.47, lon=153.03, method='nearest')
<xarray.DataArray 'tas' (time: 1980)> Size: 8kB
[1980 values with dtype=float32]
Coordinates:
lat float64 8B -27.8
lon float64 8B 152.5
time (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
Attributes: (12/19)
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
description: near-surface (usually, 2 meter) air temperature
frequency: mon
id: tas
... ...
time_label: time-mean
time_title: Temporal mean
title: Near-Surface Air Temperature
type: real
units: K
variable_id: tas
因此，数据中最接近的位置位于 lat=-27.8 , lon=152.5。
总结
★以上内容在 .ipynb 进行整理，偷懒这里直接将其转换为markdown，可能会影响部分内容的展示。打算全部写完再分享相关文件和数据，希望在暑假结束前整理完~
> https://docs.xarray.dev/en/stable/user-guide/indexing.html
> https://docs.xarray.dev/en/stable/index.html
给我点在看的人
越来越好看
