Python netcdf EOF analysis

I am performing a EOF analysis for JJAS NDVI.

Here is the NDVI data (https://drive.google.com/drive/folders/1-QBkzTlAJRq8etuD15x_7l08woFTR86Y?usp=sharing) and codes.

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)

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.