Python netcdf EOF分析

Question

我正在对 JJAS NDVI 进行 EOF 分析。

这里是 NDVI 数据（https://drive.google.com/drive/folders/1-QBkzTlAJRq8etuD15x_7l08woFTR86Y?usp=sharing）和代码。

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()

但是“solver = Eof(an4,weights=wgts)”有错误。 所以我尝试重塑 JJA4 数据，但效果不佳。

如何通过 netCDF 获得 eof 结果？

---------------------------------------------------------------------------
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

其实我不明白'Eof()'应该由什么函数组成。 我发现JJ4（去趋势异常）应该有时间和空间维度。 添加和传输尺寸数据，仍然出现错误打扰我！

我想绘制 PC1 和 PC2（或 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

您标记了您的问题 cdo-climate，所以这里有一个使用 cdo 执行此操作的快速方法：

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

其中 n 是您想要作为输出的 EOF 数量。更多详细信息这里，您还可以使用 cdo 包直接从 python 中调用 cdo。