Python xarray、numpy、matplotlib netcdf 屏蔽海洋？

Question

我正在尝试绘制 NDVI 的线性回归图。 但使用Countour Map，它也充满了海洋。 我想移除没有价值观的海洋。

使用 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')

使用 Nan 或 maskoceans

但是 maskoceans 效果不佳..

Answer 1

您可以使用 maskoceans 和 meshgrid 来解决这个问题。

lon_gridded, lat_gridded = np.meshgrid(dn.lon, dn.lat)
fig, ax = plt.subplots(1, 1, figsize = (16, 8), subplot_kw={'projection': ccrs.PlateCarree()})
trends_masked = maskoceans(lon_gridded, lat_gridded, trends)
cs = plt.contourf(dn.lon, dn.lat, trends_masked, 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')

剧情

注：大陆内部的白斑是由掩蔽海洋侵蚀内陆湖泊造成的。如果您不想要这个，可以通过

inlands=False

。