我正在尝试绘制归一化植被指数的线性回归图。但是使用等高线图,它也充满了海洋。我想删除没有值的海洋。
使用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并不好用。
1条答案
按热度按时间2fjabf4q1#
你可以使用maskoceans和meshgrid来解决这个问题。
字符串
情节
的数据
注意:大陆内部的白色斑点是由maskoceans移除内陆湖泊造成的。如果你不想这样,你可以通过
inlands=False
。