根据特定索引屏蔽 xarray 或数据集答案

【问题标题】：mask an xarray or dataset based on specific indices根据特定索引屏蔽 xarray 或数据集
【发布时间】：2021-05-21 08:38:58
【问题描述】：

我有一个网格数据集（例如气温数据），我想在除海岸线以外的任何地方进行屏蔽（例如，离岸 3 个网格单元，陆上 3 个网格单元，包含海岸线的 6 个网格单元缓冲区）。解决这个问题的最好方法（我认为）是使用包含陆地/海洋信息的掩码数组。

我的陆地/海洋面具在这里。不幸的是，我无法弄清楚如何通过代码重新创建此文件（我已经下载了文件并在此处提供了文件详细信息）。

>>> mask.shape
(1, 81, 1440)
>>> mask
masked_array(
  data=[[[0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ...,
          0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
         [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ...,
          0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
         [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ...,
          0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
         ...,
         [9.4939685e-01, 7.2325826e-01, 4.7175312e-01, ...,
          9.9737060e-01, 9.8022592e-01, 9.5705426e-01],
         [7.6317883e-01, 5.6850266e-01, 2.9611492e-01, ...,
          9.9931705e-01, 9.8200846e-01, 9.4740820e-01],
         [3.4933090e-01, 1.7199135e-01, 3.2305717e-05, ...,
          9.5507932e-01, 7.9238594e-01, 5.7203436e-01]]],
  mask=False,
  fill_value=1e+20,
  dtype=float32)

我发现了如何获取掩码在 0（海洋）和 1（陆地）之间的索引：

indices_btwn_0_1 = np.where(np.logical_and(mask>0,mask<1))
lat_indices_btwn_0_1 = indices_btwn_0_1[1] #
lon_indices_btwn_0_1 = indices_btwn_0_1[2]

现在我认为我必须以某种方式使用这些指数来创建下面数据集的子集，这样我就有了海岸线周围的温度子集。

ds = xr.open_dataset(ifile_climate_var_data)
>>> ds
<xarray.Dataset>
Dimensions:    (latitude: 81, longitude: 1440, time: 8760)
Coordinates:
  * longitude  (longitude) float32 -180.0 -179.75 -179.5 -179.25 -179.0 ...
  * latitude   (latitude) float32 85.0 84.75 84.5 84.25 84.0 83.75 83.5 ...
  * time       (time) datetime64[ns] 2019-01-01 2019-01-01T01:00:00 ...
Data variables:
    t2m        (time, latitude, longitude) float32 253.0178 252.98686 ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-02-15 19:05:54 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...

temp_2m = ds.t2m

我不明白如何从这里走。我已经查看了有关高级索引的 xarray 文档，但我仍然不确定如何使用它来解决我的问题，感谢任何帮助。

【问题讨论】：

标签： python mask netcdf python-xarray

【解决方案1】：

我认为您使用掩码的方法是正确的，但我不认为 xarray 的高级索引是您正在寻找的答案。

我建议看看 scipy.ndimage 中的 binary_dilation 和 binary_erosion： https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.binary_dilation.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.binary_erosion.html

如果你例如用 1 标记陆地，用 0 标记海洋（或作为布尔值），您可以看到与 scipy 示例的对应关系：

from scipy import ndimage
a = np.zeros((7,7), dtype=int)
a[1:6, 2:5] = 1
a

所以原来的数组是这样的：

array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])

侵蚀的单次迭代：

ndimage.binary_erosion(a).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])

所以在你的例子中，你可以这样做：

is_land = (mask > 0)  # this might not be the right condition
eroded = is_land.copy(data=ndimage.binary_erosion(is_land.values, iterations=3)
dilated = is_land.copy(data=ndimage.binary_dilation(is_land.values, iterations=3)
coastal_zone = eroded | dilated

【讨论】：