Как замаскировать районы океана из моего набора данных, используя xarray, numpy и Basemap?

Question

У меня есть файл netCDF "VWC_95percentile_by_seasons.nc".Он содержит три измерения seasons (DJF, JJA, MAM, SON), latitude и longitude и имя переменной VWC 95TH PERCENTILE.

Переменная распечатывается следующим образом:

<xarray.DataArray 'VWC 95TH PERCENTILE' (season: 4, latitude: 106, longitude: 193)>
dask.array<shape=(4, 106, 193), dtype=float64, chunksize=(4, 106, 193)>
Coordinates:
  * latitude   (latitude) float32 -39.2 -39.149525 ... -33.950478 -33.9
  * longitude  (longitude) float32 140.8 140.84792 140.89584 ... 149.95209 150.0
  * season     (season) object 'DJF' 'JJA' 'MAM' 'SON'

Мне нужно будет нарисовать один график для каждого сезона на одной фигуре.Каждый участок должен быть замаскирован из районов океана.

Мой код:

import os
import time, datetime
import xarray as xr
import numpy as np
import pandas as pd
from glob import glob
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, maskoceans

shapefile1=r'{PATH}\shp\regions'

ds_path = r'{PATH}\output\VWC_95percentile_by_seasons.nc'
variable_name = 'VWC 95TH PERCENTILE'
lable = 'VWC'
title = 'VWC 95TH PERCENTILE'


print("### LOAD DATA")
ds = xr.open_mfdataset(ds_path)
da = ds[variable_name]
print(da)

print("### PLOTTING")

# define the colormap
cmap = plt.cm.jet
# extract all colors from the .jet map
cmaplist = [cmap(i) for i in range(cmap.N)]
# create the new map
cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)

# colour levels
levels=[0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]

# Plot
p = da.plot(levels=10, cmap=cmap, x='longitude', y='latitude', col='season', cbar_kwargs={'label': lable})        # p is an instance of <class 'xarray.plot.facetgrid.FacetGrid'>
plt.suptitle(title)

for i, ax in enumerate(p.axes.flatten()):

    ax.set_xlabel('')
    ax.set_ylabel('')
    map = Basemap(llcrnrlat=-39.2,urcrnrlat=-33.9,llcrnrlon=140.8,urcrnrlon=150.0,resolution='i',ax=ax)
    map.readshapefile(shapefile1, 'regions', linewidth=0.5)    # draw disctict boundaries
plt.show()

В настоящее время он нарисован:

Как я могу замаскировать все океанские районы, используя xarray, numpy и Basemap?