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

编程问答 2022-05-25

问题描述

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

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

>>> mask.shape
(1,81,1440)
>>> mask
masked_array(
  data=[[[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 文档，但我仍然不确定如何使用它来解决我的问题，我们不胜感激。

解决方法

我认为您使用掩码是正确的，但我不认为 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,1,0]])

侵蚀的单次迭代：

ndimage.binary_erosion(a).astype(a.dtype)
array([[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
masknetcdfpythonpython-xarray