mask an xarray or dataset based on specific indices

user_6610946 0
Viewed 538

I have a gridded dataset (air temperature data for example) that I would like to mask everywhere except along the coastlines (3 grid cells offshore, 3 grid cells onshore for example, a 6 grid cell buffer containing the coastline). The best way to go about this (I think) is to use a masked array containing land/ocean information.

My land/ocean mask is here. Unfortunately I cannot figure out how to recreate this file via code (I have downloaded the file and provided the file details here).

>>> 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)

I have found out how to get the indices where the mask is between 0 (ocean) and 1 (land):

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]

And now I think I have to use these indices somehow in order to create a subset of the dataset below, so that I have a subset of temperatures around the coastlines.

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

What I do not understand is how to go from here. I have looked at the xarray documentation for advanced indexing, but I still am not sure how to use that to solve my problem and any help is appreciated.

python-xarray netcdf mask python
1 Answers

I think you're on the right track with using a mask, but I don't think xarray's advanced indexing is the answer you're looking for.

I'd suggest having a look at binary_dilation and binary_erosion from scipy.ndimage: 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

If you e.g. mark the land by 1s, and the sea with 0s (or as boolean), you can see the correspondence with the scipy example:

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

So the original array looks like:

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]])

A single iteration of erosion:

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]])

So in your example, you could do:

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
user_151656000
Related