Python netcdf EOF analysis

Question

I am performing a EOF analysis for JJAS NDVI.

enter image description here

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
import scipy.signal
from scipy import signal
import numpy.polynomial.polynomial as poly
from eofs.standard import Eof
mpl.rcParams['figure.figsize'] = [8., 6.]
filename = 'D:/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
def lon_jjas(lon):
 return (lon >= -15) & (lon <= 20)
JJ_lon = dn.sel(lon=lon_jjas(dn['lon']))
JJ2 = JJ_lon.mean(dim='lon', keep_attrs=True)
JJ2
def lat_jjas(lat):
 return (lat >= 11) & (lat <= 20)
JJ_lat = JJ2.sel(lat=lat_jjas(JJ2['lat']))
JJ3 = JJ_lat.mean(dim='lat', keep_attrs=True)
JJ3 = JJ3.groupby('time.year').mean('time')
JJ3
import scipy.signal
JJ4=scipy.signal.detrend(JJ3)
JJ4
JJ4m= np.mean(JJ4, axis=0)
an4 = JJ4 - JJ4m
lons = ds['lon'].values
lats = ds['lat'].values
# Create an EOF solver to do the EOF analysis. Square-root of cosine of
# latitude weights are applied before the computation of EOFs.
coslat = np.cos(np.deg2rad(lats)).clip(0., 1.)
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(an4, weights=wgts)
eof1 = solver.eofs(neofs=10)
pc1 = solver.pcs(npcs=10, pcscaling=0)
varfrac = solver.varianceFraction()
lambdas = solver.eigenvalues()

However 'solver = Eof(an4, weights=wgts)' has an error. So I tried to reshape JJA4 data but it is not working well.

How can I get eof results by netCDF?

---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[27], <a href='vscode-notebook-cell:?execution_count=27&line=1'>line 1</a>
----> <a href='vscode-notebook-cell:?execution_count=27&line=1'>1</a> solver = Eof(an4, weights=wgts)
File ~\anaconda3\Lib\site-packages\eofs\standard.py:112, in Eof.__init__(self, dataset, weights, center, ddof)
 <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:110'>110</a> # Store the input data in an instance variable.
 <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:111'>111</a> if dataset.ndim < 2:
--> <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:112'>112</a> raise ValueError('the input data set must be '
 <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:113'>113</a> 'at least two dimensional')
 <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:114'>114</a> self._data = dataset.copy()
 <a href='~\anaconda3\Lib\site-packages\eofs\standard.py:115'>115</a> # Check if the input is a masked array. If so fill it with NaN.
ValueError: the input data set must be at least two dimensional

Actually I don't understand what function 'Eof()' should be consisted of. I found out that JJ4(detrended anomaly) should have time and space dimension. Adding and Transferring dimensional data, still error is interrupting me!

I would like to plot PC1 & PC2 (or PC1, PC2, PC3)...

https://github.com/royalosyin/Python-Practical-Application-on-Climate-Variability-Studies/blob/master/ex18-EOF%20analysis%20global%20SST.ipynb https://ajdawson.github.io/eofs/latest/userguide/index.html

# an4 as DataArray
an4_da = xr.DataArray(an4, dims=["time", "lat", "lon"], coords={"time": ds['time'], "lat": ds['lat'], "lon": ds['lon']})
# Transfer into 2D (time x space)
an4_stacked = an4_da.stack(space=('lat', 'lon'))
# NaN
an4_stacked = an4_stacked.fillna(0)

Answer 1

You flagged your question cdo-climate so here is a quick way to do this with cdo is this:

export CDO_WEIGHT_MODE=off
export MAX_JACOBI_ITER=100
cdo eof,n input.nc eigval.nc pattern.nc
cdo eofcoeff pattern.nc input.nc pc.nc
cdo -div eigval.nc -timsum eigval.nc /explvar.nc

where n is the number of EOFs you want as output. More details here, and you can also call cdo directly from within python using the cdo package.

• Mincheol: Thank you so much!! I would try the solution. If I want only JJAS NDVI, then I should slice them first. right?
• ClimateUnboxed: yes you need to extract the months you want, cdo selmon,6/9 in out I think
• Mincheol: Great! I am sorry for late accept