|
@@ -0,0 +1,1982 @@
|
|
|
|
+CFU_PCA<-function(data, cl = 0.98, ofiles = FALSE) {
|
|
|
|
+#
|
|
|
|
+# Computes a Principal Component Analysis (PCA/EOF) of data.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns an array of EOFs, an array of PCs, a vector of FVARs (fractions of explained variance),
|
|
|
|
+# an array of CORRs (correlation maps), and the threshold for the correlation coefficient based
|
|
|
|
+# on a two-tailed Student's t-test.
|
|
|
|
+# If applied, also returns ascii files for PCs and FVARs.
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_PCA(data, sl = 0.98, ofiles = FALSE)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# data: Array containing the anomalies field for the analysis
|
|
|
|
+#
|
|
|
|
+# cl: Value of the confidence level for the correlation map.
|
|
|
|
+# It is based on a two-tailed Student's t-test.
|
|
|
|
+# Default is cl=0.98 (98% confidence level; 2% significance level).
|
|
|
|
+#
|
|
|
|
+# ofiles: Logical.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# $PCs: Array of principal components (nt, nm)
|
|
|
|
+# $FVARs: Vector of fractions of explained variance (nm)
|
|
|
|
+# $EOFs: Array of empirical orthogonal functions (c(dim), nm)
|
|
|
|
+# $CORRs: Array of correlation maps associated with EOFs (c(dim), nm)
|
|
|
|
+# $Rtest: Value of the statistically significant corr.coefficient at cl
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 17 September 2010
|
|
|
|
+# Modified by jgarcia@ic3.cat (29 August 2011): The calculation is now done by a
|
|
|
|
+# singular value decomposition of the data matrix (prcomp), not by using 'eigen'
|
|
|
|
+# on the covariance matrix (princomp).
|
|
|
|
+#
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+d=dim(data)
|
|
|
|
+if (length(d) == 3){
|
|
|
|
+ nlon=d[1]; nlat=d[2]; ns=nlon*nlat; nt=d[3]
|
|
|
|
+} else {
|
|
|
|
+ nlon=1; nlat=d[1]; ns=nlon*nlat; nt=d[2]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+M=array(dim=c(ns,nt))
|
|
|
|
+for (i in 1:nlat){
|
|
|
|
+ M[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]=data[,i,]
|
|
|
|
+}
|
|
|
|
+not=is.na(M[,1])
|
|
|
|
+yes=array(dim=c(ns))
|
|
|
|
+counter=0
|
|
|
|
+for (i in 1:ns){
|
|
|
|
+ if (not[i]=='FALSE'){counter=counter+1; yes[counter]=i}
|
|
|
|
+}
|
|
|
|
+yes=yes[1:counter]
|
|
|
|
+Myes=M[yes,] #(ns_with_data,nt)
|
|
|
|
+
|
|
|
|
+#covariance matrix
|
|
|
|
+
|
|
|
|
+#covM=(Myes%*%t(Myes))/(nt-1)
|
|
|
|
+
|
|
|
|
+#PCA analysis
|
|
|
|
+
|
|
|
|
+#PCA=princomp(covmat=covM)
|
|
|
|
+PCA=prcomp(t(Myes))
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+nm=10 #number of modes to be retained (aprox 80% of FVAR)
|
|
|
|
+
|
|
|
|
+std=PCA$sdev
|
|
|
|
+lambda=std*std; fvar=lambda/sum(lambda)
|
|
|
|
+fvar=round(fvar,digits=4); fvar=fvar[1:nm]*100
|
|
|
|
+
|
|
|
|
+#U=PCA$loadings #(ns_with_data,ns_with_data)
|
|
|
|
+U=PCA$rotation
|
|
|
|
+Un=U[,1:nm]
|
|
|
|
+
|
|
|
|
+PC=t(Myes)%*%Un #(nt,nm)
|
|
|
|
+PCstd=array(dim=c(nt,nm))
|
|
|
|
+for (i in 1:nm){
|
|
|
|
+ S=sd(PC[,i]); PCstd[,i]=PC[,i]/S
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+EOFyes=(Myes%*%PCstd)/nt #(ns_with_data,nm)
|
|
|
|
+eof=array(dim=c(ns,nm))
|
|
|
|
+eof[yes,]=EOFyes
|
|
|
|
+EOF=array(dim=c(nlon,nlat,nm))
|
|
|
|
+for (i in 1:nlat){
|
|
|
|
+ EOF[,i,]=eof[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CORRyes=array(dim=c(counter,nm))
|
|
|
|
+for (i in 1:nm){
|
|
|
|
+ for (j in 1:counter){
|
|
|
|
+ CORRyes[j,i]=cor(Myes[j,],PCstd[,i])
|
|
|
|
+ }
|
|
|
|
+}
|
|
|
|
+corr=array(dim=c(ns,nm))
|
|
|
|
+corr[yes,]=CORRyes
|
|
|
|
+CORR=array(dim=c(nlon,nlat,nm))
|
|
|
|
+for (i in 1:nlat){
|
|
|
|
+ CORR[,i,]=corr[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#t-Test
|
|
|
|
+
|
|
|
|
+cl=cl; cl2=(1-cl)/2
|
|
|
|
+t=qt(cl+cl2, nt-2); R=sqrt((t*t)/((t*t)+nt-2))
|
|
|
|
+
|
|
|
|
+#options
|
|
|
|
+
|
|
|
|
+#if(ofiles) {
|
|
|
|
+# save(PCstd, ascii=TRUE, file="$plotsdir/PCs_of_data.dat")
|
|
|
|
+# save(fvar, ascii=TRUE, file="$plotsdir/FVARs_of_data.dat")
|
|
|
|
+#}
|
|
|
|
+
|
|
|
|
+invisible(list(PCs=PCstd,FVARs=fvar,EOFs=EOF,CORRs=CORR,Rtest=R))
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_MCA<-function(datax, datay, cl = 0.98, ofiles = FALSE) {
|
|
|
|
+#
|
|
|
|
+# Computes a Maximum Covariance Analysis (MCA) between datax and datay.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns a vector of squared covariance fraction (SCFs) explained by each pair of
|
|
|
|
+# covariability modes, a vector of correlation coefficient (RUVs) between expansion
|
|
|
|
+# coefficients (ECs) that measures their linear relationship, and a set of regression (MCAs)
|
|
|
|
+# and correlation (CORRs) maps associated with the covariability modes. Note that both,
|
|
|
|
+# MCAs and CORRs, are 'homogeneous' patterns obtained as regression/correlation between
|
|
|
|
+# each field (predictor, predictand) and its expansion coefficient.
|
|
|
|
+# The threshold for the correlation coefficient, based on a two-tailed Student's t-test,
|
|
|
|
+# is also produced.
|
|
|
|
+# If applied, also returns ascii files for SCFs, RUVs, and ECs.
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_MCA(datax, datay, sl = 0.98, ofiles = FALSE)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# datax: Array containing the anomalies field for the predictor
|
|
|
|
+#
|
|
|
|
+# datay: Array containing the anomalies field for the predictand
|
|
|
|
+#
|
|
|
|
+# cl: Value of the confidence level for the correlation map.
|
|
|
|
+# It is based on a two-tailed Student's t-test.
|
|
|
|
+# Default is cl=0.98 (98% confidence level; 2% significance level).
|
|
|
|
+#
|
|
|
|
+# ofiles: Logical.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# $SCFs: Vector of squared covariance fractions (nm)
|
|
|
|
+# $RUVs: Vector of correlations between expansion coefficients (nm)
|
|
|
|
+# $ECs_U: Array of expansion coefficients of predictor field (nt, nm)
|
|
|
|
+# $MCAs_U: Array of covariability patterns of predictor field (c(dim), nm)
|
|
|
|
+# $CORRs_U:Array of corr. maps associated with predictor field (c(dim), nm)
|
|
|
|
+# $ECs_V: Array of expansion coefficients of predictand field (nt, nm)
|
|
|
|
+# $MCAs_V: Array of covariability patterns of predictand field (c(dim), nm)
|
|
|
|
+# $CORRs_V:Array of corr. maps associated with predictand field (c(dim), nm)
|
|
|
|
+# $Rtest: Value of the statitically significant corr. coefficient at cl
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 24 September 2010
|
|
|
|
+#
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+dx=dim(datax)
|
|
|
|
+if (length(dx) == 3){
|
|
|
|
+ nlonx=dx[1]; nlatx=dx[2]; nsx=nlonx*nlatx; nt=dx[3]
|
|
|
|
+} else {
|
|
|
|
+ nlonx=1; nlatx=dx[1]; nsx=nlonx*nlatx; nt=dx[2]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+X=array(dim=c(nsx,nt))
|
|
|
|
+for (i in 1:nlatx){
|
|
|
|
+ X[(nlonx*(i-1)+1):(nlonx*(i-1)+nlonx),]=datax[,i,]
|
|
|
|
+}
|
|
|
|
+not=is.na(X[,1])
|
|
|
|
+yesx=array(dim=c(nsx))
|
|
|
|
+counterx=0
|
|
|
|
+for (i in 1:nsx){
|
|
|
|
+ if (not[i]=='FALSE'){counterx=counterx+1; yesx[counterx]=i}
|
|
|
|
+}
|
|
|
|
+yesx=yesx[1:counterx]
|
|
|
|
+Xyes=X[yesx,] #(nsx_with_data,nt)
|
|
|
|
+
|
|
|
|
+dy=dim(datay)
|
|
|
|
+if (length(dy) == 3){
|
|
|
|
+ nlony=dy[1]; nlaty=dy[2]; nsy=nlony*nlaty; nt=dy[3]
|
|
|
|
+} else {
|
|
|
|
+ nlony=1; nlaty=dy[1]; nsy=nlony*nlaty; nt=dy[2]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+Y=array(dim=c(nsy,nt))
|
|
|
|
+for (i in 1:nlaty){
|
|
|
|
+ Y[(nlony*(i-1)+1):(nlony*(i-1)+nlony),]=datay[,i,]
|
|
|
|
+}
|
|
|
|
+not=is.na(Y[,1])
|
|
|
|
+yesy=array(dim=c(nsy))
|
|
|
|
+countery=0
|
|
|
|
+for (i in 1:nsy){
|
|
|
|
+ if (not[i]=='FALSE'){countery=countery+1; yesy[countery]=i}
|
|
|
|
+}
|
|
|
|
+yesy=yesy[1:countery]
|
|
|
|
+Yyes=Y[yesy,] #(nsy_with_data,nt)
|
|
|
|
+
|
|
|
|
+#covariance matrix
|
|
|
|
+
|
|
|
|
+covXY=(Xyes%*%t(Yyes))/(nt-1)
|
|
|
|
+
|
|
|
|
+#MCA analysis
|
|
|
|
+
|
|
|
|
+MCA=svd(covXY)
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+nm=10 #number of modes to be retained
|
|
|
|
+
|
|
|
|
+SV=MCA$d #singular values
|
|
|
|
+SV2=SV*SV; S=sum(SV2)
|
|
|
|
+SCF=SV2/S #squared covariance fraction
|
|
|
|
+SCF=SCF[1:nm]*100
|
|
|
|
+
|
|
|
|
+Ux=MCA$u
|
|
|
|
+Vy=MCA$v
|
|
|
|
+Un=Ux[,1:nm]; Vn=Vy[,1:nm]
|
|
|
|
+ECu=t(Xyes)%*%Un; ECv=t(Yyes)%*%Vn #(nt,nm)
|
|
|
|
+ECu_std=array(dim=c(nt,nm))
|
|
|
|
+ECv_std=array(dim=c(nt,nm))
|
|
|
|
+for (i in 1:nm){
|
|
|
|
+ Su=sd(ECu[,i]); ECu_std[,i]=ECu[,i]/Su
|
|
|
|
+ Sv=sd(ECv[,i]); ECv_std[,i]=ECv[,i]/Sv
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+RUV=array(dim=c(1,nm))
|
|
|
|
+for (i in 1:nm){
|
|
|
|
+ RUV[i]=cor(ECu_std[,i],ECv_std[,i])
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+Uyes=(Xyes%*%ECu_std)/nt #(ns_with_data,nm)
|
|
|
|
+Vyes=(Yyes%*%ECv_std)/nt
|
|
|
|
+mcaU=array(dim=c(nsx,nm)); mcaU[yesx,]=Uyes
|
|
|
|
+mcaV=array(dim=c(nsy,nm)); mcaV[yesy,]=Vyes
|
|
|
|
+MCAU=array(dim=c(nlonx,nlatx,nm))
|
|
|
|
+MCAV=array(dim=c(nlony,nlaty,nm))
|
|
|
|
+for (i in 1:nlatx){
|
|
|
|
+ MCAU[,i,]=mcaU[(nlonx*(i-1)+1):(nlonx*(i-1)+nlonx),]
|
|
|
|
+}
|
|
|
|
+for (i in 1:nlaty){
|
|
|
|
+ MCAV[,i,]=mcaV[(nlony*(i-1)+1):(nlony*(i-1)+nlony),]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CORRUyes=array(dim=c(counterx,nm))
|
|
|
|
+CORRVyes=array(dim=c(countery,nm))
|
|
|
|
+for (i in 1:nm){
|
|
|
|
+ for (j in 1:counterx){
|
|
|
|
+ CORRUyes[j,i]=cor(Xyes[j,],ECu_std[,i])
|
|
|
|
+ }
|
|
|
|
+ for (j in 1:countery){
|
|
|
|
+ CORRVyes[j,i]=cor(Yyes[j,],ECv_std[,i])
|
|
|
|
+ }
|
|
|
|
+}
|
|
|
|
+corrU=array(dim=c(nsx,nm)); corrU[yesx,]=CORRUyes
|
|
|
|
+corrV=array(dim=c(nsy,nm)); corrV[yesy,]=CORRVyes
|
|
|
|
+CORRU=array(dim=c(nlonx,nlatx,nm))
|
|
|
|
+CORRV=array(dim=c(nlony,nlaty,nm))
|
|
|
|
+for (i in 1:nlatx){
|
|
|
|
+ CORRU[,i,]=corrU[(nlonx*(i-1)+1):(nlonx*(i-1)+nlonx),]
|
|
|
|
+}
|
|
|
|
+for (i in 1:nlaty){
|
|
|
|
+ CORRV[,i,]=corrV[(nlony*(i-1)+1):(nlony*(i-1)+nlony),]
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#t-Test
|
|
|
|
+
|
|
|
|
+cl=cl; cl2=(1-cl)/2
|
|
|
|
+t=qt(cl+cl2, nt-2); R=sqrt((t*t)/((t*t)+nt-2))
|
|
|
|
+
|
|
|
|
+#options
|
|
|
|
+
|
|
|
|
+#if(ofiles) {
|
|
|
|
+# save(SCF, ascii=TRUE, file="$plotsdir/SCFs_of_data.dat")
|
|
|
|
+# save(RUV, ascii=TRUE, file="$plotsdir/RUVs_of_data.dat")
|
|
|
|
+# save(ECu_std, ascii=TRUE, file="$plotsdir/ECsU_of_data.dat")
|
|
|
|
+# save(ECv_std, ascii=TRUE, file="$plotsdir/ECsV_of_data.dat")
|
|
|
|
+#}
|
|
|
|
+
|
|
|
|
+invisible(list(SCFs=SCF,RUVs=RUV,ECs_U=ECu_std,MCAs_U=MCAU,CORRs_U=CORRU,ECs_V=ECv_std,MCAs_V=MCAV,CORRs_V=CORRV,Rtest=R))
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_lonlat2ns<-function(data, miss.val = FALSE){
|
|
|
|
+#
|
|
|
|
+# Reshape a 3D array into a 2D array handling only spatial dimensions
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns a 2D array of 3D data placing the spatial dimensions into a vector according
|
|
|
|
+# to GrADS protocol, and keeps the temporal dimension with no change.
|
|
|
|
+# It localizes missing values and returns an (spatial) index with no-missing values positions.
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_latlon2ns(data, miss.val = FALSE)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# data: 3-dimensional array, e.g. (longitude, latitude, time)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# matrix: 2-dimensional array with 'ns' spatial dimensions, e.g. (ns, time)
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 30 September 2010
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+d=dim(data)
|
|
|
|
+nlon=d[1]; nlat=d[2]; ns=nlon*nlat; nt=d[3]
|
|
|
|
+
|
|
|
|
+M=array(dim=c(ns,nt))
|
|
|
|
+for (i in 1:nlat){
|
|
|
|
+ M[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]=data[,i,]
|
|
|
|
+}
|
|
|
|
+not=is.na(M[,1])
|
|
|
|
+yes=array(dim=c(ns))
|
|
|
|
+counter=0
|
|
|
|
+for (i in 1:ns){
|
|
|
|
+ if (not[i]=='FALSE'){counter=counter+1; yes[counter]=i}
|
|
|
|
+}
|
|
|
|
+yes=yes[1:counter]
|
|
|
|
+Myes=M[yes,] #(ns_with_data,nt)
|
|
|
|
+
|
|
|
|
+#printed in screen
|
|
|
|
+
|
|
|
|
+print(ns)
|
|
|
|
+print(counter)
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+if(miss.val) { Myes <- M }
|
|
|
|
+
|
|
|
|
+invisible(list(matrix=Myes,nt=nt,ns_original=ns,ns_noNaN=counter,index_noNaN=yes))
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_ns2lonlat<-function(data, nlon=1, nlat=1, yes=NULL){
|
|
|
|
+#
|
|
|
|
+# Reshape a 2D array into a 3D array handling only spatial dimensions
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns a 2D array of 3D data placing the spatial dimensions into a vector according
|
|
|
|
+# to GrADS protocol, and keeps the temporal dimension with no change.
|
|
|
|
+# It localizes missing values and returns an (spatial) index with no-missing values positions.
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_latlon2ns(data, miss.val = FALSE)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# data: 3-dimensional array, e.g. (longitude, latitude, time)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# matrix: 2-dimensional array with 'ns' spatial dimensions, e.g. (ns, time)
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 06 October 2010
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+nlon=max(abs(nlon), 1); nlat=max(abs(nlat), 1); ns=nlon*nlat
|
|
|
|
+d=dim(data)
|
|
|
|
+nt=d[2]
|
|
|
|
+
|
|
|
|
+miss.val=is.null(yes)
|
|
|
|
+if(miss.val=='TRUE'){
|
|
|
|
+ M=array(dim=c(nlon,nlat,nt))
|
|
|
|
+ for (i in 1:nlat){
|
|
|
|
+ M[,i,]=data[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]
|
|
|
|
+ }
|
|
|
|
+} else {
|
|
|
|
+ Myes=array(dim=c(ns,nt)); Myes[yes,]=data
|
|
|
|
+ M=array(dim=c(nlon,nlat,nt))
|
|
|
|
+ for (i in 1:nlat){
|
|
|
|
+ M[,i,]=Myes[(nlon*(i-1)+1):(nlon*(i-1)+nlon),]
|
|
|
|
+ }
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+matrix=M
|
|
|
|
+matrix
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_detrend<-function(vector){
|
|
|
|
+#
|
|
|
|
+# Remove the linear trend in a time-series
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns a vector with the same length as the original (nt) but in which
|
|
|
|
+# the linear trend has been removed.
|
|
|
|
+# No change is produced if vector is a vector of missing values.
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_detrend(vector)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# vector: 1-dimensional array, time (nt dimension)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# 1-dimensional array without the linear trend (nt)
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 05 October 2010
|
|
|
|
+# Modified by jgarcia & vguemas@ic3.cat (27 April 2011) : to convert output into a vector (before, a list)
|
|
|
|
+
|
|
|
|
+nt=length(vector)
|
|
|
|
+
|
|
|
|
+if (is.na(vector[1])=='TRUE'){
|
|
|
|
+ print("vector contains missing values")
|
|
|
|
+ detrend=array(dim=nt)
|
|
|
|
+} else {
|
|
|
|
+detrend=resid(lm(vector ~ seq(vector)))
|
|
|
|
+}
|
|
|
|
+detrend=as.vector(detrend)
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+detrend
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_testRR<-function(R1, N1, R2, N2, cl = 0.95){
|
|
|
|
+#
|
|
|
|
+# Testing the equality of correlation coefficients from independent samples
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Statistical test for the equality of correlation coefficients (R1, R2)
|
|
|
|
+# from two independent samples. These samples should belong to two populations,
|
|
|
|
+# e.g. dynamical prediction (GCM) and statistical prediction; and the hypothesis
|
|
|
|
+# will test the equality of the population correlation (bilateral test).
|
|
|
|
+#
|
|
|
|
+# H0: R1==R2 vs. H1: R1=/=R2
|
|
|
|
+#
|
|
|
|
+# The statistic Z follows a N(0,1) distribution; if |Z| > Z_alfa/2, H0 is rejected.
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_testRR(R1, N1, R2, N2, cl)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# R1: number, correlation of samples #1
|
|
|
|
+#
|
|
|
|
+# N1: number, nº of elements in the sample (length #1)
|
|
|
|
+#
|
|
|
|
+# R2: number, correlation of samples #2
|
|
|
|
+#
|
|
|
|
+# N2: number, nº of elements in the sample (length #2)
|
|
|
|
+#
|
|
|
|
+# cl: Value of the confidence level.
|
|
|
|
+# Default is cl=0.95 (95% confidence level; 5% significance level).
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# value: binary output (0,1) / 0 = H0 accepted, 1 = H0 rejected
|
|
|
|
+#
|
|
|
|
+# p_value: area under the N(0,1) density curve to the right of endpoint |Z|
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 02 December 2010
|
|
|
|
+# Modified by jgarcia & lrodrigues@ic3.cat (13 December 2011) : Correcting p-values
|
|
|
|
+# in two-tailed test (alfa is now multiplied by 2).
|
|
|
|
+
|
|
|
|
+Y1=0.5*log((1+R1)/(1-R1)); Y2=0.5*log((1+R2)/(1-R2))
|
|
|
|
+
|
|
|
|
+B=(1/(N1-3)) + (1/(N2-3)); Z=abs(Y1-Y2)/sqrt(B) #abs() -> two-tailed test
|
|
|
|
+
|
|
|
|
+cl=cl; cl2=(1-cl)/2; Z_cl=qnorm(cl+cl2)
|
|
|
|
+
|
|
|
|
+if (is.na(Z)==FALSE){
|
|
|
|
+ if (Z > Z_cl){ value=1 } else { value=0 }
|
|
|
|
+} else {
|
|
|
|
+ value=NA
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#
|
|
|
|
+#p=pnorm(Z)-pnorm(-Z) #two-tailed / bilateral --> confidence level
|
|
|
|
+#p=1-p # --> significance level
|
|
|
|
+#p=p/2 # alfa/2 --> p_value
|
|
|
|
+#
|
|
|
|
+
|
|
|
|
+if (is.na(Z)==FALSE){
|
|
|
|
+ p=2*pnorm(Z,lower.tail=FALSE)
|
|
|
|
+} else {
|
|
|
|
+ p=NA
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+invisible(list(value=value, p_value=p))
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_testR<-function(R, N, cl = 0.95, test="two"){
|
|
|
|
+#
|
|
|
|
+# Testing the Null-Hypothesis of population correlation equal to zero
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Statistical test for sample correlation coefficient (R) according to
|
|
|
|
+# the Null Hypothesis of population correlation coefficient (RO) equal to zero.
|
|
|
|
+#
|
|
|
|
+# H0: RO==0 vs. H1: RO=/=0
|
|
|
|
+#
|
|
|
|
+# The statistic t follows a Student's (T-test) distribution:
|
|
|
|
+#
|
|
|
|
+# if t > t_alfa,N-2 (one-tailed) H0 is rejected
|
|
|
|
+#
|
|
|
|
+# if |t| > t_alfa/2,N-2 (two-tailed) H0 is rejected
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_testR(R, N, cl, test=c("one","two"))
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# R: number, sample correlation
|
|
|
|
+#
|
|
|
|
+# N: number, nº of elements in the sample
|
|
|
|
+#
|
|
|
|
+# cl: Value of the confidence level.
|
|
|
|
+# Default is cl=0.95 (95% confidence level; 5% significance level).
|
|
|
|
+#
|
|
|
|
+# test: string of characters that can be "one" or "two" (by default) according
|
|
|
|
+# to required test; one-tailed, unilateral (the former) o two-tailed,
|
|
|
|
+# bilateral (the latter).
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# value: binary output (0,1) / 0 = H0 accepted, 1 = H0 rejected
|
|
|
|
+#
|
|
|
|
+# p_value: area under the Student's T density curve to the right of endpoint |t|
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 08 December 2010
|
|
|
|
+#
|
|
|
|
+# Modified by jgarcia & lrodrigues@ic3.cat (13 July 2011) : Correcting p-values
|
|
|
|
+# in two-tailed test (alfa is now multiplied by 2).
|
|
|
|
+
|
|
|
|
+t=R*sqrt((N-2)/(1-R**2))
|
|
|
|
+
|
|
|
|
+cl=cl
|
|
|
|
+
|
|
|
|
+if (test=="one"){
|
|
|
|
+ t_cl=qt(cl, N-2)
|
|
|
|
+ if (abs(t) > t_cl){ value=1 } else { value=0 }
|
|
|
|
+#
|
|
|
|
+ if (t > 0){
|
|
|
|
+ p=pt(t, N-2, lower.tail=FALSE)
|
|
|
|
+ } else { p=pt(t, N-2, lower.tail=TRUE) }
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+if (test=="two"){
|
|
|
|
+ cl2=(1-cl)/2; t_cl=qt(cl+cl2, N-2)
|
|
|
|
+ if (abs(t) > t_cl){ value=1 } else { value=0 }
|
|
|
|
+#
|
|
|
|
+ if (t > 0){
|
|
|
|
+ p=pt(t, N-2, lower.tail=FALSE)*2
|
|
|
|
+ } else { p=pt(t, N-2, lower.tail=TRUE)*2 }
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#ouputs
|
|
|
|
+
|
|
|
|
+invisible(list(value=value, p_value=p))
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_AA<-function(data, lat, yes=NULL){
|
|
|
|
+#
|
|
|
|
+# Compute the area-average of a 2D (longitude x latitude) array
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Area-average (AA) of "data" over the spatial domain described by 'longitude' and
|
|
|
|
+# 'latitude'.
|
|
|
|
+# The weights for the field-mean (AA) are based on the cosine of the latitude at each
|
|
|
|
+# gridpoint, i.e. all data in a latitudinal circle (all longitudes) have the same
|
|
|
|
+# weight, cosine("lat").
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_AA(data, lat, yes)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# data: 2-dimensional array (longitude, latitude; in this order)
|
|
|
|
+#
|
|
|
|
+# lat: vector containing the latitude coordinates
|
|
|
|
+#
|
|
|
|
+# yes: vector containing the position of non-NaN values, so 'good' values
|
|
|
|
+# (in a ns=nlon*nlat, GrADS format)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# data_ave; value of the area-average
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Javier Garcia-Serrano <jgarcia@ic3.cat> 03 December 2010
|
|
|
|
+# Modified by Virginie Guemas <vguemas@ic3.cat> (March 2011) : Weight = 0 for NA grid point
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+ycos=cos((lat*pi)/180); nlat=dim(data)[2]; nlon=dim(data)[1]; ns=nlon*nlat
|
|
|
|
+weight=array(dim=c(nlon,nlat))
|
|
|
|
+for (i in 1:nlon){
|
|
|
|
+ weight[i,]=ycos
|
|
|
|
+}
|
|
|
|
+#weight=weight/sum(weight) #normalization
|
|
|
|
+
|
|
|
|
+data=data*weight #weightening
|
|
|
|
+
|
|
|
|
+miss.val=is.null(yes)
|
|
|
|
+if(miss.val=='TRUE'){
|
|
|
|
+ M=array(dim=c(ns))
|
|
|
|
+ weight2=array(dim=c(ns))
|
|
|
|
+ for (i in 1:nlat){
|
|
|
|
+ M[(nlon*(i-1)+1):(nlon*(i-1)+nlon)]=data[,i]
|
|
|
|
+ weight2[(nlon*(i-1)+1):(nlon*(i-1)+nlon)]=weight[,i]
|
|
|
|
+ }
|
|
|
|
+ not=is.na(M)
|
|
|
|
+ yes=array(dim=c(ns))
|
|
|
|
+ counter=0
|
|
|
|
+ for (i in 1:ns){
|
|
|
|
+ if (not[i]=='FALSE'){counter=counter+1; yes[counter]=i}
|
|
|
|
+ }
|
|
|
|
+ yes=yes[1:counter]
|
|
|
|
+ tmp=sum(M[yes])
|
|
|
|
+ sweight=sum(weight2[yes])
|
|
|
|
+ data_ave=tmp/sweight
|
|
|
|
+ #data_ave=sum(M[yes])
|
|
|
|
+} else {
|
|
|
|
+ yes=yes
|
|
|
|
+ M=array(dim=c(ns))
|
|
|
|
+ weight2=array(dim=c(ns))
|
|
|
|
+ for (i in 1:nlat){
|
|
|
|
+ M[(nlon*(i-1)+1):(nlon*(i-1)+nlon)]=data[,i]
|
|
|
|
+ weight2[(nlon*(i-1)+1):(nlon*(i-1)+nlon)]=weight[,i]
|
|
|
|
+ }
|
|
|
|
+ tmp=sum(M[yes])
|
|
|
|
+ sweight=sum(weight2[yes])
|
|
|
|
+ data_ave=tmp/sweight
|
|
|
|
+ #data_ave=sum(M[yes])
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+data_ave
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_neff<-function(vector, fig=FALSE){
|
|
|
|
+#
|
|
|
|
+# Compute the effective degrees of freedom of a time series
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Number of effective degrees of freedom of "vector".
|
|
|
|
+# Based on: Zieba, A. (2010): Effective number of observations and unbiased estimators
|
|
|
|
+# of variance for autocorrelated data - an overview. Metrol. Meas. Syst. Vol XVII, No. 1,
|
|
|
|
+# pp. 3-16; index 330930, ISSN 0860-8229. www.metrology.pg.gda.pl
|
|
|
|
+#
|
|
|
|
+# Usage: CFU_neff(vector)
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# vector: 1-dimensional array, time (nt dimension)
|
|
|
|
+#
|
|
|
|
+# fig: logical, if plot in the autocorrelation of "vector" is desired
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# neff; number of effective degrees of freedom
|
|
|
|
+#
|
|
|
|
+# Authors:
|
|
|
|
+#
|
|
|
|
+# Created by Caio A. S. Coelho <caio.coelho@cptec.inpe.br>
|
|
|
|
+# Implemented in CFU_Rfunc by Javier García-Serrano <jgarcia@ic3.cat> (August 2011)
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+dof=length(vector)
|
|
|
|
+
|
|
|
|
+a=acf(vector,lag.max=dof-1,plot=fig)$acf[2:dof,1,1]
|
|
|
|
+
|
|
|
|
+s=0
|
|
|
|
+for (k in 1:(dof-1)){
|
|
|
|
+ s=s+(((dof-k)/dof)*(a[k]**2))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+neff=(dof/(1+(2*s)))-1
|
|
|
|
+
|
|
|
|
+#if (neff>dof){neff=dof}
|
|
|
|
+
|
|
|
|
+#outputs
|
|
|
|
+
|
|
|
|
+#list(neff=neff)
|
|
|
|
+neff
|
|
|
|
+
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_BS <- function(obs,pred,thresholds = seq(0,1,0.1)) {
|
|
|
|
+#
|
|
|
|
+# Compute the BS and its decomposition given a set of probability predictions of binary event
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns the values of the BS and its standard decomposition as well as the addition of the
|
|
|
|
+# two winthin-bin extra components (Stephenson et al., 2008). It also solves the bias-corrected
|
|
|
|
+# decomposition of the BS (Ferro and Fricker, 2012). BSS having the climatology as the reference
|
|
|
|
+# forecast.
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_BS(obs,pred,thresholds = seq(0,1,0.1))
|
|
|
|
+#
|
|
|
|
+# Input:
|
|
|
|
+#
|
|
|
|
+# obs: Vector of binary observations (1 or 0)
|
|
|
|
+#
|
|
|
|
+# pred: Vector of probablistic predictions [0,1]
|
|
|
|
+#
|
|
|
|
+# thresholds: Values used to bin the forecasts. By default the bins are
|
|
|
|
+# {[0,0.1), [0.1, 0.2), ... [0.9, 1]}.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# $rel: standard reliability
|
|
|
|
+# $res: standard resolution
|
|
|
|
+# $unc: standard uncertainty
|
|
|
|
+# $bs: Brier score
|
|
|
|
+# $bs_check_res: rel-res+unc
|
|
|
|
+# $bss_res: res-rel/unc
|
|
|
|
+# $gres: generalized resolution
|
|
|
|
+# $bs_check_gres: rel-gres+unc
|
|
|
|
+# $bss_gres: gres-rel/unc
|
|
|
|
+# $rel_bias_corrected: bias-corrected rel
|
|
|
|
+# $gres_bias_corrected: bias-corrected gres
|
|
|
|
+# $unc_bias_corrected: bias-corrected unc
|
|
|
|
+# $bss_bias_corrected: gres_bias_corrected-rel_bias_corrected/unc_bias_corrected
|
|
|
|
+# $nk: number of forecast in each bin
|
|
|
|
+# $fkbar: average probability of each bin
|
|
|
|
+# $okbar: relative frequency that the observed event occurred
|
|
|
|
+# $bins: bins used
|
|
|
|
+# $pred: values with which the forecasts are verified
|
|
|
|
+# $obs: probability forecasts of the event
|
|
|
|
+#
|
|
|
|
+# References:
|
|
|
|
+#
|
|
|
|
+# Wilks (2006) Statistical Methods in the Atmospheric Sciences.
|
|
|
|
+# Stephenson et al. (2008). Two extra components in the Brier score decomposition. Weather
|
|
|
|
+# and Forecasting, 23: 752-757.
|
|
|
|
+# Ferro and Fricker (2012). A bias-corrected decomposition of the BS. Quarterly Journal of
|
|
|
|
+# the Royal Meteorological Society, DOI: 10.1002/qj.1924.
|
|
|
|
+#
|
|
|
|
+# Example:
|
|
|
|
+#
|
|
|
|
+# a=runif(10)
|
|
|
|
+# b=round(a)
|
|
|
|
+# x=CFU_BS(b,a)
|
|
|
|
+# x$bs-x$bs_check_res
|
|
|
|
+# x$bs-x$bs_check_gres
|
|
|
|
+# x$rel_bias_corrected-x$gres_bias_corrected+x$unc_bias_corrected
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Luis Rodrigues <lrodrigues@ic3.cat> 16 April 2012
|
|
|
|
+#
|
|
|
|
+ if (max(pred) > 1 | min(pred) < 0) {
|
|
|
|
+ cat("Predictions outside [0,1] range. \n Are you certain this is a probability forecast? \n")
|
|
|
|
+ } else if (max(obs) != 1 & min(obs) != 0) {
|
|
|
|
+ cat("Binary events must be either 0 or 1. \n Are you certain this is a binary event? \n")
|
|
|
|
+ } else {
|
|
|
|
+ nbins=length(thresholds)-1 # Number of bins
|
|
|
|
+ n=length(pred)
|
|
|
|
+ bins=as.list(paste("bin",1:nbins,sep=""))
|
|
|
|
+ for (i in 1:nbins) {
|
|
|
|
+ if (i == nbins) {
|
|
|
|
+ bins[[i]]=list(which(pred>=thresholds[i] & pred<=thresholds[i+1]))
|
|
|
|
+ } else {
|
|
|
|
+ bins[[i]]=list(which(pred>=thresholds[i] & pred<thresholds[i+1]))
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+#
|
|
|
|
+ fkbar=okbar=nk=array(0,dim=nbins)
|
|
|
|
+ for (i in 1:nbins) {
|
|
|
|
+ nk[i]=length(bins[[i]][[1]])
|
|
|
|
+ fkbar[i]=sum(pred[bins[[i]][[1]]])/nk[i]
|
|
|
|
+ okbar[i]=sum(obs[bins[[i]][[1]]])/nk[i]
|
|
|
|
+ }
|
|
|
|
+#
|
|
|
|
+ obar=sum(obs)/length(obs)
|
|
|
|
+ relsum=ressum=term1=term2=0
|
|
|
|
+ for (i in 1:nbins) {
|
|
|
|
+ if (nk[i] > 0) {
|
|
|
|
+ relsum=relsum+nk[i]*(fkbar[i]-okbar[i])^2
|
|
|
|
+ ressum=ressum+nk[i]*(okbar[i]-obar)^2
|
|
|
|
+ for (j in 1:nk[i]) {
|
|
|
|
+ term1=term1+(pred[bins[[i]][[1]][j]]-fkbar[i])^2
|
|
|
|
+ term2=term2+(pred[bins[[i]][[1]][j]]-fkbar[i])*(obs[bins[[i]][[1]][j]]-okbar[i])
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ rel=relsum/n
|
|
|
|
+ res=ressum/n
|
|
|
|
+ unc=obar*(1-obar)
|
|
|
|
+ bs=sum((pred-obs)^2)/n
|
|
|
|
+ bs_check_res=rel-res+unc
|
|
|
|
+ bss_res=(res-rel)/unc
|
|
|
|
+ gres=res-term1*(1/n)+term2*(2/n) # Generalized resolution
|
|
|
|
+ bs_check_gres=rel-gres+unc # BS using GRES
|
|
|
|
+ bss_gres=(gres-rel)/unc # BSS using GRES
|
|
|
|
+#
|
|
|
|
+# Estimating the bias-corrected components of the BS
|
|
|
|
+#
|
|
|
|
+ term3=array(0,nbins)
|
|
|
|
+ for (i in 1:nbins) {
|
|
|
|
+ term3[i]=(nk[i]/(nk[i]-1))*okbar[i]*(1-okbar[i])
|
|
|
|
+ }
|
|
|
|
+ term_a=sum(term3,na.rm=T)/n
|
|
|
|
+ term_b=(obar*(1-obar))/(n-1)
|
|
|
|
+ rel_bias_corrected=rel-term_a
|
|
|
|
+ gres_bias_corrected=gres-term_a+term_b
|
|
|
|
+ if (rel_bias_corrected < 0 || gres_bias_corrected < 0 ) {
|
|
|
|
+ rel_bias_corrected2=max(rel_bias_corrected,rel_bias_corrected-gres_bias_corrected,0)
|
|
|
|
+ gres_bias_corrected2=max(gres_bias_corrected,gres_bias_corrected-rel_bias_corrected,0)
|
|
|
|
+ rel_bias_corrected=rel_bias_corrected2
|
|
|
|
+ gres_bias_corrected=gres_bias_corrected2
|
|
|
|
+ }
|
|
|
|
+ unc_bias_corrected=unc+term_b
|
|
|
|
+ bss_bias_corrected=(gres_bias_corrected-rel_bias_corrected)/unc_bias_corrected
|
|
|
|
+#
|
|
|
|
+ if (round(bs,8)==round(bs_check_gres,8) & round(bs_check_gres,8)==round((rel_bias_corrected-gres_bias_corrected+unc_bias_corrected),8)) {
|
|
|
|
+# cat("No error found \n")
|
|
|
|
+# cat("BS = REL-GRES+UNC = REL_lessbias-GRES_lessbias+UNC_lessbias \n")
|
|
|
|
+ }
|
|
|
|
+#
|
|
|
|
+ invisible(list(rel=rel,res=res,unc=unc,bs=bs,bs_check_res=bs_check_res,bss_res=bss_res,gres=gres,bs_check_gres=bs_check_gres,bss_gres=bss_gres,rel_bias_corrected=rel_bias_corrected,gres_bias_corrected=gres_bias_corrected,unc_bias_corrected=unc_bias_corrected,bss_bias_corrected=bss_bias_corrected,nk=nk,fkbar=fkbar,okbar=okbar,bins=bins,pred=pred,obs=obs))
|
|
|
|
+#
|
|
|
|
+ } # end of if
|
|
|
|
+#
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_plotmapmostliktercile <- function(lon,lat,ptc,
|
|
|
|
+ tit='Probability forecast of most likely tercile (%)',file='G_PFC_map.ps') {
|
|
|
|
+#
|
|
|
|
+# Produce probabilistic forecast map of most likely tercile of a forecast.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Generates vector graphics file of a probabilistic forecast map showing the most likely tercile.
|
|
|
|
+# Based on an adapted version of Caio Coelho's EURO-BRISA function.
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_plotmapmostliktercile(lon,lat,ptc,tit='Probability forecast of most likely tercile (%)',file='G_PFC_map.ps')
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ptc: Array of forecast terciles information for each ensemble member, output of CFU_PB().
|
|
|
|
+#
|
|
|
|
+# tit: Title of the Figure.
|
|
|
|
+#
|
|
|
|
+# file: Name of the postscript output file.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# Vector graphics file of a probabilistic forecast map showing the most likely tercile.
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 29 Oct 2012.
|
|
|
|
+#
|
|
|
|
+ source('/cfu/pub/scripts/R/rclim.txt')
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nfc <- dim(ptc)[1]
|
|
|
|
+ nlonr <- nlon-ceiling(length(lon[lon<180 & lon > 0]))
|
|
|
|
+ # Probability of most likely tercile.
|
|
|
|
+ PMLT <- (ptc[1,1,,]*NA)
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ lo <- sum(ptc[,1,j,i],na.rm=T)/nfc
|
|
|
|
+ mid <- sum(ptc[,2,j,i],na.rm=T)/nfc
|
|
|
|
+ up <- sum(ptc[,3,j,i],na.rm=T)/nfc
|
|
|
|
+ if (all(up > c(mid,lo))) {
|
|
|
|
+ PMLT[j,i] <- up
|
|
|
|
+ } else if (all(mid > c(up,lo))) {
|
|
|
|
+ PMLT[j,i] <- 0
|
|
|
|
+ } else if (all(lo > c(mid,up))) {
|
|
|
|
+ PMLT[j,i] <- -lo
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ # EURO-BRISA figure. Modified Caio Coelho's function.
|
|
|
|
+ plotmapmostliktercile <- function(lon,lat,data,maintit='',legtit='',equi=TRUE,bw=FALSE,
|
|
|
|
+ cont=FALSE,reg=FALSE,...,lonlim=c(min(floor(lon)),max(ceiling(lon))),latlim=c(min(floor(lat)),max(ceiling(lat))),orientation=NULL,
|
|
|
|
+ mapdat="world",xmaplim=c(min(floor(lon)),max(ceiling(lon))),ymaplim=c(min(floor(lat)),max(ceiling(lat))),
|
|
|
|
+ longrds=NULL,latgrds=NULL,
|
|
|
|
+ breaks=NULL,n=11,colours=NULL,
|
|
|
|
+ roundleg=1,invert=FALSE) {
|
|
|
|
+ rg<-range(data, na.rm=T)
|
|
|
|
+ lowerlim <- rg[1]; upperlim <- rg[2]
|
|
|
|
+ maximum <- max(abs(c(lowerlim,upperlim)))
|
|
|
|
+ if(is.null(breaks)){
|
|
|
|
+ if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n)
|
|
|
|
+ else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n)
|
|
|
|
+ }
|
|
|
|
+ if(is.null(colours)){
|
|
|
|
+ if(!bw){
|
|
|
|
+ # Check if range includes negative values to use appropriate colour scale.
|
|
|
|
+ if (rg[1] <0) {
|
|
|
|
+ if (!invert) {
|
|
|
|
+ colours <- bluered(seq(0,n-1,by=1), "linear", yellow =TRUE)
|
|
|
|
+ } else {
|
|
|
|
+ colours <- rev(bluered(seq(0,n-1,by=1), "linear", yellow =TRUE))
|
|
|
|
+ }
|
|
|
|
+ } else {
|
|
|
|
+ colours <- bluered(seq(0,n-1,by=1), "linear", white=0, invert=T)
|
|
|
|
+ }
|
|
|
|
+ } else {
|
|
|
|
+ if (rg[1] <0) {
|
|
|
|
+ colours <- grey(seq(0, 1, length = length(breaks)-1))
|
|
|
|
+ colours <- c(colours[1:((n-1)/2)],rev(colours[(((n-1)/2)+1):n]))
|
|
|
|
+ } else {
|
|
|
|
+ colours <- grey(seq(0, 1, length = length(breaks)-1))
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ if (equi) {
|
|
|
|
+ layout(matrix(1:2,ncol=1,nrow=2),heights=c(5,1))
|
|
|
|
+ par(mar=c(6,2,6.5,1),cex.axis=1.0)
|
|
|
|
+ if (reg) {
|
|
|
|
+ layout(matrix(1:2,ncol=1,nrow=2),heights=c(5,1))
|
|
|
|
+ par(mar=c(2,3,4,1),las=1,cex=1.4)
|
|
|
|
+ }
|
|
|
|
+ # Create probabilistic map.
|
|
|
|
+ image(lon,lat,data,axes=F,col=colours,xlab='',ylab='',breaks=breaks)
|
|
|
|
+ par(cex=2.1)
|
|
|
|
+ if (cont) {contour(lon,lat,data,levels=round(breaks,1),labels=round(breaks,1),add=T,...)}
|
|
|
|
+ # Check if lon is from 0 to 360 or -180 to 180 to use appropriate world map.
|
|
|
|
+ if (min(lon)<0) {
|
|
|
|
+ map('world',interior = F,add = T, lwd=2) # Low resolution world map (lon -180 to 180).
|
|
|
|
+ } else {
|
|
|
|
+ map('world2',interior = F,add = T, lwd=2) # Low resolution world map (lon 0 to 360).
|
|
|
|
+ }
|
|
|
|
+ box()
|
|
|
|
+ map.axes()
|
|
|
|
+ par(font.main=1,cex=1.4)
|
|
|
|
+ title(maintit)
|
|
|
|
+ # Adding colorbar.
|
|
|
|
+ par(mar=c(4.5,2,0,1),mgp=c(1.5,0.3,0),las=1)
|
|
|
|
+ if (reg) {
|
|
|
|
+ par(mar=c(3,0.3,1.5,0.3),mgp=c(1.5,0.3,0),las=1)
|
|
|
|
+ }
|
|
|
|
+ image(c(1:n),1,t(t(c(1:n))),axes=F,col=colours,xlab='',ylab='')
|
|
|
|
+ box()
|
|
|
|
+ par(font.main=1,cex=1.3)
|
|
|
|
+ if (maximum > 1) {
|
|
|
|
+ axis(1,at=seq(1,9),tick=F,pos=0.8,cex.axis=1.5,
|
|
|
|
+ labels = c("85-100","70-85","55-70","40-55"," ","40-55","55-70","70-85","85-100"))
|
|
|
|
+ axis(3,at=5,tick=F,cex.axis=1.5,pos=-1.3,labels=c("White=central tercile most likely"))
|
|
|
|
+ axis(3,at=2.5,tick=F,cex.axis=1.5,pos=1.3,labels=c("Below normal"))
|
|
|
|
+ axis(3,at=7.5,tick=F,cex.axis=1.5,pos=1.3,labels=c("Above normal"))
|
|
|
|
+ } else {
|
|
|
|
+ axis(1,at=seq(0.5, length(breaks)-0.5),tick=F,
|
|
|
|
+ labels=round(breaks[1:(length(breaks))],roundleg))
|
|
|
|
+ }
|
|
|
|
+ # Redefine font size.
|
|
|
|
+ par(cex=1.4)
|
|
|
|
+ # Add the title.
|
|
|
|
+ title(xlab = legtit,cex.lab=1.4)
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ postscript(paste(file,sep=''),paper="special",width=12,height=9.25,horizontal=F)
|
|
|
|
+ plotmapmostliktercile(c(lon[c((nlonr+1):nlon)]-360,lon[1:nlonr]),rev(lat),t(PMLT[nlat:1,c((nlonr+1):nlon,1:nlonr)]*100),
|
|
|
|
+ reg=T,maintit=tit,invert=F,breaks=c(-100,-85,-70,-55,-40,40,55,70,85,100),n=9)
|
|
|
|
+ dev.off()
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_PB <- function(lon,lat,ano,fcavt=2:4,fcyr,fcsys=1,qt=3,ncpu=4) {
|
|
|
|
+#
|
|
|
|
+# Compute probabilistic bins (qt) of a forecast (fcyr) relative to the forecast climatology of the residual
|
|
|
|
+# forecasts excluding the selected forecast year (fcyr).
|
|
|
|
+# The maximal values of the lowest/highest quantiles is +/-10^+06 due to limitation in the function hist().
|
|
|
|
+# By default (qt=3) CFU_PB() computes terciles that can be plotted with CFU_plotmapmostliktercile().
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns probabilistic bins of a forecast (fcyr) relative to the forecast climatology without selected forecast.
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_PB(lon,lat,ano,fcavt=2:4,fcyr,fcsys=1,qt=3,ncpu=4)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ano: Array of forecast anomalies in CFU_load format.
|
|
|
|
+#
|
|
|
|
+# fcavt: A vector of time steps to average across.
|
|
|
|
+#
|
|
|
|
+# fcyr: Number of the forecast start year of the forecast.
|
|
|
|
+#
|
|
|
|
+# fcsys: Number of the forecast system.
|
|
|
|
+#
|
|
|
|
+# qt: Number of probabilistic bins 3: terciles, 4: quartiles, etc.
|
|
|
|
+#
|
|
|
|
+# ncpu: Number of CPUs to be used.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# Probabilistic forecast information (1 or 0) at each grid point with
|
|
|
|
+# dimensions (number of ensemble members, qt, number of longitudes, latitudes).
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 29 Oct 2012.
|
|
|
|
+#
|
|
|
|
+ library(doMC)
|
|
|
|
+ registerDoMC(ncpu)
|
|
|
|
+ # Probabilistic forecast.
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nyr <- dim(ano)[3]; nfc <- dim(ano)[2]; lt <- dim(ano)[4]
|
|
|
|
+ PTC <- array(NA, dim=c(nfc,3,nlat,nlon))
|
|
|
|
+ for (im in 1:nfc) {
|
|
|
|
+ XY.p <- foreach(i=1:nlon) %dopar% {
|
|
|
|
+ PL <- array(NA, dim=c(qt,nlat))
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ # Mean over forecast range.
|
|
|
|
+ f.p <- apply(ano[fcsys,im,-fcyr,fcavt,j,i],c(1),mean)
|
|
|
|
+ # PDF of model climatology without selected forecast.
|
|
|
|
+ qum <- quantile(f.p,probs=c((1:(qt-1))/qt),prob,na.rm=F,names=F,type=8)
|
|
|
|
+ # Bins of selected forecast.
|
|
|
|
+ f.p <- mean(ano[fcsys,im,fcyr,fcavt,j,i])
|
|
|
|
+ if (all(c(qum,f.p)==0)) {
|
|
|
|
+ f0 <- rep(0,qt)
|
|
|
|
+ if (!is.integer(qt/2)) { f0[median(1:qt)] <- 1 }
|
|
|
|
+ PL[,j] <- f0
|
|
|
|
+ } else {
|
|
|
|
+ PL[,j] <- hist(f.p,breaks=c(-1e+06,qum,1e+06),plot=F)$counts
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ PL
|
|
|
|
+ }
|
|
|
|
+ for (i in 1:nlon) { PTC[im,,,i] <- XY.p[[i]] }
|
|
|
|
+ }
|
|
|
|
+ PTC
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_MODE <- function(ano,eof,mode=1) {
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Project anomalies onto modes to get temporal evolution of the EOF mode selected.
|
|
|
|
+# Returns principal components (PCs) by area-weighted projection onto EOF pattern (from CFU_EOF()).
|
|
|
|
+# Able to handle NAs.
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# ano: Array of observed or forecast anomalies from CFU_ano() with dimensions
|
|
|
|
+# (number of forecast systems, ensemble members, start years, forecast months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# eof: EOF object from CFU_EOF()
|
|
|
|
+#
|
|
|
|
+# mode: mode number in EOF object onto which to project
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# PC.ver: Array of PCs in verification format
|
|
|
|
+# (number of forecast systems, ensemble members, start years, forecast months).
|
|
|
|
+# PC.ver.em: Array of ensemble-mean PCs in verification format
|
|
|
|
+# (number of forecast systems, start years, forecast months).
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 19 Nov 2012.
|
|
|
|
+#
|
|
|
|
+# Lauriane Batte <lauriane.batte@ic3.cat> March 2014. Bug-fixes: 1-Extra weighting of the anomalies before projection.
|
|
|
|
+# 2-Reversion of the anomalies along latitudes.
|
|
|
|
+# 3-Extra-normalisation not necessary.
|
|
|
|
+#
|
|
|
|
+# Virginie Guemas <virginie.guemas@ic3.cat> 17 March 2014. Bug-fixes:1-Another extra-normalisation.
|
|
|
|
+# 2-15 lines to compute the em reduced to 1.
|
|
|
|
+# Lauriane Batte <lauriane.batte@ic3.cat> 18 March 2014. Normalization by std before returning PCs to be coherent with CFU_EOF
|
|
|
|
+#
|
|
|
|
+# Virginie Guemas <virginie.guemas@ic3.cat> 11 April 2014. 1 - Removal of lon, lat, ncpu and neofs argument unused
|
|
|
|
+# 2 - Security checks ano and eof consistency
|
|
|
|
+# 3 - Removal of the mask which is already contained in the EOFs
|
|
|
|
+# 4 - Removal of the PC normalization since we have chosen in CFU_EOF
|
|
|
|
+# to normalize the EOFs and multiply the PCs by the normalization
|
|
|
|
+# factor and the eigenvalue so that the restitution of the original
|
|
|
|
+# field is done simply by PC * EOFs
|
|
|
|
+# 5 - The new convention in CFU_EOF is to divide by the weights so that the
|
|
|
|
+# reconstruction of the original field rather than the weighted field
|
|
|
|
+# is obtained by PC * EOFs. The EOFs need therefore to be multiplied
|
|
|
|
+# back by the weights before projection so that EOF * t(EOF) = 1
|
|
|
|
+# 6 - Since W *X = PC * EOF if EOF is multiplied back by the weights,
|
|
|
|
+# PC = W * X * t(EOF) and X the input field to be projected (X)
|
|
|
|
+# needs to be multiplied by W.
|
|
|
|
+#
|
|
|
|
+# library(s2dverification) # To simplify some lines. Virginie
|
|
|
|
+# Getting input dimensions. Fabian
|
|
|
|
+ nlon<- dim(ano)[6]; nlat <- dim(ano)[5]
|
|
|
|
+ nyr <- dim(ano)[3]; lt <- dim(ano)[4]
|
|
|
|
+ nfc <- dim(ano)[2]; nmod <- dim(ano)[1]
|
|
|
|
+ # Security checks. Virginie
|
|
|
|
+ if ( dim(eof$EOFs)[2] != nlat ) {
|
|
|
|
+ stop("Inconsistent number of latitudes between eof and input field")
|
|
|
|
+ }
|
|
|
|
+ if ( dim(eof$EOFs)[3] != nlon ) {
|
|
|
|
+ stop("Inconsistent number of longitudes between eof and input field")
|
|
|
|
+ }
|
|
|
|
+ # Initialization of pc.ver. Fabian.
|
|
|
|
+ pc.ver <- array(NA,dim=c(nmod,nfc,nyr,lt))
|
|
|
|
+ # Weights are already in EOF and shouldn't be accounted for twice. Lauriane
|
|
|
|
+ #e.1 <- eof$EOFs[mode,,]; mask <- eof$mask; wght <- eof$wght
|
|
|
|
+ # Now the EOF outputs from CFU_EOF are masked . Virginie
|
|
|
|
+ #e.1 <- eof$EOFs[mode,,]; mask <- eof$mask; wght <- rep(1,nlon*nlat)
|
|
|
|
+ #dim(wght) <- c(nlon,nlat)
|
|
|
|
+ # Since we have chosen as a new convention in CFU_EOF to divide by the
|
|
|
|
+ # weights so that the restitution of the original field is only obtained
|
|
|
|
+ # by PCs*EOFs, I multiply here back by the weigths to ensure that the
|
|
|
|
+ # projection is correct. Virginie
|
|
|
|
+ e.1 <- eof$EOFs[mode,,]*eof$wght
|
|
|
|
+ # Since W *X = PC * EOF after such multiplication by the weights,
|
|
|
|
+ # the PC is obtained by W * X * t(EOF) so we multiply the field to project
|
|
|
|
+ # by the weights. Virginie
|
|
|
|
+ ano <- ano*InsertDim(InsertDim(InsertDim(InsertDim(eof$wght,1,nmod),2,nfc),3,nyr),4,lt)
|
|
|
|
+ # The 100 lines below can be simplified in a few lines as done below.
|
|
|
|
+ # Indeed, we can loop on nmod, nfc, nyr and lt even if the length of these dimentions is only 1.
|
|
|
|
+ # Furthermore, the PC of the ensemble mean anomaly equals the ensemble means of the PC of single
|
|
|
|
+ # members. [ 1/n sum(Xi) ] * t(EOF) = 1/n sum [ Xi * t(EOF) ]. Virginie
|
|
|
|
+ # # For one realization only. Fabian
|
|
|
|
+ # if (nfc==1) {
|
|
|
|
+ # # Case for only one member should not be separate from several members. We can loop over 1 member only. Virginie.
|
|
|
|
+ # pc.ver.em <- FALSE
|
|
|
|
+ # # Why is there no loop over the models ? Virginie.
|
|
|
|
+ # for (iy in 1:nyr) {
|
|
|
|
+ # for (il in 1:lt) {
|
|
|
|
+ # if (!all(is.na(ano[1,1,iy,il,,]))) {
|
|
|
|
+ # #t1 <- (t(ano[nmod,nfc,iy,il,nlat:1,])*e.1)
|
|
|
|
+ # # Latitudes should not be reversed to stay coherent with CFU_EOF. Lauriane
|
|
|
|
+ # t1 <- (ano[nmod,nfc,iy,il,,]*e.1)
|
|
|
|
+ # # The division by sum(mask) is an additional extra-normalisation. Virginie
|
|
|
|
+ # #nominator <- (sum(t1*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # # The EOFs are already masked now in CFU_EOF. Virginie
|
|
|
|
+ # #nominator <- sum(t1*mask,na.rm=TRUE)
|
|
|
|
+ # nominator <- sum(t1,na.rm=TRUE)
|
|
|
|
+ # # Normalization should not be necessary. Projection is a simple scalar product. Lauriane
|
|
|
|
+ # #denominator <- (sum(e.1^2*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # denominator <- 1.
|
|
|
|
+ # pc.ver[1,1,iy,il] <- nominator/denominator
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # # Normalization by standard deviation. Lauriane
|
|
|
|
+ # #for (il in 1:lt) {
|
|
|
|
+ # # pc.ver[1,1,,il] <- pc.ver[1,1,,il]/sd(pc.ver[1,1,,il])
|
|
|
|
+ # #}
|
|
|
|
+ # # I have removed this normalization in agreement with the conventions chosen for CFU_EOF.
|
|
|
|
+ # # The EOFs are normalized to 1 and the simple product EOF*PC should restitute the original field.
|
|
|
|
+ # # Virginie
|
|
|
|
+ # } else {
|
|
|
|
+ # pc.ver.em <- array(NA,dim=c(nmod,nyr,lt))
|
|
|
|
+ # # Virginie. This single line below gives the same result as l. 1162-1176
|
|
|
|
+ # ano.em <- apply(ano,c(1,3,4,5,6),mean,na.rm=TRUE)
|
|
|
|
+ # for (imo in 1:nmod) {
|
|
|
|
+ # # Ensemble mean anomalies.
|
|
|
|
+ # #if (lt==1 | nyr==1 ) {
|
|
|
|
+ # # if (lt==1 & !nyr==1 ) {
|
|
|
|
+ # # ano.em <- apply(ano[imo,,,,,],c(2,3,4),mean,na.rm=TRUE)
|
|
|
|
+ # # dim(ano.em) <- c(nyr,1,nlat,nlon)
|
|
|
|
+ # # } else if (!lt==1 & nyr==1 ) {
|
|
|
|
+ # # ano.em <- apply(ano[imo,,,,,],c(2,3,4),mean,na.rm=TRUE)
|
|
|
|
+ # # dim(ano.em) <- c(1,lt,nlat,nlon)
|
|
|
|
+ # # } else if (lt==1 & nyr==1 ) {
|
|
|
|
+ # # ano.em <- apply(ano[imo,,,,,],c(2,3),mean,na.rm=TRUE)
|
|
|
|
+ # # dim(ano.em) <- c(1,1,nlat,nlon)
|
|
|
|
+ # # }
|
|
|
|
+ # #} else {
|
|
|
|
+ # # ano.em <- apply(ano[imo,,,,,],c(2,3,4,5),mean,na.rm=TRUE)
|
|
|
|
+ # #}
|
|
|
|
+ # for (iy in 1:nyr) {
|
|
|
|
+ # for (il in 1:lt) {
|
|
|
|
+ # if (!all(is.na(ano.em[imo,iy,il,,]))) {
|
|
|
|
+ # #t1 <- (t(ano.em[imo,iy,il,nlat:1,])*e.1)
|
|
|
|
+ # # Latitudes should not be reversed to stay coherent with CFU_EOF. Lauriane.
|
|
|
|
+ # t1 <- (ano.em[imo,iy,il,,]*e.1)
|
|
|
|
+ # # The division by sum(mask) is an additional extra-normalisation. Virginie
|
|
|
|
+ # #nominator <- (sum(t1*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # # The EOFs are already masked now in CFU_EOF. Virginie
|
|
|
|
+ # #nominator <- sum(t1*mask,na.rm=TRUE)
|
|
|
|
+ # nominator <- sum(t1,na.rm=TRUE)
|
|
|
|
+ # # Normalization should not be necessary. Projection is a simple scalar product. Lauriane
|
|
|
|
+ # #denominator <- (sum(e.1^2*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # denominator <- 1.
|
|
|
|
+ # pc.ver.em[imo,iy,il] <- nominator/denominator
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # # Normalization by standard deviation: done for each model independently. Lauriane
|
|
|
|
+ # #for (il in 1:lt) {
|
|
|
|
+ # # pc.ver.em[imo,,il] <- pc.ver.em[imo,,il]/sd(pc.ver.em[imo,,il])
|
|
|
|
+ # #}
|
|
|
|
+ # # I have removed this normalization in agreement with the conventions chosen for CFU_EOF.
|
|
|
|
+ # # The EOFs are normalized to 1 and the simple product EOF*PC should restitute the original field.
|
|
|
|
+ # # Virginie
|
|
|
|
+ # # Ensemble member anomalies.
|
|
|
|
+ # for (im in 1:nfc) {
|
|
|
|
+ # for (iy in 1:nyr) {
|
|
|
|
+ # for (il in 1:lt) {
|
|
|
|
+ # if (!all(is.na(ano[imo,im,iy,il,,]))) {
|
|
|
|
+ # #t1 <- (t(ano[imo,im,iy,il,nlat:1,])*e.1)
|
|
|
|
+ # # Latitudes should not be reversed to stay coherent with CFU_EOF. Lauriane.
|
|
|
|
+ # t1 <- (ano[imo,im,iy,il,,]*e.1)
|
|
|
|
+ # # The division by sum(mask) is an additional extra-normalisation. Virginie
|
|
|
|
+ # #nominator <- (sum(t1*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # # The EOFs are already masked now in CFU_EOF. Virginie
|
|
|
|
+ # #nominator <- sum(t1*mask,na.rm=TRUE)
|
|
|
|
+ # nominator <- sum(t1,na.rm=TRUE)
|
|
|
|
+ # # Normalization should not be necessary. Projection is a simple scalar product. Lauriane
|
|
|
|
+ # #denominator <- (sum(e.1^2*wght*mask,na.rm=TRUE)/sum(wght*mask,na.rm=TRUE))
|
|
|
|
+ # denominator <- 1.
|
|
|
|
+ # pc.ver[imo,im,iy,il] <- nominator/denominator
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # # Normalization by standard deviation: done for each model and ensemble member independently. Lauriane
|
|
|
|
+ # #for (il in 1:lt) {
|
|
|
|
+ # # pc.ver[imo,im,,il] <- pc.ver[imo,im,,il]/sd(pc.ver[imo,im,,il])
|
|
|
|
+ # #}
|
|
|
|
+ # # I have removed this normalization in agreement with the conventions chosen for CFU_EOF.
|
|
|
|
+ # # The EOFs are normalized to 1 and the simple product EOF*PC should restitute the original field.
|
|
|
|
+ # # Virginie
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ # }
|
|
|
|
+ for (imo in 1:nmod) {
|
|
|
|
+ for (im in 1:nfc) {
|
|
|
|
+ for (iy in 1:nyr) {
|
|
|
|
+ for (il in 1:lt) {
|
|
|
|
+ if (!all(is.na(ano[imo,im,iy,il,,]))) {
|
|
|
|
+ pc.ver[imo,im,iy,il] <- sum(ano[imo,im,iy,il,,]*e.1,na.rm=TRUE)
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ # For compatibility with previous version, I keep the output of the ensemble mean although it is useless
|
|
|
|
+ # because the operation can be done easily outside of CFU_MODE if necessary. This should be removed in
|
|
|
|
+ # future. Virginie
|
|
|
|
+ pc.ver.em=Mean1Dim(pc.ver,2)
|
|
|
|
+
|
|
|
|
+ return(list(PC.ver=pc.ver,PC.ver.em=pc.ver.em))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_DTF <- function(lon,lat,ano,gmt,ncpu=4) {
|
|
|
|
+#
|
|
|
|
+# Linearly regress out forecast global mean temperature trend (GMT) from forecast monthly anomalies
|
|
|
|
+# of ten-year hindcasts started each year. The regression depends on the forecast lead month across all forecasts.
|
|
|
|
+# NAs are excluded from the linear model fit, but kept in the output.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns linearly detrended forecast anomalies in the verification format (from CFU_ano()).
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_DTF(lon,lat,ano,gmt)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ano: Array of forecast anomalies from CFU_ano() with dimensions
|
|
|
|
+# (number of forecast systems, ensemble members, start years, forecast months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# gmt: Array of forecast global mean temperature anomalies with dimensions
|
|
|
|
+# (number of forecast systems, ensemble members, start years, forecast months).
|
|
|
|
+#
|
|
|
|
+# ncpu: Number of CPUs to be used.
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# ano.dt.ver: Array of detrended forecast anomalies from CFU_ano() with dimensions
|
|
|
|
+# (number of forecast systems, ensemble members, start years, forecast months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 16 Nov 2012.
|
|
|
|
+#
|
|
|
|
+ library(doMC)
|
|
|
|
+ registerDoMC(ncpu)
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nyr <- dim(ano)[3]; lt <- dim(ano)[4]
|
|
|
|
+ nfc <- dim(ano)[2]; nmod <- dim(ano)[1]
|
|
|
|
+ ano.dt.ver <- array(NA,dim=c(nmod,nfc,nyr,lt,nlat,nlon))
|
|
|
|
+ ## Map of regression slopes.
|
|
|
|
+ r.XY <- array(NA,dim=c(nmod,lt,nlat,nlon))
|
|
|
|
+ for (imo in 1:nmod) {
|
|
|
|
+ XY.p <- foreach(i=1:nlon, .verbose=F) %dopar% {
|
|
|
|
+ PL1 <- array(NA,dim=c(nfc,nyr,lt,nlat))
|
|
|
|
+ PL2 <- array(NA,dim=c(lt,nlat))
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ for (iy in 1:lt) {
|
|
|
|
+ # Forecast anomaly.
|
|
|
|
+ fctpr <- as.vector(ano[imo,,,iy,j,i])
|
|
|
|
+ # Forecast trend wrt. 1961-2011.
|
|
|
|
+ trend.f <- as.vector(gmt[imo,,,iy])
|
|
|
|
+ trend <- trend.f-mean(trend.f)
|
|
|
|
+ # Regress forecast anomaly onto trend and remove it from the forecast. Keep residual forecast anomaly.
|
|
|
|
+ if (!all(is.na(fctpr))) {
|
|
|
|
+ fit <- lm(fctpr~trend,na.action=na.exclude)
|
|
|
|
+ fctpr.dt <- residuals(fit)
|
|
|
|
+ dim(fctpr.dt) <- c(nfc,nyr)
|
|
|
|
+ PL1[,,iy,j] <- fctpr.dt
|
|
|
|
+ # Regression coefficient map.
|
|
|
|
+ PL2[iy,j] <- coefficients(fit)[[2]]
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ return(list(DFC=PL1,REG=PL2))
|
|
|
|
+ }
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ ano.dt.ver[imo,,,,,i] <- XY.p[[i]]$DFC
|
|
|
|
+ r.XY[imo,,,i] <- XY.p[[i]]$REG
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ gc()
|
|
|
|
+ return(list(ano.dt.ver=ano.dt.ver,reg=r.XY))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_DTO <- function(lon,lat,ano,yr1,yr2,yrint,gmt=F) {
|
|
|
|
+#
|
|
|
|
+# Linearly regress out global mean temperature trend (GMT) from observed monthly anomalies
|
|
|
|
+# that coincide with hindcasts started on November 1st every 1,2,3, etc. years.
|
|
|
|
+# The anomalies are selected from the first year that coincide with each hindcast.
|
|
|
|
+# NAs are excluded from the linear model fit, but kept in the output.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns linearly detrended and original anomalies in a format to be used for CFU_EOF().
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_DTO(lon,lat,ano,yr1,yr2,yrint,gmt=F)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ano: Array of anomalies from CFU_ano() with dimensions
|
|
|
|
+# (1, 1, number of start years, forecast months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# yr1: First start year.
|
|
|
|
+#
|
|
|
|
+# yr2: Last start year.
|
|
|
|
+#
|
|
|
|
+# yrint: Interval between hindcasts in years.
|
|
|
|
+#
|
|
|
|
+# gmt: Global mean temperature time series. (optional)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# ano.dt.ver: An array of the linearly detrended anomalies in CFU verification format.
|
|
|
|
+#
|
|
|
|
+# ano.dt: An array of the linearly detrended anomalies with dimensions
|
|
|
|
+# (number of months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# ano: An array of the original input anomalies with dimensions
|
|
|
|
+# (number of months, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# reg: An array containing a map of regression coefficients
|
|
|
|
+# (number of latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 19 Apr 2013.
|
|
|
|
+#
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nyr <- dim(ano)[3]; lt <- dim(ano)[4]
|
|
|
|
+ print(nyr); print(lt)
|
|
|
|
+ if (all(gmt==F)) {
|
|
|
|
+ # Global mean temperature GMT from HadCRUT, climate explorer.
|
|
|
|
+ # http://climexp.knmi.nl/data/icrutem3_hadsst2_0-360E_-90-90N_n.dat
|
|
|
|
+ T1 <- read.table("/cfu/data/ukmo/hadcrut/icrutem3_hadsst2_0-360E_-90-90N_n.dat",row.names=1)
|
|
|
|
+ XT <- array(NA,dim=c(12,113))
|
|
|
|
+ for (im in 1:12) {
|
|
|
|
+ XT[im,] <- T1[[im]]
|
|
|
|
+ }
|
|
|
|
+ XT[XT==-9.999000e+02] <- NA
|
|
|
|
+ # GMT anomaly wrt 1356 mo, 113 yr.
|
|
|
|
+ tgl.obs <- ts(c(as.vector(XT),rep(NA,120)),deltat=1/12,start=c(1900,1),end=c(2022,12))
|
|
|
|
+ # GMT anomaly wrt verification period.
|
|
|
|
+ tgl.vp.nov.obs <- window(tgl.obs,start=c(yr1,11),end=c((yr2+lt/12),10))
|
|
|
|
+ tgl.vp.nov.obs <- tgl.vp.nov.obs*((yr2-yr1+lt/12)/113)
|
|
|
|
+ gmt <- tgl.vp.nov.obs
|
|
|
|
+ }
|
|
|
|
+ o.XY <- array(NA,dim=c(((yr2-yr1+lt/12)*12),nlat,nlon))
|
|
|
|
+ TP1 <- ano[1,1,,1:(yrint*12),nlat:1,]
|
|
|
|
+ TP1 <- aperm(TP1, c(2,1,3,4))
|
|
|
|
+ dim(TP1) <- c(nyr*yrint*12,nlat,nlon)
|
|
|
|
+ o.XY[1:(nyr*yrint*12),,] <- TP1
|
|
|
|
+ o.XY[((nyr*yrint*12)+1):((yr2-yr1+lt/12)*12),,] <- ano[1,1,nyr,((yrint*12)+1):lt,nlat:1,]
|
|
|
|
+ if (length(gmt)!=((yr2-yr1+lt/12)*12)) {
|
|
|
|
+ print(paste('GMT time series has not the same length than the array of anomalies.',length(gmt),((yr2-yr1+lt/12)*12)))
|
|
|
|
+ stop()
|
|
|
|
+ }
|
|
|
|
+ # Remove trend from the obs for the whole period.
|
|
|
|
+ o.XY.dt <- array(NA,dim=c(((yr2-yr1+lt/12)*12),nlat,nlon))
|
|
|
|
+ r.XY <- (o.XY.dt[1,,]*NA)
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ if (!is.na(o.XY[1,j,i])) {
|
|
|
|
+ # Regress anomaly onto GMT and remove it from the anomaly.
|
|
|
|
+ a1 <- o.XY[,j,i]
|
|
|
|
+ # Exclude NAs from linear model estimate, but keep them in the residual time series. (!=$residuals)
|
|
|
|
+ fit <- lm(a1~as.vector(gmt),na.action=na.exclude)
|
|
|
|
+ a1.dt <- residuals(fit)
|
|
|
|
+ o.XY.dt[,j,i] <- a1.dt
|
|
|
|
+ # Regression coefficient map.
|
|
|
|
+ r.XY[j,i] <- coefficients(fit)[[2]]
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ # Put detrended obs in verification format.
|
|
|
|
+ ano.dt.ver <- array(NA,dim=c(1,1,nyr,lt,nlat,nlon))
|
|
|
|
+ # Works for forecast start date in each year.
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ d.all <- ts(o.XY.dt[,j,i],deltat=1/12,start=c(yr1,11),end=c((yr2+lt/12),10))
|
|
|
|
+ if (!is.na(o.XY[1,j,i])) {
|
|
|
|
+ for (iy in 1:nyr) {
|
|
|
|
+ d.sl <- window(d.all,start=c((yr1+iy*yrint-yrint),11),end=c((yr1+lt/12+iy*yrint-yrint),10))
|
|
|
|
+ ano.dt.ver[1,1,iy,,((nlat+1)-j),i] <- d.sl
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ return(list(ano.dt=o.XY.dt,ano=o.XY,reg=r.XY,ano.dt.ver=ano.dt.ver,gmt=gmt,tgl.obs=tgl.obs))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_DTO_an <- function(lon,lat,ano,yr1,yr2,yrint,gmt=F) {
|
|
|
|
+#
|
|
|
|
+# Linearly regress out global mean temperature trend (GMT) from observed annual anomalies
|
|
|
|
+# that coincide with hindcasts started on November 1st every 1,2,3, etc. years.
|
|
|
|
+# The anomalies are selected from the first year that coincide with each hindcast.
|
|
|
|
+# NAs are excluded from the linear model fit, but kept in the output.
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Returns linearly detrended and original anomalies in a format to be used for CFU_EOF().
|
|
|
|
+#
|
|
|
|
+# Usage:
|
|
|
|
+#
|
|
|
|
+# CFU_DTO_an(lon,lat,ano,yr1,yr2,yrint,gmt=F)
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ano: Array of anomalies from CFU_ano() with dimensions
|
|
|
|
+# (1, 1, number of start years, forecast years, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# yr1: First start year.
|
|
|
|
+#
|
|
|
|
+# yr2: Last start year.
|
|
|
|
+#
|
|
|
|
+# yrint: Interval between hindcasts in years.
|
|
|
|
+#
|
|
|
|
+# gmt: Global mean temperature time series. (optional)
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# ano.dt.ver: An array of the linearly detrended anomalies in CFU verification format.
|
|
|
|
+#
|
|
|
|
+# ano.dt: An array of the linearly detrended anomalies with dimensions
|
|
|
|
+# (number of years, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# ano: An array of the original input anomalies with dimensions
|
|
|
|
+# (number of years, latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# reg: An array containing a map of regression coefficients
|
|
|
|
+# (number of latitudes, longitudes).
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 15 May 2013.
|
|
|
|
+#
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nyr <- dim(ano)[3]; lt <- dim(ano)[4]
|
|
|
|
+ if (all(gmt==F)) {
|
|
|
|
+ # Global mean temperature GMT from HadCRUT, climate explorer.
|
|
|
|
+ # http://climexp.knmi.nl/data/icrutem3_hadsst2_0-360E_-90-90N_n.dat
|
|
|
|
+ T1 <- read.table("/cfu/data/ukmo/hadcrut/icrutem3_hadsst2_0-360E_-90-90N_n.dat",row.names=1)
|
|
|
|
+ XT <- array(NA,dim=c(12,113))
|
|
|
|
+ for (im in 1:12) {
|
|
|
|
+ XT[im,] <- T1[[im]]
|
|
|
|
+ }
|
|
|
|
+ XT[XT==-9.999000e+02] <- NA
|
|
|
|
+ # GMT anomaly wrt 1356 mo, 113 yr.
|
|
|
|
+ tgl.obs <- ts(c(as.vector(XT),rep(NA,120)),deltat=1/12,start=c(1900,1),end=c(2022,12))
|
|
|
|
+ # GMT anomaly wrt verification period.
|
|
|
|
+ tgl.vp.nov.obs <- window(tgl.obs,start=c(yr1,11),end=c(((yr2+lt)),10))
|
|
|
|
+ tgl.vp.nov.obs <- tgl.vp.nov.obs*((yr2-yr1+lt)/113)
|
|
|
|
+ gmt <- aggregate(tgl.vp.nov.obs,FUN=mean)
|
|
|
|
+ }
|
|
|
|
+ o.XY <- array(NA,dim=c(((yr2-yr1+lt)),nlat,nlon))
|
|
|
|
+ if (yrint>1) {
|
|
|
|
+ TP1 <- ano[1,1,,1:(yrint),nlat:1,]
|
|
|
|
+ TP1 <- aperm(TP1, c(2,1,3,4))
|
|
|
|
+ } else {
|
|
|
|
+ TP1 <- ano[1,1,,1,nlat:1,]
|
|
|
|
+ }
|
|
|
|
+ dim(TP1) <- c(nyr*yrint,nlat,nlon)
|
|
|
|
+ o.XY[1:(nyr*yrint),,] <- TP1
|
|
|
|
+ o.XY[((nyr*yrint)+1):((yr2-yr1+lt)),,] <- ano[1,1,nyr,((yrint)+1):lt,nlat:1,]
|
|
|
|
+ if (length(gmt)!=((yr2-yr1+lt))) {
|
|
|
|
+ print(paste('GMT time series has not the same length than the array of anomalies.',length(gmt),((yr2-yr1+lt))))
|
|
|
|
+ stop()
|
|
|
|
+ }
|
|
|
|
+ # Remove trend from the obs for the whole period.
|
|
|
|
+ o.XY.dt <- array(NA,dim=c(((yr2-yr1+lt)),nlat,nlon))
|
|
|
|
+ r.XY <- (o.XY.dt[1,,]*NA)
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ if (!is.na(o.XY[1,j,i])) {
|
|
|
|
+ # Regress anomaly onto GMT and remove it from the anomaly.
|
|
|
|
+ a1 <- o.XY[,j,i]
|
|
|
|
+ # Exclude NAs from linear model estimate, but keep them in the residual time series. (!=$residuals)
|
|
|
|
+ fit <- lm(a1~as.vector(gmt),na.action=na.exclude)
|
|
|
|
+ a1.dt <- residuals(fit)
|
|
|
|
+ o.XY.dt[,j,i] <- a1.dt
|
|
|
|
+ # Regression coefficient map.
|
|
|
|
+ r.XY[j,i] <- coefficients(fit)[[2]]
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ # Put detrended obs in verification format.
|
|
|
|
+ ano.dt.ver <- array(NA,dim=c(1,1,nyr,lt,nlat,nlon))
|
|
|
|
+ # Works for forecast start date in each year.
|
|
|
|
+ for (i in 1:nlon) {
|
|
|
|
+ for (j in 1:nlat) {
|
|
|
|
+ d.all <- ts(o.XY.dt[,j,i],deltat=1,start=c(yr1,11),end=c((yr2+lt),10))
|
|
|
|
+ if (!is.na(o.XY[1,j,i])) {
|
|
|
|
+ for (iy in 1:nyr) {
|
|
|
|
+ d.sl <- window(d.all,start=c((yr1+iy*yrint-yrint),11),end=c((yr1+lt+iy*yrint-yrint),10))
|
|
|
|
+ ano.dt.ver[1,1,iy,,((nlat+1)-j),i] <- d.sl
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ }
|
|
|
|
+ return(list(ano.dt=o.XY.dt,ano=o.XY,reg=r.XY,ano.dt.ver=ano.dt.ver,gmt=gmt))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_EOF <- function(lon,lat,ano,neofs=15,corr=FALSE) {
|
|
|
|
+#
|
|
|
|
+# Description:
|
|
|
|
+#
|
|
|
|
+# Performs an area-weighted EOF analysis using SVD based on a covariance matrix by default,
|
|
|
|
+# based on a correlation matrix if corr argument is set to TRUE.
|
|
|
|
+#
|
|
|
|
+# Arguments:
|
|
|
|
+#
|
|
|
|
+# lon: A vector of longitudes.
|
|
|
|
+#
|
|
|
|
+# lat: A vector of latitudes.
|
|
|
|
+#
|
|
|
|
+# ano: Array of anomalies with dimensions
|
|
|
|
+# (number of timesteps, number of latitudes, number of longitudes).
|
|
|
|
+#
|
|
|
|
+# neofs: Number of modes to be kept. Default = 15
|
|
|
|
+#
|
|
|
|
+# corr: based on correlation matrix rather than covariance matrix. Default = FALSE
|
|
|
|
+#
|
|
|
|
+# Output:
|
|
|
|
+#
|
|
|
|
+# EOFs: An array of EOF patterns normalized to 1 with dimensions
|
|
|
|
+# (number of modes, number of latitudes, number of longitudes)
|
|
|
|
+#
|
|
|
|
+# PCs: An array of principal components with dimensions
|
|
|
|
+# (number of timesteps, number of modes)
|
|
|
|
+#
|
|
|
|
+# Var: Percentage (%) of variance fraction of total variance explained by each mode (number of modes).
|
|
|
|
+#
|
|
|
|
+# mask: Mask with dimensions (number of latitudes, number of longitudes)
|
|
|
|
+#
|
|
|
|
+# wght: Weights with dimensions (number of latitudes, number of longitudes)
|
|
|
|
+#
|
|
|
|
+# Author:
|
|
|
|
+#
|
|
|
|
+# Fabian Lienert <flienert@ic3.cat> 29 Oct 2012. Inspired by R. Benestad's EOF() in R-package clim.pact.
|
|
|
|
+#
|
|
|
|
+# Lauriane Batte <lauriane.batte@ic3.cat> March 2014. Bug-fixes : 1-reversion of latitudes in the weights
|
|
|
|
+# 2-correlation matrix was used instead of covariance
|
|
|
|
+# 3-double use of the weights
|
|
|
|
+#
|
|
|
|
+# Virginie Guemas <virginie.guemas@ic3.cat> 17 March 2014. Bug-fixes: 1-weight computation - division by sum of cos(lat)
|
|
|
|
+# 2-shuffling of EOFs in EOF.2 intermediate vector
|
|
|
|
+# 3-crash when neofs=1 sorted out
|
|
|
|
+# 4-crash when neofs>nt sorted out
|
|
|
|
+#
|
|
|
|
+# Lauriane Batte <lauriane.batte@ic3.cat> 19 March 2014 : BIG cleanup of code and clarification
|
|
|
|
+# Reduction of the number of transpositions and bug-fixes associated
|
|
|
|
+# Remove of the obsolete LINPACK options
|
|
|
|
+#
|
|
|
|
+# Virginie Guemas <virginie.guemas@ic3.cat> 11 April 2014.
|
|
|
|
+# 1 - Bug-fix in dimensions handling
|
|
|
|
+# EOF composition restitutes now the original field in all cases
|
|
|
|
+# 2 - Simplification of the convention transpose
|
|
|
|
+# 3 - Options to use the correlation matrix rather than the covariance matrix
|
|
|
|
+# 4 - Security checks
|
|
|
|
+# 5 - New normalization of PCs so that PC*EOF only reconstruct the original file
|
|
|
|
+# 6 - Weights = sqrt(cos(lat)) for ano so that covariance matrice weighted by cos(lat)
|
|
|
|
+# 7 - Division of EOF by weights so that the reconstruction is simply EOF * PC
|
|
|
|
+#
|
|
|
|
+ library(s2dverification)
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ # The two lines below have been moved up for the security checks. Virginie
|
|
|
|
+ dim.dat <- dim(ano)
|
|
|
|
+ ny <- dim.dat[2]; nx <- dim.dat[3]; nt <- dim.dat[1]
|
|
|
|
+ # Security check. Virginie
|
|
|
|
+ if ( ny != nlat ) {
|
|
|
|
+ stop("Inconsistent number of latitudes and input field dimensions")
|
|
|
|
+ }
|
|
|
|
+ if ( nx != nlon ) {
|
|
|
|
+ stop("Inconsistent number of longitudes and input field dimensions")
|
|
|
|
+ }
|
|
|
|
+ # Buildup of the mask. Fabian.
|
|
|
|
+ mask <- ano[1,,]; mask[!is.finite(mask)] <- NA
|
|
|
|
+ mask[is.finite(mask)] <- 1
|
|
|
|
+ # Replace mask of NAs with 0s for EOF analysis. Fabian
|
|
|
|
+ ano[!is.finite(ano)] <- 0
|
|
|
|
+ # Area weighting. Weights for EOF. Fabian. Bug below. Lauriane
|
|
|
|
+ # wght <- array(cos(rev(lats)*pi/180),dim=(c(nlat,nlon)))
|
|
|
|
+ # The reversion of latitudes is not suitable. Lauriane
|
|
|
|
+ wght <- array(cos(lat*pi/180),dim=(c(nlat,nlon)))
|
|
|
|
+ # The sum of the weights should equal to 1. Virginie
|
|
|
|
+ wght <- wght/sum(wght,na.rm=T)
|
|
|
|
+ # We want the covariance matrix to be weigthed by the grid cell area
|
|
|
|
+ # so the anomaly field should be weighted by its square root since
|
|
|
|
+ # the covariance matrix equals transpose(ano) times ano. Virginie.
|
|
|
|
+ wght <- sqrt(wght)
|
|
|
|
+ # Fabian was initially dividing by stdv to get the correlation matrix.
|
|
|
|
+ # Weighting of the anomalies that corresponds to the use of a covariance matrix. Lauriane.
|
|
|
|
+ #for (it in 1:nt) {ano[it,,] <- ano[it,,]*wght/stdv}
|
|
|
|
+ ano <- ano*InsertDim(wght,1,nt)
|
|
|
|
+ # The use of the correlation matrix is done under the option corr. Virginie.
|
|
|
|
+ if (corr==TRUE) {
|
|
|
|
+ stdv <- sd(ano,na.rm=T)
|
|
|
|
+ ano <- ano/InsertDim(stdv,1,nt)
|
|
|
|
+ }
|
|
|
|
+ # Time/space matrix for SVD. Fabian
|
|
|
|
+ dim(ano) <- c(nt,ny*nx)
|
|
|
|
+ dim.dat <- dim(ano)
|
|
|
|
+
|
|
|
|
+ # "transposed" means already transposed. Fabian's convention.
|
|
|
|
+ # This is very misleading because the user has actually not transposed its array.
|
|
|
|
+ # Only by chance, the spatial dimension is smaller than the time dimension.
|
|
|
|
+ # I therefore change the convention to:
|
|
|
|
+ # "transpose" means the array needs to be transposed before calling La.svd for
|
|
|
|
+ # computational efficiency because the spatial dimension is larger than the time
|
|
|
|
+ # dimension. This goes with transposing the outputs of LA.svd also. Virginie
|
|
|
|
+ #print(c('dim (nt,npts):',dim.dat))
|
|
|
|
+ #if (dim.dat[2] < dim.dat[1]) {transposed <- TRUE} else {transposed <- FALSE}
|
|
|
|
+ if (dim.dat[2] > dim.dat[1]) {transpose <- TRUE} else {transpose <- FALSE}
|
|
|
|
+ if (transpose) {
|
|
|
|
+ pca <- La.svd(t(ano))
|
|
|
|
+ } else {
|
|
|
|
+ pca <- La.svd(ano)
|
|
|
|
+ }
|
|
|
|
+ # Remove of the obsolete LINPACK cases. Lauriane.
|
|
|
|
+
|
|
|
|
+ # The number of output eofs might be lower than 15 if less than 15 grid points are used. Fabian.
|
|
|
|
+ #if (neofs > dim(ano)[2]) {neofs <- dim(ano)[2]}
|
|
|
|
+ # Whatever the spatial or temporal dimension is smaller than neofs, neofs should be bounded. Virginie.
|
|
|
|
+ neofs=min(dim(ano),neofs)
|
|
|
|
+ # The line below by Fabian is not used. Virginie
|
|
|
|
+ #dim.v <- dim(pca$v); dim.u <- dim(pca$v); dim.x <- dim(ano)
|
|
|
|
+ #
|
|
|
|
+ # Lauriane :
|
|
|
|
+ # La.svd conventions: decomposition X = U D t(V)
|
|
|
|
+ # La.svd$u returns U
|
|
|
|
+ # La.svd$d returns diagonal values of D
|
|
|
|
+ # La.svd$v returns t(V) !!
|
|
|
|
+ # The usual convention is PC=U and EOF=V.
|
|
|
|
+ # If La.svd is called for ano (transpose=FALSE case):
|
|
|
|
+ # EOFs: $v
|
|
|
|
+ # PCs: $u
|
|
|
|
+ # If La.svd is called for t(ano) (transposed=TRUE case):
|
|
|
|
+ # EOFs: t($u)
|
|
|
|
+ # PCs: t($v)
|
|
|
|
+ if (transpose) {
|
|
|
|
+ pca.EOFs <- t(pca$u)
|
|
|
|
+ pca.PCs <- t(pca$v)
|
|
|
|
+ #print('transpose')
|
|
|
|
+ #print(c('Dim(pca.EOFs) :',dim(pca.EOFs)))
|
|
|
|
+ #print(c('Dim(pca.PCs) :',dim(pca.PCs)))
|
|
|
|
+ } else {
|
|
|
|
+ pca.EOFs <- pca$v
|
|
|
|
+ pca.PCs <- pca$u
|
|
|
|
+ #print('not transpose')
|
|
|
|
+ #print(c('Dim(pca.EOFs) :',dim(pca.EOFs)))
|
|
|
|
+ #print(c('Dim(pca.PCs) :',dim(pca.PCs)))
|
|
|
|
+ }
|
|
|
|
+ # The numbers of transposition has been reduced. Virginie
|
|
|
|
+ #PC <- pca.PCs[1:neofs,]
|
|
|
|
+ PC <- pca.PCs[,1:neofs]
|
|
|
|
+ EOF <- pca.EOFs[1:neofs,]
|
|
|
|
+ # The lines seem to be at the origin of the lack of recomposition in last version. Virginie.
|
|
|
|
+ #dim(EOF) <- c(neofs,nx,ny)
|
|
|
|
+ #EOF <- aperm(EOF, c(1,3,2))
|
|
|
|
+ dim(EOF) <- c(neofs,ny,nx)
|
|
|
|
+ # To sort out crash when neofs=1. Virginie
|
|
|
|
+ if (neofs==1) {
|
|
|
|
+ PC=InsertDim(PC,2,1)
|
|
|
|
+ #EOF=InsertDim(EOF,1,1) # This is not necessary with line dim(EOF) above. Lauriane
|
|
|
|
+ }
|
|
|
|
+ # Computation of the % of variance associated with each mode. Fabian.
|
|
|
|
+ W <- pca$d[1:neofs]
|
|
|
|
+ tot.var <- sum(pca$d^2)
|
|
|
|
+ Var.eof <- 100 * pca$d[1:neofs]^2/tot.var
|
|
|
|
+ for (e in 1:neofs) {
|
|
|
|
+ # sg is initialized here. It detects the sign of the pattern and latter the pattern is
|
|
|
|
+ # multiply by -1 if necessary to obtain the nominal IPO or PDO. Fabian
|
|
|
|
+ # CFU_EOF should be generic enough to compute any type of EOF. The change of sign of
|
|
|
|
+ # the pattern should be done outside CFU_EOF. Lauriane & Virginie.
|
|
|
|
+ #sg <- 1
|
|
|
|
+ # Factor to normalize the PC.
|
|
|
|
+ eof.pc.sd <- sd(PC[,e])
|
|
|
|
+ # Factor to normalize the EOF.
|
|
|
|
+ eof.patt.nn <- EOF[e,,]*mask
|
|
|
|
+ # Weights have already been taken into account in covariance matrix. Lauriane.
|
|
|
|
+ #eof.patt.ms <- sum(eof.patt.nn^2*t(wght*mask),na.rm=TRUE)/sum(t(wght*mask),na.rm=TRUE)
|
|
|
|
+ eof.patt.ms <- sum(eof.patt.nn^2,na.rm=TRUE)
|
|
|
|
+ # Normalize PC and EOF
|
|
|
|
+ eof.patt <- eof.patt.nn/eof.patt.ms
|
|
|
|
+ # I would rather multiply eof.pc by eof.patt.ms and W[e] so that the user can reconstruct
|
|
|
|
+ # the original field by EOF * PC only (*weight) rather than EOF * PC * (normalization
|
|
|
|
+ # factor that the user does not know) * (eigenvalue that is not straightforward to get back
|
|
|
|
+ # from the % of variance explained) (*weight). The division by sd(pc) is very easy to do
|
|
|
|
+ # outside CFU_EOF in case the user needs a normalized PC but it is more difficult to get
|
|
|
|
+ # back the amplitude of the variability represented by a given mode outside CFU_EOF if
|
|
|
|
+ # we do not output such a muliplied PC multiplied. Virginie.
|
|
|
|
+ #eof.pc <- PC[,e]/eof.pc.sd
|
|
|
|
+ eof.pc <- PC[,e]*eof.patt.ms*W[e]
|
|
|
|
+ # I am also wondering if we should not output the EOF divided by the weight
|
|
|
|
+ # so that the reconstruction is only EOF * PC only rather than EOF * PC * weight
|
|
|
|
+ # since we have multiplied ano by weight. Virginie
|
|
|
|
+ eof.patt <- eof.patt/wght
|
|
|
|
+ # Change the sign of first EOF pattern for PDO, IPO, etc. Fabian
|
|
|
|
+ #if (e==1) if (eof.patt[round(nlon/2),round(nlat/2)] > 0) sg <- -1
|
|
|
|
+ EOF[e,,] <- eof.patt
|
|
|
|
+ PC[,e] <- eof.pc
|
|
|
|
+ }
|
|
|
|
+ return(list(EOFs=EOF,PCs=PC,Var=Var.eof,mask=mask,wght=wght))
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+
|
|
|
|
+CFU_plotbox <- function(NAO,yr1,yr2,tit='NAO index',tar='DJF',pr.qu=expression('NAO index (PC1) SLP'~(sigma)),obs='ERAInterim',fcsys='EC-Earth 3',file='NAO_fc.ps') {
|
|
|
|
+##
|
|
|
|
+## Produce time series plot showing the distribution of the forecast members with box-and-whisker plot.
|
|
|
|
+##
|
|
|
|
+## Description:
|
|
|
|
+##
|
|
|
|
+## Generates vector graphics file of a forecast box-and-whisker plot time series vs. the observed evolution.
|
|
|
|
+## Only works for re-forecasts started each year.
|
|
|
|
+##
|
|
|
|
+## Usage:
|
|
|
|
+##
|
|
|
|
+## CFU_plotbox(lon,lat,NAO,yr1,yr2,tit='NAO index',tar='DJF',pr.qu='NAO index (PC1) SLP',
|
|
|
|
+## obs='ERAInterim',fcsys='EC-Earth 3',file='NAO_fc.ps')
|
|
|
|
+## Arguments:
|
|
|
|
+##
|
|
|
|
+## NAO: List of obs, forecast arrays of, e.g., the NAO index for each ensemble member, output of CFU_NAO().
|
|
|
|
+## yr1: Year of first forecast target period.
|
|
|
|
+##
|
|
|
|
+## yr2: Year of last forecast target period.
|
|
|
|
+##
|
|
|
|
+## tit: Title of the Figure.
|
|
|
|
+##
|
|
|
|
+## tar: Forecast target period.
|
|
|
|
+##
|
|
|
|
+## pr.qu: Prognostic quantity.
|
|
|
|
+##
|
|
|
|
+## obs: Observational dataset.
|
|
|
|
+##
|
|
|
|
+## fcsys: Climate forecast system.
|
|
|
|
+##
|
|
|
|
+## file: Name of the postscript output file.
|
|
|
|
+##
|
|
|
|
+## Output:
|
|
|
|
+##
|
|
|
|
+## Vector graphics file of a forecast box-and-whisker plot time series vs. the observed evolution.
|
|
|
|
+##
|
|
|
|
+## Author:
|
|
|
|
+##
|
|
|
|
+## Fabian Lienert <flienert@ic3.cat> 09 Aug 2013.
|
|
|
|
+##
|
|
|
|
+## Modifications:
|
|
|
|
+##
|
|
|
|
+## Lauriane Batte <lauriane.batte@ic3.cat> 18 Mar 2014
|
|
|
|
+## 1. No normalization of the indices (already taken care of
|
|
|
|
+## in previous functions)
|
|
|
|
+## 2. Taking into account NA cases (missing data)
|
|
|
|
+##
|
|
|
|
+##
|
|
|
|
+ nyr <- dim(NAO$NAOO.ver)[3]
|
|
|
|
+ postscript(file,horizontal=T)
|
|
|
|
+ ## Observed time series.
|
|
|
|
+ pc.o <- ts(NAO$NAOO.ver[1,1,],deltat=1,start=yr1,end=yr2)
|
|
|
|
+ ## Normalization of obs, forecast members.
|
|
|
|
+ ## This has already been done in CFU_MODE
|
|
|
|
+ # pc.o <- pc.o/sd(pc.o)
|
|
|
|
+ # sd.fc <- apply(NAO$NAOF.ver[1,,],c(1),sd)
|
|
|
|
+ # NAO$NAOF.ver[1,,] <- NAO$NAOF.ver[1,,]/sd.fc
|
|
|
|
+ ## Produce plot.
|
|
|
|
+ par(mar = c(5, 6, 4, 2))
|
|
|
|
+ boxplot(NAO$NAOF.ver[1,,],add=F,main=paste(tit,tar),ylab='',xlab='',
|
|
|
|
+ col='red',lwd=2,t='b',axes=F,cex.main=2,
|
|
|
|
+ ylim=c(-max(abs(c(NAO$NAOF.ver[1,,],pc.o)),na.rm=TRUE),max(abs(c(NAO$NAOF.ver[1,,],pc.o)),na.rm=TRUE)))
|
|
|
|
+ lines(1:nyr,pc.o,lwd=3,col='blue'); abline(h=0,lty=1)
|
|
|
|
+ legend('bottomleft',c(obs,fcsys),lty=c(1,1),lwd=c(3,3),pch=c(NA,NA),
|
|
|
|
+ col=c('blue','red'),horiz=T,bty='n',cex=2.2)
|
|
|
|
+ axis(1,c(1:nyr),NA,cex.axis=2.0)
|
|
|
|
+ axis(1,seq(1,nyr,by=1),seq(start(pc.o)[1],end(pc.o)[1],by=1),cex.axis=2.0)
|
|
|
|
+ mtext(1,line=3,text=tar,cex=1.9)
|
|
|
|
+ mtext(3,line=-2,text=paste(' AC =',round(cor(pc.o,apply(NAO$NAOF.ver[1,,],c(2),mean,na.rm=TRUE)),2)),cex=1.9,adj=0)
|
|
|
|
+ axis(2,cex.axis=2.0)
|
|
|
|
+ mtext(2,line=3,text=pr.qu,cex=1.9)
|
|
|
|
+ box()
|
|
|
|
+ dev.off()
|
|
|
|
+}
|
|
|
|
+
|
|
|
|
+CFU_NAO <- function(lon,lat,ano.f,ano.o,fcavt=2:4,fcsys=1) {
|
|
|
|
+##
|
|
|
|
+## Compute the North Atlantic Oscillation (NAO) index based on the leading EOF of sea level pressure (SLP) anomalies.
|
|
|
|
+## The PCs are obtained by projecting the forecast and oberved anomalies onto the observed EOF pattern.
|
|
|
|
+## By default (fcavt=2:4) CFU_NAO() computes the NAO index for 1-month lead seasonal forecasts that can be
|
|
|
|
+## plotted with CFU_boxplot().
|
|
|
|
+##
|
|
|
|
+## Description:
|
|
|
|
+##
|
|
|
|
+## Returns cross-validated PCs of the NAO index for forecast (ano.f) and observations (ano.o) based on the
|
|
|
|
+## observed leading EOF pattern.
|
|
|
|
+##
|
|
|
|
+## Usage:
|
|
|
|
+##
|
|
|
|
+## CFU_NAO(lon,lat,ano.f,ano.o,fcavt=2:4,fcsys=1)
|
|
|
|
+##
|
|
|
|
+## Arguments:
|
|
|
|
+##
|
|
|
|
+## lon: A vector of longitudes.
|
|
|
|
+##
|
|
|
|
+## lat: A vector of latitudes.
|
|
|
|
+##
|
|
|
|
+## ano.f: Array of North Atlantic SLP (0-60N, 80W-0) forecast anomalies from CFU_ano() with dimensions
|
|
|
|
+## (number of forecast systems, ensemble members, start years, forecast months, latitudes, longitudes).
|
|
|
|
+##
|
|
|
|
+## ano.o: Array of North Atlantic SLP (0-60N, 80W-0) forecast anomalies from CFU_ano() with dimensions
|
|
|
|
+## (number of forecast systems, ensemble members, start years, forecast months, latitudes, longitudes).
|
|
|
|
+##
|
|
|
|
+## fcavt: A vector of time steps to average across defining the target period.
|
|
|
|
+##
|
|
|
|
+## fcsys: Number of the forecast system.
|
|
|
|
+##
|
|
|
|
+##
|
|
|
|
+## Output:
|
|
|
|
+##
|
|
|
|
+## NAOF.ver: Array of forecast NAO index in verification format
|
|
|
|
+## (number of forecast systems, ensemble members, start years, forecast months).
|
|
|
|
+##
|
|
|
|
+## NAOO.ver: Array of observed NAO index in verification format
|
|
|
|
+## (1, 1, number of start years, forecast months).
|
|
|
|
+##
|
|
|
|
+## Author:
|
|
|
|
+##
|
|
|
|
+## Fabian Lienert <flienert@ic3.cat> 08 Aug 2013.
|
|
|
|
+##
|
|
|
|
+## Lauriane Batte <lauriane.batte@ic3.cat> March 2013.
|
|
|
|
+##
|
|
|
|
+## Virginie Guemas <virginie.guemas@ic3.cat> 17 March 2013. Removing the rotation.
|
|
|
|
+##
|
|
|
|
+ nlon <- length(lon); nlat <- length(lat)
|
|
|
|
+ nyr <- dim(ano.f)[3]; nfc <- dim(ano.f)[2]; lt <- dim(ano.f)[4]
|
|
|
|
+ # A security check should be added nlon==nx & nlat==ny. Virginie.
|
|
|
|
+ nlonr <- ceiling(length(lon[lon<180 & lon>=0]))
|
|
|
|
+
|
|
|
|
+ ## Rotate field to start west of 0 longitude. Fabian.
|
|
|
|
+ # This rotation should be done after the projection to output the
|
|
|
|
+ # NAO pattern only. Virginie.
|
|
|
|
+ #ano.o.sh <- (ano.o*NA)
|
|
|
|
+ #ano.o.sh[1,1,,,,] <- ano.o[1,1,,,,c((nlonr+1):nlon,1:nlonr)]
|
|
|
|
+ #ano.f.sh <- (ano.f*NA)
|
|
|
|
+ #ano.f.sh[1,,,,,] <- ano.f[fcsys,,,,,c((nlonr+1):nlon,1:nlonr)]
|
|
|
|
+ #remove(ano.o,ano.f);
|
|
|
|
+ ano.f.sh <- ano.f; ano.o.sh=ano.o
|
|
|
|
+ gc()
|
|
|
|
+
|
|
|
|
+ ## Target period mean. Fabian.
|
|
|
|
+ # The same operations can be done in half the lines using s2dverification. Virginie
|
|
|
|
+ # Also this operation could be more general using CFU_season. Virginie
|
|
|
|
+ o1 <- apply(ano.o.sh[1,1,,fcavt,,],c(1,3,4),mean,na.rm=T)
|
|
|
|
+ dim(o1) <- c(1,1,nyr,1,nlat,nlon)
|
|
|
|
+ f1 <- apply(ano.f.sh[1,,,fcavt,,],c(1,2,4,5),mean,na.rm=T)
|
|
|
|
+ dim(f1) <- c(1,nfc,nyr,1,nlat,nlon)
|
|
|
|
+
|
|
|
|
+ ## Cross-validated PCs. Fabian.
|
|
|
|
+ # This should be extended to nmod and nlt by simple loops. Virginie
|
|
|
|
+ NAOF.ver <- array(NA,c(1,nfc,nyr))
|
|
|
|
+ NAOO.ver <- array(NA,c(1,1,nyr))
|
|
|
|
+ for (iy in 1:nyr) {
|
|
|
|
+ ## Observed EOF excluding one forecast start year. Fabian
|
|
|
|
+ # Without the rotation and only 1 EOF. Virginie
|
|
|
|
+ # OEOF <- CFU_EOF(c(lon[c((nlonr+1):nlon)]-360,lon[1:nlonr]),lat,o1[1,1,-iy,1,,],neofs=3)
|
|
|
|
+ OEOF <- CFU_EOF(lon,lat,o1[1,1,-iy,1,,],neofs=1)
|
|
|
|
+ ## Correct polarity of pattern. Fabian
|
|
|
|
+ sign <- 1; if ( 0 < mean(OEOF$EOFs[1,,round(2*nlat/3)],na.rm=T) ) {sign <- -1}
|
|
|
|
+ OEOF$EOFs <- OEOF$EOFs*sign; OEOF$PCs <- OEOF$PCs*sign
|
|
|
|
+ # There should be an option that allows compute the model EOFs instead of only projecting on the observed EOFS. Virginie
|
|
|
|
+ # for (ifc in 1:nfc) {
|
|
|
|
+ # FEOF <- CFU_EOF(lon,lat,f1[1,ifc,,1,,],neofs=1)
|
|
|
|
+ # NAOF.ver[1,ifc,]=FEOF$PCs
|
|
|
|
+ # }
|
|
|
|
+ ## Project forecast anomalies. Fabian
|
|
|
|
+ # Without the rotation. Virginie
|
|
|
|
+ #PCF <- CFU_MODE(c(lon[c((nlonr+1):nlon)]-360,lon[1:nlonr]),lat,f1,OEOF,mode=1)
|
|
|
|
+ PCF <- CFU_MODE(lon,lat,f1,OEOF,mode=1)
|
|
|
|
+ # for (imo in 1:nmod) { for (il in 1:lt) {NAOF.ver[imo,,iy,il] <- PCF$PC.ver[imo,,iy,imo]}}
|
|
|
|
+ NAOF.ver[1,,iy] <- PCF$PC.ver[1,,iy,1]
|
|
|
|
+ ## Project observed anomalies. Fabian.
|
|
|
|
+ # Without the rotation. Virginie.
|
|
|
|
+ #PCO <- CFU_MODE(c(lon[c((nlonr+1):nlon)]-360,lon[1:nlonr]),lat,o1,OEOF,mode=1)
|
|
|
|
+ PCO <- CFU_MODE(lon,lat,o1,OEOF,mode=1)
|
|
|
|
+ ## Keep PCs of excluded forecast start year. Fabian.
|
|
|
|
+ NAOO.ver[1,1,iy] <- PCO$PC.ver[1,1,iy,1]
|
|
|
|
+ }
|
|
|
|
+ # Here there could be a rotation of the pattern. Virginie
|
|
|
|
+ # OEOF$EOFs[1,,] <- OEOF$EOFs[1,c((nlonr+1):nlon,1:nlonr),]
|
|
|
|
+ # The problem with this rotation is that the user does not have the corresponding longitudes anymore.
|
|
|
|
+ # I don't think if should be rotated at all then. Virginie
|
|
|
|
+ return(list(NAOF.ver=NAOF.ver,NAOO.ver=NAOO.ver,OEOF=OEOF))
|
|
|
|
+}
|