Python xarray, numpy, matplotlib netcdf masking ocean?

Question

I am trying to draw linear regression map of NDVI. But using countour map, it filled the ocean too. I would like to remove the ocean without values.

using Nan or maskoceans

but maskoceans is not working well..

I added the code. NDVI netcdf file is here. ( https://drive.google.com/file/d/1r5N8lEQe6HP02cSz_m4edJE3AJi7lTcf/view?usp=drive_link )

I used cdo to mask ocean for netCDF file but using np.polyfit to calculate linear regression made the Nan to 0 (np.isnan). That's why the Ocean colored in the countour map.

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

using Nan or maskoceans

but maskoceans is not working well..

Answer 1

You can solve this using maskoceans and 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')

Plot

Note: the white spots inside the continent are caused by maskoceans removing inland lakes. You can pass inlands=False if you do not want this.