matlab经验正交分解函数EOF的实现—基于Climate Data Toolbox操作

已于 2023-04-14 19:28:09 修改
阅读量1.3w 收藏 162
点赞数 52
分类专栏： matlab 文章标签： matlab 经验正交分解
于 2022-03-11 20:35:43 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_43637490/article/details/123418759
版权

经验正交函数分析方法(Empirical Orthogonal Function,缩写为EOF)，也称特征向量分析(Eigenvector Analysis)，或者主成分分析(Principal Component Analysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。EOF分析方法能够把随时间变化的变量场分解为不随时间变化的空间函数部分以及只依赖时间变化的时间函数部分。空间函数部分概括场的地域分布特点，而时间函数部分则是由场的空间点的变量线性组合所构成，称为主要分量。
百度百科

以下内容未经允许请不要随意转载

参考视频：【李华的研究室】经验正交方程EOF-原理，代码，气象学实例，文献解读
这位UP主关于EOF的原理及用法解释的比较详细，关于这个视频做了一些笔记如下：
在这里插入图片描述
在这里插入图片描述
视频里UP主说关于纬度权重的问题，他考虑到格网的距离是随纬度变化的，所以给原始数据乘上了纬度的权重。这个看个人吧。关于EOF单位量纲的处理很到位。

在网上找了很多代码，基本都是直接根据EOF思路编写的，总是感觉有些欠缺，所以找到Climate Data Toolbox（CDT）工具箱里的EOF函数来做。
工具箱下载地址：Climate Data Toolbox for Matlab
EOF使用文档：eof documentation
eof documentation
The eof function gives eigenmode maps of variability and corresponding principal component time series for spatiotemporal data analysis. This function is designed specifically for 3D matricies of data such as sea surface temperatures where dimensions 1 and 2 are spatial dimensions (e.g., lat and lon; lon and lat; x and y, etc.), and the third dimension represents different slices or snapshots of data in time.

Syntax

eof_maps = eof(A)
eof_maps = eof(A,n)
eof_maps = eof(…,‘mask’,mask)
[eof_maps,pc,expvar] = eof(…)

Description

最简单的：
eof_maps = eof(A) calculates all modes of variability in A, where A is a 3D matrix whose first two dimensions are spatial, the third dimension is temporal(要求输入的矩阵为三维矩阵，前两维是空间，第三维时间), and data are assumed to be equally spaced in time. Output eof_maps have the same dimensions as A, where each map along the third dimension represents a mode of variability order of importance.

自定义导出多少模态的：
eof_maps = eof(A,n) only calculates the first n modes of variability. For large datasets, it’s computationally faster to only calculate the number of modes you’ll need. If n is not specified, all EOFs are calculated (one for each time slice).

自定义掩膜区域的：
eof_maps = eof(…,‘mask’,mask) only performs EOF analysis on the grid cells represented by ones in a logical mask whose dimensions correspond to dimensions 1 and 2 of A. This option is provided to prevent solving for things that don’t need to be solved, or to let you do analysis on one region separately from another. By default, any grid cells in A which contain any NaNs are masked out.

常用下面这个
[eof_maps,pc,expvar] = eof(…) returns the principal component time series pc whose rows each represent a different mode from 1 to n and columns correspond to time steps. For example, pc(1,:) is the time series of the first (dominant) mode of varibility. The third output expvar is the percent of variance explained by each mode. See the note below on interpreting expvar.

CDT工具箱EOF简单介绍

在直接使用之前希望能把下面这一部分阅读一下，这决定了能否做出EOF以及做出结果与他人结果之间的差别。

这个工具箱直接集成了EOF计算过程，其实代码很少比较简单，但是他考虑了一些细节部分，比起直接根据思路编写的代码比起来很有用，下面列举一些细节处理部分。这些部分也决定了其做出的结果与其他算法结果不一样。

1、首先是无数据区域的掩膜部分。我这里处理的是海洋部分，所以下载的nc数据陆地部分表示为Nan。如果直接对带着Nan的数据进行协方差矩阵计算，得到的协方差矩阵也会全变成Nan。所以这个EOF函数直接在开头就考虑到这一点，提取了数据中Nan的区域，让它作为“Mask”。然后在后续处理中将无数据区域掩膜掉。
在这里插入图片描述在这里插入图片描述
但是这个mask也存在一些问题。可以看到生成mask的代码是~any(isnan(A),3)，这就意味着，一旦一个格点在某一时间点是无数据的，那么这个格点整个时间的数据在后期都会被忽略，也就没法做出EOF结果，可以参考我下面这个结果图。
在这里插入图片描述
可以看到大部分区域都是没结果的状态，然后我去看了他生成的mask，发现那些区域的格点在某个时间点没数据，结果这个格点都被忽略了。

2、这个函数里关于数据去趋势，也即求距平，他用的是detrend函数。
在这里插入图片描述
这个时候的A已经被转换为二维数组，时间×空间，也就是每一行代表某一时间点的所有空间数据点，每一列代表某一空间点所有时间的数据。查询detrend的使用方法可知，它是按照列的顺序去趋势，那也就是对该点所有时间数据去趋势。我看到有些代码写的是对某一时间所有空间内数据去趋势，所以这里我也不知道到底该怎么做。这是需要注意的一点。

3、关于求协方差矩阵，这个函数很细心的考虑到了时间和空间的维度。但是他没有除以n时间数，在视频里UP主说除以或者不除以只会影响最后的特征值（而且会影响EOF的数据范围，我试了一下，如果除以n，EOF会很大而特征值很小）。
在这里插入图片描述在这里插入图片描述

以上差不多是这个工具箱函数所涉及到的一些细节问题，了解之后便可以开始使用了。

CDT工具箱EOF函数使用

根据EOF工具可以直接进行分析处理，但是还涉及到数据预处理的一部分，所以下面的代码很大一部分是关于数据处理和画图的。提一句画图用的是m_map，但是感觉效果不好，最后结果分析部分会写。
以下是一些自己写代码过程中遇到的问题记录：

1、注意对原始数据的处理。读取nc文件的时候先看一下他各个变量的名称，有的会简写有的会全拼，有的会大写有的会小写。他这个EOF函数需要三维矩阵，空间和时间。我下载的nc数据是四维的，包含深度维度，所以我先需要把深度处理一下，用的reshape函数。我的数据还涉及到一个问题，就是直接读取的nc文件和真实的地理分布情况不一样。可以看一下下图（我直接复制数据到excel里看一下效果）。
在这里插入图片描述这很明显处于一个颠倒的状态，所以要根据自己的实际区域去做调整，比如下面代码里我直接逆时针旋转90度rot90即可。有的数据诡异到需要左右上下翻转（有个数据需要flipud(rot90())）之后才对得上。下图是翻转之后的结果。
在这里插入图片描述2、关于North检验
在这里插入图片描述
因为后续我要做North检验，所以我修改了工具箱EOF函数的一部分，我把EOF计算的特征值作为函数输出了，也就是变成了[eofs_CHL,pcs_CHL,expvar_CHL,lamda] = eof(CHL_True);，这个需要在EOF函数里进行修改，在EOF的最后添加lamda = diag(D);，然后在对应的各个函数输出里添加上即可。
在这里插入图片描述在这里插入图片描述
在这里插入图片描述

3、m_grid画边界的时候，起始和终止经纬度一定得是整数，如果是小数会报错。
4、代码中涉及到EOF结果导出tif，后续会提到，如果不需要可以直接去掉。

clear all
clc
data = 'F:\0TCLI\Marine Data Analysis\Data\CHL\cmems_mod_glo_bgc_my_0.25_P1M-m_20112020CHLetcl.nc';  %10年月均数据
ncdisp(data);
ncinfo(data);
lon = ncread(data,'longitude');  % 经度longitude
lat = ncread(data,'latitude');  % 纬度latitude
time = ncread(data,'time');  % 时间
CHL_data = ncread(data,'chl');  % 叶绿素a Chlorophyll-a
ntimes = length(time); nlats = length(lat); nlons = length(lon);

%% CHL数据处理
CHL_3D = reshape(CHL_data,size(CHL_data,1),size(CHL_data,2),[]);  % 四维数组变三维，深度这一维度不起作用
CHL_True = rot90(CHL_3D);  % 逆时针旋转90度变成真正符合实际情况的形式
boundary = [106 120 3 24];  % [起始经度 终止经度 起始纬度 终止纬度]


%% 做EOF，及结果导出
% 这个考虑了维度变化权重
% lat_1 = linspace(boundary(4),boundary(3),nlats)';
% [eofs_CHL,pcs_CHL,expvar_CHL,lamda] = eof(CHL_True.*sqrt(cosd(lat_1))); 

[eofs_CHL,pcs_CHL,expvar_CHL,lamda] = eof(CHL_True);  % 使用CDT工具箱里的EOF函数
EOF_1_CHL = eofs_CHL(:,:,1).*std(pcs_CHL(1,:));  % 第一模态标准化
PC_1_CHL = pcs_CHL(1,:)./std(pcs_CHL(1,:));  % 第一模态时间系数
[A,GeoReference] = geotiffread('F:\0TCLI\Marine Data Analysis\Data\CHL\CHL_TEST.tif');  % 我要把EOF结果导出
geotiffwrite('CHL_EOF_1_Origin.tif',eofs_CHL,GeoReference);  % EOF结果导出tif

EOF_2_CHL = eofs_CHL(:,:,2).*std(pcs_CHL(2,:));  % 第二模态
PC_2_CHL = pcs_CHL(2,:)./std(pcs_CHL(2,:));

% North检验
mask = ~any(isnan(CHL_True),3);  % 获取Nan区域
N = sum(sum(double(mask) == 1));  % 统计有效数据点数
d = lamda.*sqrt(2/N);  % 误差
delta_lamda(:,1) = lamda-d;  % 特征值下限
delta_lamda(:,2) = lamda+d;  % 特征值上限

%% 绘图
m_proj('Mercator','lat',[boundary(3) boundary(4)],'lon',[boundary(1) boundary(2)]);  % 定义投影
% 生成网格
lat_2 = linspace(boundary(4),boundary(3),nlats);  % 要从高纬度到低纬度生成（我是北半球）
lon_2 = linspace(boundary(1),boundary(2),nlons); 
[plon,plat] = meshgrid(lon_2,lat_2);
% 绘制图形
subplot(1,2,1)
m_pcolor(plon,plat,EOF_1_CHL)             % 添加要画的内容
colormap('jet');  % m_colmap('jet')  % 设置色带的颜色
% shading interp  % 这个可以让结果图很平滑
% caxis([-0.1 0.1])  % 这个设置显示范围
hold on
% 添加边界，我这里导入了外部的矢量文件边界
% 可以直接用m_coast('color',[0 0 0],'linewidth',2) 绘制海岸线，填充陆地
Border = shaperead('E:\0Data\世界国家\世界国家.shp');
bor_x = [Border(:).X];
bor_y = [Border(:).Y];
m_plot(bor_x,bor_y,'linewidth',2,'color',[0 0 0]);
m_grid('box','fancy')                    %添加边框
hold on
%添加标题
title('EOF First','FontSize',12,'FontName','Times New Roman') 
%添加色带
h = colorbar('eastoutside','FontSize',12,'FontName','Times New Roman');
% 下面的可以为色带添加小标题，写上单位，还能修改色带位置让它离图远一点
% set(get(h,'title'),'string','mmol m^{-3}','FontSize',12,'FontName','Times New Roman');
% position = get(h,'Position');
% set(h,'Position', [position(1)+0.15,position(2),position(3),position(4)]);

% 请注意这里的hours和1950这些要根据nc文件里的time来确定
% 如果设置不对则坐标轴显示的时间不对
time_true = hours(time) + datetime('1950-01-01 00:00:00');
formatOut = 'yyyy-mmm';
time_label = string(datestr(time_true,formatOut));  % 设置显示的时间标签
time_str = datevec(time_true);
time_year = unique(time_str(:,1));

subplot(1,2,2)
anomaly(1:numel(time_label),PC_2_CHL)  % 前面已经变成了负的
xticks([1:12:numel(time_label)])
xticklabels(time_label(1:12:end,:))  % 设置显示的x标签
xtickangle(45)  % 标签斜45度
xticks([1:12:numel(time_label)])
% plot(date,PC_2_CHL)  % 也能直接画

20230414更新
原始代码中画的x轴标签不正确，修改了x轴的相关代码，效果如下，可以根据自己需要自行修改
在这里插入图片描述

EOF结果分析

下图即为第一模态和时间系数结果，我用的是-EOF和-PC，为啥加上负号可以参考开头提到的视频。
在这里插入图片描述
各个模态方差贡献率和特征值误差：
在这里插入图片描述

可以看到，整个海域的Chl-a变化好像并不太明显，除了靠近陆地的那个区域变化剧烈，还有中西部有一点变化。时间系数也有很明显的季节变化规律。前五个模态都通过了North检验。
但是这个结果图让我看着太难受了，给人一种整个海域都没多大变化的感觉。所以要对结果进行导出处理。代码里也涉及到了导出，导出为tif文件需要地理参考，这个参考从哪来？就从原始nc文件里来，在arcgis中先打开nc文件，然后导出为tif文件，再导入matlab中获取它的地理参考，即可当作EOF文件的地理信息。geotiffread、geotiffwrite

导出之后在arcgis中打开，做直方图均值化处理，效果直接提升了一个档次。
在这里插入图片描述

范围都是一样的，但是给人感觉不一样。可以看到南海中西部很明显有剧烈变化，这跟已有文献的结果差不多（虽然我的范围有点太大了，但是感觉还行吧）。
在这里插入图片描述
在这里插入图片描述
那这里就涉及到一个结果显示的问题，matlab也能做直方图均衡化啊，直接在matlab里histeq不完了。
首先我们的EOF结果中是存在Nan值的，直接histeq会显示：
在这里插入图片描述
他会自动把Nan变成0值，然后去做。0也作为了像元的值参与了直方图均值化，直接影响了结果。而且陆地被均值化之后还要用之前的mask去处理掉不让他显示
在这里插入图片描述
显示结果如下图所示。左边是matlab直方图均值化，右边是arcgis直方图均值化。可以看到中部区域在左图显得变化很剧烈，但其实只有中西部一小部分区域是剧烈变化的。
在这里插入图片描述

再回到EOF结果导出这，导出为tif之后直接在arcgis中打开会成这样：
在这里插入图片描述
也就是陆地那些Nan值会变成0然后被赋予颜色。那我不想让他显示，需要在ENVI里处理一下。ENVI APP STORE中找个插件，或者参考这个ENVI扩展工具：利用波段运算修改NaN方法总结
在这里插入图片描述
然后用Band Math工具把Nan赋成0，再修改metadata，设置Data Ignore Value为0，然后导出dat文件，再在arcgis中打开即可忽略背景。

EOF还原

CDT提供了还原为原始数据的函数：reof
在这里插入图片描述只是写的时候第三个参数一定写1：模态总数，如果只写一个那就只还原那个模态。
如果EOF里去趋势了，那么还原出来的也是去趋势的数据。

% Detrend_CHL是从EOF函数里输出出来的
mask = ~any(isnan(CHL_True),3);  % 去除Nan值
Detrend_CHL_maps = nan(size(Detrend_CHL,1),numel(mask));   
Detrend_CHL_maps(:,mask(:)) = Detrend_CHL;
Detrend_CHL_maps1 = reshape(Detrend_CHL_maps,size(Detrend_CHL,1),size(mask,1),size(mask,2));
Origin_CHL_De = permute(Detrend_CHL_maps1,[2 3 1]);

到此，EOF处理算是差不多了。

PeanutbutterBoh
关注
专栏目录
Matlab实现EOF分析
dm1901的博客
04-02 2893
这里使用了Matlab的eig函数，计算出协方差矩阵的特征值和特征向量。这里使用了Matlab的mean函数，计算了数据的平均值，同时使用repmat函数将平均值扩展为一个与数据矩阵相同大小的矩阵，然后使用减法得到去平均化的数据矩阵data_anom。一般而言，需要计算出数据的平均值，并将原始数据减去平均值，得到去平均化的数据。这里使用了Matlab的矩阵乘法运算符，将去平均化的数据矩阵和选定的特征向量矩阵相乘，得到EOF模态矩阵eof_mode。这里选择前3个特征向量。首先，需要导入待分析的数据。
MATLAB eof用法,[转载]基于Matlab软件进行EOF分解、回归趋势分析
weixin_29475313的博客
03-16 3044
clearcloseclcncload('E:/data/mon/ndvi79-06.nc','time','lat','lon','ndvi');%第一部分：数据处理，剔除缺失值，求距平，并修正为等权重%下面均针对温娜的sst进行标注%选取区域：纬度：-9.5-59.5 ；经度：60.5-149.5ilat=find(lat>=-10 & lat<60);nlat=lengt...
经验正交函数 (EOF) 分析：经验正交函数 (EOF) 分析常用于气象学和气候学-matlab开发
05-30
在统计和信号处理中，经验正交函数 (EOF) 分析方法是根据从数据确定的正交基函数对信号或数据集进行分解。 它与对数据执行主成分分析相同，只是 EOF 方法可以找到时间序列和空间模式。 该术语也可与地球物理学中的地理加权 PCA 互换。
EOF.zip_EOF正交分解_EOF经验正交_eof_物理海洋_经验正交函数
07-15
经验正交分解EOF,经验正交函数现已成为大气，海洋及一切地球物理科学中数据分析的有力工具。
基于matlab的各种EOF分解程序
04-02
基于MATLAB分析代码源程序，分析程序 ，
EOF的MATLAB程序
05-04
function [E,V,A,C]=eeof(X, M, convert) % Syntax: [E,V,A,C]=eeof(X, M); [E,V,A,C]=eeof(X, M, 1); % This function performs an extended empirical orthogonal % function (EEOF) analysis of matrix 'X', for embedding dimension 'M'. % Each of the L columns of X is a time series of length N. % % Returns: E - eigenfunction matrix. (LM by LM) % V - vector containing variances (unnormalized eigenvalues). % A - matrix of principal components. % C - lag-covariance matrix.
经验正交函数分析法（EOF）在matlab上的实现
热门推荐
03-05 5万+
原理与计算步骤 数据准备 数据集为福建省28个气象站1960-2013年的年际降水量数据，行为站点编号（数值格式），列为年份时间（数值格式）。部分数据如下： matlab代码 需要代码(包括REOF)的请根据文章最后面的信息联系我。 运行结果及分析 空间分布特征分析 前5个特征向量特征值的累积贡献率达到85.4%，但只有前两个特征根的误差范围不重叠通过...
matlab EOF程序
毛发浓密的猿工
11-02 5802
% eof第一模态图用变量eof的第一列，以此类推；相应的时间系数用pc第一行 clear clc x = [2 6 1 5 2;9 4 0 5 4;12 2 55 9 10;4 55 78 2 13]; %原始数据，列为站点或格点，行为时间序列上的值 x(1,:) = x(1,:) - mean(x(1,:)); x(2,:) = x(2,:) - mean(x(2,:)); %距平处理 c = x*x'/5; %协方差矩阵 [eof,eigenvalue] = ei
EOF分析matlab代码
10-21
气象统计EOF的MATLAB程序，可以直接使用，比较容易看懂的
EOF_EOG_经验正交函数分解_eof_
10-04
总的来说，"EOF_EOG_经验正交函数分解_eof_"是一个涉及EOG数据处理和EOF分解技术的主题，通过FORTRAN程序实现，有助于科学家和研究人员从复杂的眼球运动数据中提取关键信息，进一步理解和分析视觉行为。
EOF的matlab程序
01-15
EOF简单的matlab程序，实用易懂。
用于MATLAB的Climate Data Toolbox：用于分析和显示地球科学数据的一组底层函数和教程。-matlab开发
05-28
该工具箱包含一组标准的MATLAB函数，用于分析和显示气候数据。 这些函数计算效率高，易于使用，并且附带了许多教程，这些教程不仅描述了如何使用CDT函数，而且还提供了有关如何在地球科学过程的背景下解释结果的指南。 CDT徽标插图，作者：Adam S. Nelsen（ http://www.adamnelsen.com ）
eof-matlab
03-11
eof-analysis,diyicifa,kenengyouxiewenti,haiqingzhizheng
EOF 经验正交分解
03-23
eof 可以修改后再修改 matlab 中可以使用
EOF.zip_EOF；经验正交分解_matlab eof_正交分解_气象_海洋
09-24
EOF（Empirical Orthogonal Functions，经验正交函数）是一种在气象学、海洋学等领域广泛应用的数据分析方法，用于降维和模式识别。它通过对多维数据集进行正交分解，将复杂的数据结构简化为一系列主成分，这些...
eof_EOF分解_eof_经验分布_经验正交函数_
09-28
经验正交函数（EOF，Empirical Orthogonal Functions）是一种在气象学、海洋学、地球物理学以及其他领域广泛应用的数据分析方法，用于处理大型空间时间序列数据。它通过分解数据矩阵，将其转化为一组正交基，以便更...
MATLAB eof用法,经验正交函数分析法（EOF）在matlab上的实现
weixin_28871475的博客
03-16 2110
原理与计算步骤数据准备数据集为福建省28个气象站1960-2013年的年际降水量数据，行为站点编号，列为年份时间。部分数据如下：matlab代码clc;clear;data=xlsread('C:\Users\Administrator\Desktop\福建降水量.xls');[r1,c1]=size(data);x=data(2:end,2:end);[m,n]=size(x);fori=1:...
经验正交分解EOF详解及案例
同步于【自学气象人公众号】！
05-12 1万+
由于这些主分量中前几个就可能占有原空间点(即多个变量)总方差的很大部分，即通过EOF分解可以很容易地将原始要素场的变化信息浓缩在前几个主分量及其对应的空间函数上，所以，EOF分解常被用于气象要素场时空变化特征规律的提取和研究。因为所有空间模态(或时间系数)所对应的方差贡献率加起来等于 1，所以，当方差贡献率最高的几个空间模态(或时间系数)加起来得到的累计方差贡献率较大时(比如达到80%)，我们就会用它们几个来表示原始数据中所包含的主要特征。中第 t 时刻的空间场(共有m个空间点)，就可以表示为。
eof经验正交分解matlab
最新发布
08-29
EOF（Empirical Orthogonal Function）经验正交分解是一种常用的数据分析方法，用于提取数据中的主要模态和变异性。在MATLAB中，可以使用`eof`函数进行EOF分解。 `eof`函数的基本语法是： ``` [eof_maps, pc, expvar = eof(A) ``` 其中，`A`是一个数据矩阵，每行代表一个观测样本，每列代表一个变量。`eof_maps`是一个矩阵，每行代表一个EOF模态，每列对应于空间上的一个网格点，表示该模态在每个网格点上的空间分布。`pc`是一个矩阵，每行代表一个EOF模态，每列对应于时间上的一个时间步长，表示该模态在每个时间步长上的时间变化。`expvar`是一个向量，表示每个EOF模态解释的方差百分比。 `eof`函数还可以通过设置其他参数来进行更多的操作，比如设置`n`参数来指定要提取的主要EOF模态的数量，设置`mask`参数来指定要忽略的区域。 请注意，EOF分解是一种数据分析方法，适用于各种领域的数据，包括气象、海洋、地球物理等。它可以帮助我们了解数据的主要模态和变异性，并用于数据降维、模式识别等应用。<span class="em">1</span><span class="em">2</span><span class="em">3</span> #### 引用[.reference_title] - *1* *2* *3* [matlab经验正交分解函数EOF的实现—基于Climate Data Toolbox操作](https://blog.csdn.net/weixin_43637490/article/details/123418759)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_2"}}] [.reference_item style="max-width: 100%"] [ .reference_list ]
评论 24
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值