Python xarray,numpy,matplotlib netcdf屏蔽海洋?

我正试图绘制NDVI的线性回归图。 但使用计数地图,它也充满了海洋。 我想去除没有价值的海洋。

使用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的工作不太好。。