我正在尝试绘制归一化植被指数的线性回归图。但是使用等高线图，它也充满了海洋。我想删除没有值的海洋。
使用Nan或maskoceans
但是maskoceans并不好用。
我添加了代码。NDVI netcdf文件在这里。（https://drive.google.com/file/d/1r5N8lEQe6HP02cSz_m4edJE3AJi7lTcf/view?usp=drive_link）
我使用cdo来为netCDF文件屏蔽海洋，但是使用np.polyfit来计算线性回归使得Nan变为0（np.isnan）。这就是为什么海洋在轮廓图中着色。

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import calendar

import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
import cftime
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
%matplotlib inline

mpl.rcParams['figure.figsize'] = [8., 6.]

filename = 'E:/ERA5/ndvi331.nc'
ds = xr.open_dataset(filename)
ds

da = ds['NDVI']
da

def is_jjas(month):
    return (month >= 6) & (month <= 9)

dd = da.sel(time=is_jjas(da['time.month']))

def is_1982(year):
    return (year> 1981)

dn = dd.sel(time=is_1982(dd['time.year']))
dn

JJAS= dn.groupby('time.year').mean('time')

JJAS

JJAS2 = JJAS.mean(dim='year', keep_attrs=True)
JJAS2

fig, ax = plt.subplots(1, 1, figsize = (16, 8), subplot_kw={'projection': ccrs.PlateCarree()})
im = plt.pcolormesh(JJAS2.lon, JJAS2.lat, JJAS2, cmap='YlGn', vmin=0, vmax=1)
# Set the figure title, add lat/lon grid and coastlines
ax.set_title('AVHRR GIMMS NDVI Climatology (1982-2019)', fontsize=16)
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') 
ax.coastlines(color='black')
ax.set_extent([-20, 60, -10, 40], crs=ccrs.PlateCarree())

cbar = plt.colorbar(im,fraction=0.05, pad=0.04, extend='both', orientation='horizontal')

vals = JJAS.values

vals[np.nonzero(np.isnan(vals))] = 0
vals.shape

years = JJAS['year'].values
np.unique(years) 

years

vals2 = vals.reshape(len(years), -1)

vals2.shape

from scipy import polyfit, polyval

reg = np.polyfit(years, vals2, 1)

reg

trends = reg[0,:].reshape(vals.shape[1], vals.shape[2])

trends

trends.shape

vals.shape[1]

trends.ndim

trends.shape

np.max(trends)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm, maskoceans
from scipy.interpolate import griddata

plt.figure(figsize=(13,5))
ax = plt.subplot(111, projection=ccrs.PlateCarree()) #ccrs.Mollweide()
mm = ax.pcolormesh(dn.lon,
                   dn.lat,
                   trends,                   
                   vmin=-0.02,
                   vmax=0.02,
                   transform=ccrs.PlateCarree(),cmap='bwr' )
ax.set_global()
#ax.set_extent([-180, 180, -70, 70])
ax.coastlines();
cb=plt.colorbar(mm,ax=ax,fraction=0.046, pad=0.01)

fig, ax = plt.subplots(1, 1, figsize = (16, 8), subplot_kw={'projection': ccrs.PlateCarree()})

cs = plt.contourf(dn.lon, dn.lat, trends, levels=[-0.02, -0.015, -0.010, -0.005, 0, 0.005, 0.010, 0.015, 0.02],
                  vmin=-0.02, vmax=0.02, cmap='bwr', extend='both')

# Set the figure title, add lat/lon grid and coastlines
ax.set_title('AVHRR GIMMS NDVI Linear regression (1982-2019)', fontsize=16)
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') 
ax.coastlines(color='black')
ax.set_extent([-20, 60, -10, 40], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.OCEAN)

cbar = plt.colorbar(cs,fraction=0.05, pad=0.04, extend='both', orientation='horizontal')

字符串
enter image description here的
使用Nan或maskoceans
但是maskoceans并不好用。