eof文档
的eof函数给出了固有模式的地图变化和相应的主成分时空数据的时间序列分析。这个函数是专门为3 d matricies海洋表面温度等数据的维度1和2维空间(例如,纬度和经度;经度和纬度;x和y,等等),第三维度代表不同的片或快照的数据。
内容
语法
eof_maps = eof (A) eof_maps = eof (A, n) eof_maps = eof(…,“面具”,面具)[eof_maps、pc、expvar] = eof (…)
描述
eof_maps = eof (A)计算所有模式的变化一个,在那里一个是一个三维矩阵的第一个二维空间,第三维是时间,数据被认为是等距的。输出eof_maps有相同的尺寸吗一个沿着第三维度,每个地图代表一个可变性重要性顺序模式。
eof_maps = eof (n),只计算第一n模式的变化。对于大型数据集,它的计算速度只计算模式你需要的数量。如果n没有指定,所有计算eof(每个时间片)。
eof_maps = eof(…,“面具”,面具)只执行EOF分析网格细胞为代表的逻辑面具的维度对应维度1和2的一个。此选项提供了防止解决不需要解决的事情,或者让你做分析独立于另一个地区。默认情况下,任何网格细胞一个含有任何nan掩盖了。
[eof_maps、pc、expvar] = eof (…)返回主分量时间序列个人电脑的每一行代表一个不同的模式从1到n和列对应于时间的步骤。例如,电脑(1:)时间序列的第一个varibility(主导)模式。第三输出expvar是方差的百分比用每种模式来解释。请参见下面的说明解释expvar。
一个简单的例子
这里的如何使用一个简单的例子eof函数。适当的EOF分析需要消除趋势和deseasoning数据计算EOF之前,这些步骤中描述下面的教程,但是现在让我们假装这个示例数据集已经准备好进行分析。加载示例数据,然后计算第一个EOF。
这里我策划第一EOF地图使用cmoceanδcolormap (Thyng et al ., 2016)“主”参数,以确保为中心是零。
%样本数据加载:负载PacOcean.mat%计算第一个EOF海面温度及其%主分量时间序列:(eofmap, pc) = eof (sst, 1);%画出第一个EOF地图:显示亮度图像(经度、纬度、eofmap);轴xy图像从%可选:使用cmocean colormap:cmocean (“δ”,“主”,0)
这是第一个EOF的SST数据集,但由于我们没有季节性周期,第一个EOF主要是季节性变化。作为证据,上面的模式与季节性周期有关,看看对应的主成分时间序列。
图绘制(t, pc)轴紧xlim ([datenum (1990年1月1日的)datenum (1995年1月1日的)))datetick (“x”,“keeplimits”)
我看起来很季节性。如果你喜欢情节异常时间序列共同双色阴影风格,使用异常函数可以在文件交换。
异常(t, pc)
清晰的
教程:从原始气候再分析数据来ENSO, PDO等等。
的eof函数提供了一个示例数据集PacOcean.mat哈德利中心的,这是一个downsampled子集的HadISST海面温度数据集。本教程的最后有一个部分描述如何我原始NetCDF数据导入到Matlab和过程我用来子集。如果你跟随本教程从上到下你应该能够应用EOF分析,任何类似的数据集。
如果你还没有加载示例数据集,现在装入并了解其内容通过检查变量的名称和大小:
负载PacOcean.mat谁
类属性名称大小字节纬度60 x1 480双经度55 x1 440海温60双x55x802 21172800双t 802 x1 6416双
我们有一个3 d风场矩阵维度对应纬度x经度的x。什么时间范围,什么是时间步,你问?让我们看看第一个和最后一个日期,和平均时间步骤:
结束datestr (t([1])的意思是(diff (t))
ans = 2×20字符数组的15 - 1月- 1950 12:00:00的15 - 10月- 2016 12:00:00 ans = 30.4370
海洋表面平均温度
好吧,这是月度数据,集中在每个月的15日,从1950年到2016年。得到一个数据集的样子,显示平均温度。我使用imagescn自动使NaN值透明,但您可以使用显示亮度图像代替。我也在使用cmocean热colormap (Thyng et al ., 2016):
图imagescn(经度,纬度,意味着(sst, 3));轴xy从cb = colorbar;ylabel (cb、“平均气温C{\保监会}”)cmocean热
全球变暖
全球变暖是真实的吗?的趋势函数让我们容易温度的线性趋势从1950年到2016年。一定要这一趋势乘以10 * 365.25将从每十年度每天度:
imagescn(经度、纬度、10 * 365.25 *趋势(海温,t, 3))轴xy从cb = colorbar;ylabel (cb、每十年的温度趋势{\保监会}C”)cmocean (“平衡”,“主”)
消除全球变暖的信号
全球变暖趋势是有趣,但EOF分析都是瑞士,没有长期趋势,所以我们必须删除的趋势detrend3:
海温= detrend3 (sst、t);
消除季节性周期
如果你再绘制温度趋势,你会发现这都是减少到零,或许几eps数值噪音。现在这是一个SST数据集,即使安东尼·瓦会批准的。
我们现在去趋势对海温数据集(也删除了的意思),但它仍然包含相当多的季节性变化,应该删除之前EOF分析季节性信号。因为我们不感兴趣快速消除季节性周期从月度数据集是确定平均SST在每个网格单元对于任何给定的月。从几个月开始对应于每一次介入t。我们不需要一天,所以我将波浪号(~)datevec输出,只保留一个月:
[~、月~]= datevec (t);
的月向量是相同的大小t,但只有12独特的价值观:
独特的(月)'
ans = 1 2 3 4 5 6 7 8 9 10 11 12
具体地说,这意味着每个时间步长与一个12个月的一年。有多少时间步长与1月吗?
总和(月= = 1)
ans = 67
有67 1月海温地图在完整的数据集,因为它是一个67年的记录。每个月,我们可以计算该月平均SST地图找到所有相关的时间步的指标。然后删除季节性信号减去平均海温地图的67年1月1月海温地图。这就是我的意思是:
% Preallocate每月3 d矩阵的方式:monthlymeans =南(长度(lat)、长度(朗),12);%计算所有地图对应的意思是每个月,和减去%的月度意味着sst数据集:为k = 1:12%月k指数:印第安纳州=月= = k;%对月平均海温凯西:monthlymeans (:,:, k) =意味着(sst(:,:,印第安纳州)3);%扣除每月的意思是:海温(:,:,印第安纳州)= bsxfun (@minus, sst(:,:,印第安纳州)monthlymeans (:,:, k));结束
我知道,bsxfun不直观。新版本的Matlab会让你减去2 dmonthlymeans从3 d矩阵风场数据集的隐式扩张,这意味着你可以用负号代替bsxfun,但如果你的Matlab版本早于2016你必须使用bsxfun方法如上所示。
季节性周期的gif
我们只是季节性周期,但它是在一个循环中,我们没有看到什么信息是我们的风场数据集。让我们做一个gif使用我gif功能:
%第一帧:图h = imagescn(经度、纬度、monthlymeans (:,: 1));xy轴图像从cb = colorbar;-10 caxis ([10]) cmocean平衡
标题(datestr (datenum(1, 1, 1), '嗯'))
%写第一帧:gif (gcf“SeasonalTemperatureAnomalies.gif”、“框架”,“延迟时间”,0.2)
%循环通过所有其他帧:k = 2:12 %颜色数据更新:组(h, cdata, monthlymeans (:,:, k))
%更新标题:标题(datestr (datenum (1 k 1), '嗯'))
%写这个月的框架:gif
所以现在我们风场数据集已经去趋势,意味着移除和季节性周期移除。剩下的风场异常——事情改变,但不是长期趋势或短期的年度周期。这是其余的方差风场异常数据集:
图imagescn(经度、纬度、var (sst, [], 3));轴xy从colorbar标题(温度的差异)colormap(飞机)%飞机是不可原谅的,除了重建旧的情节caxis ([0 1])
和上面的地图线路与图2的一个很好Messie和查韦斯(2011)告诉我们,我们在正确的轨道上。
计算eof
EOF分析不仅告诉我们事情各有不同,但多长时间,什么地区倾向于彼此在一起或的阶段不同。与我们的去趋势,deseasoned风场数据集,这个简单的EOF分析eof功能:
[eof_maps, pc, expv] = eof (sst);%画第一个模式:图显示亮度图像(经度、纬度、eof_maps(:,: 1)轴xy图像cmocean (“卷”,“主”)标题“第一个EOF模式!”
特征向量分析有一个有趣的行为可以产生EOF地图是积极的还是消极的,和解决方案可以不同的每次使用相同的输入。金宝搏官方网站积极的和消极的解决方案同样有效,认为鼓面的振金宝搏官方网站动模式,一些地区的鼓头上升下降,而其他区域,然后开关,同样海温变化的特征值的解决方案可能是积极的还是消极的。唯一重要的是,当我们重构时间序列的EOF的解决方案,我们每一个EOF映射乘以相应的主成分(个人电脑)。
我写的eof函数产生一致的结果每次运行相同的数据,但别担心如果解决不匹配的符号的符号别人的结果,如果他们使用不同的程序计算eof——这就意味着他们的程序反号的解决方案,这是非常好。
正如EOF地图可以有积极或消极的解决方案,都是同样有效的,有一些灵活性EOF地图的大小是如何显金宝搏官方网站示的。你可以把一个EOF地图的大小乘以任何你想要的价值,只要你相应的主成分时间序列除以相同的值。让我们看看前三个模式变化的时间序列:
图绘制(t,电脑(1:3,:))盒子从轴紧datetick (“x”,“keeplimits”)传说(“pc1”,“pc2”,“生物”)
可选的主成分和EOF地图缩放
这些主成分时间序列是好方法,但是有些人更喜欢每个时间序列跨度规模所需的范围。查看图5的Messie和查韦斯(2011)看来,他们选择规模每个主分量时间序列,使其涵盖范围1到1。让我们做同样的事情,每个主分量时间序列除以它的最大价值和不要忘记乘以相应的EOF映射相同的值:
为k = 1:尺寸(pc, 1)%的时间序列中的最大值每个主要组件:[maxval,印第安纳州]= max (abs (pc (k,:)));%将时间序列的最大值:电脑(k,:) = pc (k,:) / maxval;%乘以相应的EOF地图:eof_maps (:,:, k) = eof_maps (:,:, k) * maxval;结束
厄尔尼诺南方涛动(ENSO)时间序列
第一个模式去趋势,deseasoned太平洋与ENSO assoiciated。我们可以再次画出时间序列作为一个简单的情节,但异常的情节往往填充。让我们使用异常情节第一个模式,乘以1与图5的迹象Messie和查韦斯(2011)。
图(“pos”,100 100 600 250)异常(t - pc (1:)“topcolor”,[1 .35点点),“bottomcolor”,(。51算1])%第一主成分是enso盒子从轴紧datetick (“x”,“keeplimits”)文本((724316 729713 736290)。95 .99结果),“厄尔尼诺”,“水平的”,“中心”)
果然,一些最强的厄尔尼诺事件记录发生在1982 - 1983,1997 - 1998,2014 - 2016。
ENSO在频域
有时我们听到厄尔尼诺现象的特征频率每隔五年,或五到七年,有时你听到每两到七年。很难看到在时间序列,所以我们把第一个主成分在频域plotpsd,指定每年12个样品的采样频率,标注在日志x轴,在单位的x值λ(年)而不是频率:
图plotpsd (pc(1:), 12日“计算lnx”,“λ”)包含“周期(年)”集(gca),“xtick”33,[1:7])
如您所见,ENSO信号没有大幅共振频率定义,但整个两到七年的能量范围。我也贴上33年周期性因为奈奎斯特的这个特定的数据集——任何能量比尼奎斯特更长一段(或接近)应该被认为是垃圾。
可变性的地图
eof不仅仅是对时间序列——他们对空间的变化模式。每个模式都有一个特征变化模式就像一个鼓的不同振动模式。在任何给定的时间,不同的模式可以总结为创建一个总温度异常的照片。的正交经验正交函数的一部分意味着每个模式倾向于做自己的事情,独立于其他模式。让我们看看第一个6模式的再现图4Messie和查韦斯(2011)。我的一些模式乘以- 1因为我想匹配他们的迹象,记住,我们可以这样做。
s = [1 1 1 1 1 1];%(签署乘数匹配Messie和查韦斯2011)图(“pos”,100 100 500 700)为k = 1:6次要情节(3 2 k) imagescn(经度、纬度、eof_maps (:,:, k) * (k));轴xy从标题([“模式”num2str (k),“(”num2str (expv (k),' % 0.1 f '),“%)”2)caxis ([2])结束colormap飞机
百分比差异解释为每种模式匹配Messie和查韦斯因为我们用更短的时间序列相比,我们也使用一个空间世界数据的子集。尽管如此,普遍认同模式。
的飞机colormap不是完全相同的一个使用Messie查韦斯,这解释了为什么一些模式高于Messie和查韦斯看起来略有不同。但由于我们谈论colormaps,彩虹是非常不擅长代表数值数据(Thyng et al, 2016年)。科学也是一种耻辱,我们不能完全复制情节不知道什么颜色是用于发布的版本,但是我跑题了…
鉴于这些地图代表异常,他们应该由发散colormap表示出同等重量的0。设的colormaps次要情节在这个图更平衡:
colormap (gcf cmocean (“平衡”))
制作电影的海温变化eof
在任何给定的时间,海面温度异常与ENSO相关的快照可以通过绘制的地图模式1所示,乘以相应的主成分(向量电脑(1:))。同样,你可以得到一幅全球海面温度异常在给定的时间总结所有的EOF地图,每个乘以相应的主成分。这样我们可以构建一个越来越多的海温异常在我们完整的电影包括越来越多的更多的模式变化。同样的概念,可以排除模式过滤掉不受欢迎的信号,或者我们可以使用的头几个模式作为一种过滤噪音。让我们制作一个电影的前三个模式,从1990年到2005年:
%指数电影的开始和结束日期:startind =找到(t > = datenum(1990年1月1日),1,'第一次');endind =找到(t < = datenum(1999年12月31日),1,'最后');
%的地图海温异常在开始前三个模式:地图= eof_maps (:,: 1) * pc (1, startind) +…% 1模式,1990年1月eof_maps (:,: 2) * pc (2, startind) +…%模式2,1990年1月eof_maps (:,:, 3) * pc (3 startind);%模式3,1990年1月
%创建第一帧的影片:图h = imagescn(经度、纬度、地图);xy轴图像从cb = colorbar;caxis (2 [2]) cmocean平衡标题(datestr (t (startind),“yyyy”)) %初始化gif
gif (“SSTs_1990s.gif”、“框架”,gcf)
%循环剩余的饥饿:k = (startind + 1): endind
% = eof_maps更新日期k地图映射(:,:1)* pc (k) +…%模式1 eof_maps (:: 2) * pc (2 k) +…%模式2 eof_maps (:,:, 3) * pc (3 k);%模式3
组(h, cdata,地图)caxis(2[2])标题(datestr (t (k),“yyyy”))
gif %将这个框架添加到gif
你可能注意到的第一件事是,1990年代的海温异常时间序列是由ENSO,并检查出1997 - 1998信号!难怪这是这样的一个热门话题新闻那一年。尽管如此,重要的是要记住,上面的电影并不是一个完整的重建对海温异常,而是只有前三种模式,占
sum (expv (1:3))
ans = 51.8850
…超过一半的总方差的SST数据集。重建的绝对温度场而不是异常的前三个模式,需要包括所有的EOF地图,你还需要添加在平均SST地图,趋势和季节性周期。
关于解释的方差解释道。
上面,我们解决了所有802模式包含802时间步的时间序列。在一起,这些模式解释的总方差的100%风场数据集。这意味着我们可以添加所有这802模式,乘以他们的主成分,完全重构时间序列没有任何丢失的信息。看,解释方差之和……
总和(expv)
ans = 100.0000
…是100%。这个数据集和真实数据集的其他自然模式的可变性,最初几个模式将大部分的信息和其他所有应该被忽略。有时是深刻的解释方差的函数模式数字,也很常见,累计方差解释道。让我们同时绘制:
图次要情节(2,1,1)情节(expv)轴([0 802 0 37])从ylabel方差解释说每模式”次要情节(2,1,2)情节(cumsum (expv))轴ylabel ([0 802 0 100])“累积方差解释”包含“模式数字”盒子从
上面的图表明,前三个模式占一半以上的去趋势方差,deseasoned SST数据集,一切只是增加了小改进,让我们一点点接近重建完整的原始数据集。
一般的经验法则,我倾向于认为任何过去的第三或第五模式有点罗夏测试可能欺骗你认为有意义的信号,但小心试图解释这些n秩序模式——甚至在这个简单的数据集,有802模式,这意味着有可能至少有一个模式,似乎能解释你可能有任何疯狂的理论。换句话说,留心的潜力p-hacking当解释eof。
大部分的802年模式我们求出了在上面的示例中,每个贡献少量海温变化的整体,身体可能是毫无意义的,所以没有意义上使用太多电脑所能解出所有这些模式。如果您的数据集很大,它将更快,使用更少的内存来解决你所需要的模式。只有一个小抓……
%解决所有模式:[eof_maps_all, pc_all expv_all] = eof (sst);%求出第一个10模式:[eof_maps_10, pc_10 expv_10] = eof(海温,10);
我们已经解决了第一个10模式,我们看到EOF地图和主成分与解决方案获得的时间序列在协议解决所有802个模式:金宝搏官方网站
图次要情节(1、2、1)显示亮度图像(eof_maps_all (:,: 1) -eof_maps_10 (:,: 1) colorbar cmocean (“平衡”,“主”次要情节(1、2、2)情节(pc_all (: 1) -pc_10(1:))轴紧
上面的两个图显示纯数值噪音。有一些明显的结构,但看看左边的颜色轴的凯尔(10 ^ -17)和右边的y轴(10 ^ -14)。的舍入误差,结果从Matlab数字化数字作品,这是常见的舍入误差的大小来追踪数字的大小正圆。
虽然有一点的数值噪音,解决方案在解决所有802模式一样的解决方案在解决10模式。金宝搏官方网站太好了。直到我们的解释方差:
expv_10图绘制(expv_all (1:10),“o”)举行在情节(50 [0],50 [0])%直线(两个值相等)轴平等的包含“解释方差(802解决方案”)金宝搏官方网站ylabel“解释方差(10解决方案”)金宝搏官方网站
解释的方差不匹配,因为expv提供了一个衡量比例的方差解决在给定的模式运行。也就是说,如果我们解出所有802模式,方差之和解释为前10模式将不等于100%:
sum (expv_all (1:10))
ans = 73.4559
然而,由于我们只解决了前10模式expv_10,其金额等于100%
sum (expv_10 (1:10))
ans = 100.0000
因为那些10模式解释100%的方差与这10模式。有些循环逻辑,我知道,所以有什么好解释的方差的值如果你只解决了几个模式?也许他们并不是那么好。你可以做的最好的事情就是使用一个长时间序列,求出所有的模式,然后看看解释方差的渐近线。
如果你只能解决几个模式由于计算的局限性,你仍然可以的你是否已经完整的图片通过绘制的累积和解释方差从任何方式解出:
图绘制(cumsum (expv_10),“啊——”)包含“模式数字”ylabel“累积解释方差之和(%)”
当然阴谋的最终价值是100%(我们知道意味着100%的大约73.5%这个数据集),所以100%的价值并没有告诉我们自己,虽然我们可以清楚地看到它不到的渐近线之前到100%,表明有更多的故事……
对不起我没有整理者现在解决这个问题,但希望这个小讲解员给一个更好的理解eof函数可以,不能告诉你。
我得到了样本数据如何
上面所示的示例数据集来自HadISST哈德利中心找到在这里(雷纳et al ., 2003),在超过200 MB。如果你想执行同样的在世界不同地区的分析,您可以下载HadISST_sst。数控数据集,将其导入Matlab。将采样数据集或构造子集是你:
%加载完整的SST数据集:lat =双(ncread (“HadISST_sst.nc”、“纬度”));朗=双(ncread (“HadISST_sst.nc”、“经度”));双(t = ncread (“HadISST_sst.nc”,“时间”)+ datenum(1870年,1,0));海温= ncread (“HadISST_sst.nc”、“场”);
%季度样本数据集的大小,我对每一个网格点大致downsample:海温(sst(1:2:结束,1:2:最终,);lat =纬度(1:2:结束);朗=经度(1:2:结束);
%,以进一步降低尺寸,我剪的背阔肌和经度,保持只有1950后数据:行=朗< -70;朗=朗(行);关口= lat > = -60 & lat < = 60;lat =纬度(关口);* = t > = datenum(1950年1月1日);t = t(次);海温(sst(行、关口倍);海温(sst < -50) =南;
%我发现它更容易重新排列纬度x经度x:海温=排列(风场(2 1 3));
%保存样本数据:保存(“PacOcean.mat”、“纬度”,“朗”、“t”,“场”)
引用
Monique Messie,旧金山查韦斯。“全球海洋表面温度变化模式与区域气候指标。”Journal of Climate 24.16 (2011): 4314-4331.jcli3941.1 doi: 10.1175/2011。
雷纳:A。,Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., Kent, E. C., Kaplan, A. (2003). Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res.Vol. 108, No. D14, 4407jd002670 doi: 10.1029/2002。
Thyng K.M.,C.A. Greene, R.D. Hetland, H.M. Zimmerle, and S.F. DiMarco. 2016. True colors of oceanography: Guidelines for effective and accurate colormap selection. Oceanography 29(3):9-13,doi: 10.5670 / oceanog.2016.66。
作者信息
的eof函数是由乍得a。格林德克萨斯大学的地球物理研究所(UTIG) 2017年1月,但靠Guillame迷宫的caleof函数从他PCATool的贡献。本教程是乍得格林的帮助下写的Kaustubh Thirumalai。