好奇心Log 发表于 2024-3-15 22:11:33

数据处理 | 经验正交函数(EOF)与旋转经验正交函数(REOF)

本帖最后由 好奇心Log 于 2024-3-15 22:13 编辑

导入模块
from pyEOF import *
import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as pltWarning: ecCodes 2.22.0 or higher is recommended. You are running version 2.21.0

定义绘图函数
# create a function for visualization convenience
def visualization(da, pcs, eofs_da, evf):
    fig = plt.figure(figsize=(12, 12))

    ax = fig.add_subplot(n+1, 2, 1)
    da.mean(dim=["lat", "lon"]).plot(ax=ax)
    ax.set_title("average air temp")

    ax = fig.add_subplot(n+1, 2, 2)
    da.mean(dim="time").plot(ax=ax)
    ax.set_title("average air temp")

    for i in range(1, n+1):
      pc_i = pcs["PC"+str(i)].to_xarray()
      eof_i = eofs_da.sel(EOF=i)["air"]
      frac = str(np.array(evf*100).round(2))

      ax = fig.add_subplot(n+1, 2, i*2+1)
      pc_i.plot(ax=ax)
      ax.set_title("PC"+str(i)+" ("+frac+"%)")

      ax = fig.add_subplot(n+1, 2, i*2+2)
      eof_i.plot(ax=ax,
                   vmin=-0.75, vmax=0.75, cmap="RdBu_r",
                   cbar_kwargs={'label': ""})
      ax.set_title("EOF"+str(i)+" ("+frac+"%)")

    plt.tight_layout()
    plt.show()
读取气温数据
# load the DataArray
da = xr.tutorial.open_dataset('air_temperature')["air"]
print(da)

# create a mask
mask = da.sel(time=da.time)
mask = mask.where(mask<250).isnull().drop("time")

# get the DataArray with mask
da = da.where(mask)
da.sel(time=da.time).plot()
plt.show()<xarray.DataArray 'air' (time: 2920, lat: 25, lon: 53)>

Coordinates:
* lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
* lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
* time   (time) datetime64 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:   4xDaily Air temperature at sigma level 995
    units:         degK
    precision:   2
    GRIB_id:       11
    GRIB_name:   TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:   Individual Obs
    parent_stat:   Other
    actual_range:



# convert DataArray to DataFrame
df = da.to_dataframe().reset_index() # get df from da
display(df.head(5))
print("DataFrame Shape:",df.shape)


DataFrame Shape: (3869000, 4)

# reshape the dataframe to be
df_data = get_time_space(df, time_dim = "time", lumped_space_dims = ["lat","lon"])
display(df_data.head(5))
print("DataFrame Shape:",df_data.shape)

5 rows × 1325 columns

DataFrame Shape: (2920, 1325)

REOF分解
n = 4
pca = df_eof(df_data,pca_type="varimax",n_components=n)

eofs = pca.eofs(s=2, n=n) # get eofs
eofs_da = eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization
pcs = pca.pcs(s=2, n=n) # get pcs
evfs = pca.evf(n=n) # get variance fraction

# plot
visualization(da, pcs, eofs_da, evfs)/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:847: RuntimeWarning: invalid value encountered in true_divide
updated_mean = (last_sum + new_sum) / updated_sample_count
/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:687: RuntimeWarning: Degrees of freedom <= 0 for slice.
result = op(x, *args, **kwargs, dtype=np.float64)




传统EOF分解
n = 4 # define the number of components

pca = df_eof(df_data) # implement EOF

eofs = pca.eofs(s=2, n=n) # get eofs
eofs_da = eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization
pcs = pca.pcs(s=2, n=n) # get pcs
evfs = pca.evf(n=n) # get variance fraction

# plot
visualization(da, pcs, eofs_da, evfs)/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:847: RuntimeWarning: invalid value encountered in true_divide
updated_mean = (last_sum + new_sum) / updated_sample_count
/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:687: RuntimeWarning: Degrees of freedom <= 0 for slice.
result = op(x, *args, **kwargs, dtype=np.float64)




eofs模块分解结果
from eofs.standard import Eof
from sklearn.preprocessing import StandardScaler
solver = Eof(StandardScaler().fit_transform(df_data.values))

s_pcs = pd.DataFrame(data=solver.pcs(npcs=4, pcscaling=2),
                     columns = pcs.columns,
                     index = pcs.index)


s_eofs = pd.DataFrame(data = solver.eofs(neofs=4, eofscaling=2),
                      columns = eofs.columns,
                      index = eofs.index)
s_eofs_da = s_eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization

s_evfs = solver.varianceFraction(neigs=4)

# plot
visualization(da, s_pcs, s_eofs_da, s_evfs)/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:847: RuntimeWarning: invalid value encountered in true_divide
updated_mean = (last_sum + new_sum) / updated_sample_count
/home/lqy/anaconda3/lib/python3.8/site-packages/sklearn/utils/extmath.py:687: RuntimeWarning: Degrees of freedom <= 0 for slice.
result = op(x, *args, **kwargs, dtype=np.float64)




本篇文章来自于微信公众号:好奇心Log


FrankJScott 发表于 昨天 23:49

Best Cigna Dietitian Guide

In reply to the lady inquiring about health coach job description, health and lifestyle coach, find health coach, wellness nutrition coach, exercise and nutrition plan, holistic health and nutrition coach, dietitian to lose weight, the health coach institute, fitness health coach, certified health coach certification,I highly suggest this find on weight loss nutritionist blog or becoming a holistic health coach, dietitian menu, nutritionist free, wellness and health coach certification, weight coach, in health coach, holistic health coach programs, nutrition coach job description, working as a health coach, health and life coach, which is worth considering with this here for diabetic dietitian near me tips alongside all health wellness coach, meal plans and workout programs, difference between nutrition coach and nutritionist, nutrition coaching institute reviews, find a nutrition coach, best diet plan app for weight loss free, health coach testimonials, online health coaching services, health life coach, fitness and nutrition plan, as well as this continue on blue cross blue shield nutritionist link which is also great. Also, have a look at this going here for united healthcare dietitian details as well as wellness coaching, workout food plan, best diet plan app for weight loss, wellness coaching packages, certified health coach certification, in holistic health coach, exercise and eating plan, nutrition and lifestyle coach, nutrition health coach certification, health care coaching, and don't forget this high rated blue cross blue shield dietitian tips bearing in mind best holistic health coach certification, nutrition coaching, personalized nutritionist, holistic food coach, holistic wellness coach near me,weblink for alongside all health coaching packages prices, best health coaching programs, diet chart for a normal healthy person, working with a nutritionist to lose weight, working as a health coach,for good measure. Check more @ Updated Email List Info 670e706
页: [1]
查看完整版本: 数据处理 | 经验正交函数(EOF)与旋转经验正交函数(REOF)