Assignment:

1. Read your dataset into "data table" shaped matrices: NCHAN=NLON=144 columns, NSAMP=NT=240 rows. This involves a transpose, as above. Check it with quick contour plots (you don't have to show me the plots, just make sure you are on track). Are you on track?
Yes!

2. Build the terms sequentially from left to right, to decompose your first variable into:
data(lon,mon) = timelonmean + timemean(lon) + meanseasonalcycle1(lon,moy) + anomalies1(lon,mon)
where moy = month of year (1-12, repeating 20 times) and mon = month of data series (1-240). Anomalies(lon, mon) is the residual, meaning it is equal to data minus all the others.
For precip:

For olr:


Confirm that the mean square of all the terms totals up to equal mean(data.*data), meaning that the decomposition is orthogonal.

For precip (pr):
Mean square of all the terms = 35.0759
mean(precip.*precip) = 35.0759
Confirmed!

For olr:
Mean square of all the terms = 58138.1406
mean(olr.*olr) = 58138.1211
Confirmed!
3. Build terms sequentially to decompose your first variable into:
data(lon,mon) = timelonmean + lonmean(mon) + meanseasonalcycle2(lon,moy) + anomalies2(lon,mon).
For precip:

For olr:

Does the mean square of all the terms add up to mean(data.*data)?
For precip (pr):
Mean square of all the terms = 35.5708
mean(precip.*precip) = 35.0759
Yes!

For olr:
Mean square of all the terms = 58151.1094
mean(olr.*olr) = 58138.1211
Yes!
Is meanseasonalcycle1 = meanseasonalcycle2? No, see the plots above.

Is anomalies1 = anomalies2? They are almost equal, but not quite equal. Please see the plots below!:

4. Display images of anomalies(lon,mon) and the covariance matrix C(lon,lon) = anomalies' * anomalies using a symmetric, bloe-red anomaly color table.

Indicate and describe how blobs in C relate to blobs in the anomalies array. ???

5. Display the 3-month lagged covariance matrix C3(lon,lon) = anomalies(lon, 1:237)' * anomalies(lon, 3:240). How do its blobs relate to the blobs in C? What happens at 6 month or 12 month lag? (perhaps make a code loop to display all lags).

As the lag time increases, |C| decreases by some multiple/factor.

6. Concatenate the columns of your two fields datadata = [field1,field2]. Display the 288x288 combined covariance matrix CC([lon,lon],[lon,lon]).


Does some quadrant of it correspond to part 4? Yes, the lower-left quadrant of CC is equal to the negative of C.

Is there symmetry to CC? Yes, the upper-left quadrant and the lower-right quadrant of CC are very similar to each other.

My codes:

clear;
clc;
 
% 1.
precip = nc_varget('Eq_timelon_sections.nc','precip');
olr = nc_varget('Eq_timelon_sections.nc','olr');
 
% 2.
 
'2.'
 
% For precip (pr):
 
'For precip (pr):'
 
pr_timelonmean=nanmean(nanmean(precip,2),1);
pr_timelonmean=repmat(pr_timelonmean,240,144);
 
pr_timemean=nanmean(precip,1);
pr_timemean=repmat(pr_timemean,240,1);
pr_timemean=pr_timemean-pr_timelonmean;
 
pr_meanseasonalcycle1=zeros(12,144);
for moy=1:12
    for y=1:20
        pr_meanseasonalcycle1(moy,:)=pr_meanseasonalcycle1(moy,:)+precip((y-1)*12+moy,:);
    end
end
pr_meanseasonalcycle1=(1/20)*pr_meanseasonalcycle1;
pr_meanseasonalcycle1=repmat(pr_meanseasonalcycle1,20,1);
pr_meanseasonalcycle1=pr_meanseasonalcycle1-pr_timemean-pr_timelonmean;
 
pr_anomalies1=precip-pr_meanseasonalcycle1-pr_timemean-pr_timelonmean;
 
pr_mean_term1Exp2=nanmean(nanmean(pr_timelonmean.*pr_timelonmean+pr_timemean.*pr_timemean...
    +pr_meanseasonalcycle1.*pr_meanseasonalcycle1+pr_anomalies1.*pr_anomalies1));
mean_precipPrecip=nanmean(nanmean(precip.*precip));
 
['Mean square of all the terms = ',num2str(pr_mean_term1Exp2)]
['mean(precip.*precip) = ',num2str(mean_precipPrecip)]
 
lon = nc_varget('Eq_timelon_sections.nc','lon');
time=1:240;
 
pr1=zeros(5,240,144);
pr1(1,:,:)=precip;
pr1(2,:,:)=pr_timelonmean;
pr1(3,:,:)=pr_timemean;
pr1(4,:,:)=pr_meanseasonalcycle1;
pr1(5,:,:)=pr_anomalies1;
pr_titles1={'precip';'pr timelonmean';'pr timemean';'pr meanseasonalcycle1';'pr anomalies1'};
 
figure('Color',[1 1 1]);
for p=1:5
    subplot(1,5,p); imagesc(lon,time,squeeze(pr1(p,:,:))); title({char(pr_titles1(p,1)),'(mm / day)'});
    xlabel('lon (degrees east)'); ylabel('time (month #)'); colorbar; set(gca,'YDir','normal');
end
colormap(redblue);
 
% For olr (olr):
 
'For olr:'
 
olr_timelonmean=nanmean(nanmean(olr,2),1);
olr_timelonmean=repmat(olr_timelonmean,240,144);
 
olr_timemean=nanmean(olr,1);
olr_timemean=repmat(olr_timemean,240,1);
olr_timemean=olr_timemean-olr_timelonmean;
 
olr_meanseasonalcycle1=zeros(12,144);
for moy=1:12
    for y=1:20
        olr_meanseasonalcycle1(moy,:)=olr_meanseasonalcycle1(moy,:)+olr((y-1)*12+moy,:);
    end
end
olr_meanseasonalcycle1=(1/20)*olr_meanseasonalcycle1;
olr_meanseasonalcycle1=repmat(olr_meanseasonalcycle1,20,1);
olr_meanseasonalcycle1=olr_meanseasonalcycle1-olr_timemean-olr_timelonmean;
 
olr_anomalies1=olr-olr_meanseasonalcycle1-olr_timemean-olr_timelonmean;
 
olr_mean_term1Exp2=nanmean(nanmean(olr_timelonmean.*olr_timelonmean+olr_timemean.*olr_timemean...
    +olr_meanseasonalcycle1.*olr_meanseasonalcycle1+olr_anomalies1.*olr_anomalies1));
mean_olrOlr=nanmean(nanmean(olr.*olr));
 
['Mean square of all the terms = ',num2str(olr_mean_term1Exp2)]
['mean(olr.*olr) = ',num2str(mean_olrOlr)]
 
olr1=zeros(5,240,144);
olr1(1,:,:)=olr;
olr1(2,:,:)=olr_timelonmean;
olr1(3,:,:)=olr_timemean;
olr1(4,:,:)=olr_meanseasonalcycle1;
olr1(5,:,:)=olr_anomalies1;
olr_titles1={'olr';'olr timelonmean';'olr timemean';'olr meanseasonalcycle1';'olr anomalies1'};
 
lon = nc_varget('Eq_timelon_sections.nc','lon');
time=1:240;
figure('Color',[1 1 1]);
for p=1:5
    subplot(1,5,p); imagesc(lon,time,squeeze(olr1(p,:,:))); title({char(olr_titles1(p,1)),'(W m^-2)'});
    xlabel('lon (degrees east)'); ylabel('time (month #)'); colorbar; set(gca,'YDir','normal');
end
colormap(redblue);
 
% 3.
 
'3.'
 
% For precip (pr):
 
'For precip (pr):'
 
pr_lonmean=nanmean(precip,2);
pr_lonmean=repmat(pr_lonmean,1,144);
pr_lonmean=pr_lonmean-pr_timelonmean;
 
pr_meanseasonalcycle2=zeros(12,144);
for moy=1:12
    for y=1:20
        pr_meanseasonalcycle2(moy,:)=pr_meanseasonalcycle2(moy,:)+precip((y-1)*12+moy,:);
    end
end
pr_meanseasonalcycle2=(1/20)*pr_meanseasonalcycle2;
pr_meanseasonalcycle2=repmat(pr_meanseasonalcycle2,20,1);
pr_meanseasonalcycle2=pr_meanseasonalcycle2-pr_lonmean-pr_timelonmean;
 
pr_anomalies2=precip-pr_meanseasonalcycle2-pr_lonmean-pr_timelonmean;
 
pr_mean_term2Exp2=nanmean(nanmean(pr_timelonmean.*pr_timelonmean+pr_lonmean.*pr_lonmean...
    +pr_meanseasonalcycle2.*pr_meanseasonalcycle2+pr_anomalies2.*pr_anomalies2));
 
['Mean square of all the terms = ',num2str(pr_mean_term2Exp2)]
['mean(precip.*precip) = ',num2str(mean_precipPrecip)]
 
pr2=zeros(5,240,144);
pr2(1,:,:)=precip;
pr2(2,:,:)=pr_timelonmean;
pr2(3,:,:)=pr_lonmean;
pr2(4,:,:)=pr_meanseasonalcycle2;
pr2(5,:,:)=pr_anomalies2;
pr_titles2={'precip';'pr timelonmean';'pr lonmean';'pr meanseasonalcycle2';'pr anomalies2'};
 
figure('Color',[1 1 1]);
for p=1:5
    subplot(1,5,p); imagesc(lon,time,squeeze(pr2(p,:,:))); title({char(pr_titles2(p,1)),'(mm / day)'});
    xlabel('lon (degrees east)'); ylabel('time (month #)'); colorbar; set(gca,'YDir','normal');
end
colormap(redblue);
 
% For olr:
 
'For olr:'
 
olr_lonmean=nanmean(olr,2);
olr_lonmean=repmat(olr_lonmean,1,144);
olr_lonmean=olr_lonmean-olr_timelonmean;
 
olr_meanseasonalcycle2=zeros(12,144);
for moy=1:12
    for y=1:20
        olr_meanseasonalcycle2(moy,:)=olr_meanseasonalcycle2(moy,:)+olr((y-1)*12+moy,:);
    end
end
olr_meanseasonalcycle2=(1/20)*olr_meanseasonalcycle2;
olr_meanseasonalcycle2=repmat(olr_meanseasonalcycle2,20,1);
olr_meanseasonalcycle2=olr_meanseasonalcycle2-olr_lonmean-olr_timelonmean;
 
olr_anomalies2=olr-olr_meanseasonalcycle2-olr_lonmean-olr_timelonmean;
 
olr_mean_term2Exp2=nanmean(nanmean(olr_timelonmean.*olr_timelonmean+olr_lonmean.*olr_lonmean...
    +olr_meanseasonalcycle2.*olr_meanseasonalcycle2+olr_anomalies2.*olr_anomalies2));
 
['Mean square of all the terms = ',num2str(olr_mean_term2Exp2)]
['mean(olr.*olr) = ',num2str(mean_olrOlr)]
 
olr2=zeros(5,240,144);
olr2(1,:,:)=olr;
olr2(2,:,:)=olr_timelonmean;
olr2(3,:,:)=olr_lonmean;
olr2(4,:,:)=olr_meanseasonalcycle2;
olr2(5,:,:)=olr_anomalies2;
olr_titles2={'olr';'olr timelonmean';'olr lonmean';'olr meanseasonalcycle2';'olr anomalies2'};
 
figure('Color',[1 1 1]);
for p=1:5
    subplot(1,5,p); imagesc(lon,time,squeeze(olr2(p,:,:))); title({char(olr_titles2(p,1)),'(W m^-2)'});
    xlabel('lon (degrees east)'); ylabel('time (month #)'); colorbar; set(gca,'YDir','normal');
end
colormap(redblue);
 
pr_anom_diff=pr_anomalies1-pr_anomalies2; olr_anom_diff=olr_anomalies1-olr_anomalies2;
figure('Color',[1 1 1]); subplot(1,2,1); imagesc(lon,lon,pr_anom_diff); title('pr anom1 - pr anom2');
xlabel('lon (degrees east)'); ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal');
subplot(1,2,2); imagesc(lon,lon,olr_anom_diff); title('olr anom1 - olr anom2');
xlabel('lon (degrees east)'); ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal');
colormap(redblue);
 
% 4.
 
C1=pr_anomalies1' * olr_anomalies1; C2=pr_anomalies2' * olr_anomalies2;
 
figure('Color',[1 1 1]); subplot(1,2,1); imagesc(lon,lon,C1); title('C1');
xlabel('lon (degrees east)'); ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal');
subplot(1,2,2); imagesc(lon,lon,C2); title('C2'); xlabel('lon (degrees east)');
ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal'); colormap(redblue);
 
% 5.
 
C1_lag=zeros(5,144,144); C2_lag=zeros(5,144,144);
l=1;
p=1;
figure('Color',[1 1 1]);
for lag=0:3:12
    C1_lag(l,:,:)=pr_anomalies1(1:240-lag,:)' * olr_anomalies1(lag+1:240,:);
    C2_lag(l,:,:)=pr_anomalies2(1:240-lag,:)' * olr_anomalies2(lag+1:240,:);
 
    subplot(5,2,p); imagesc(lon,lon,squeeze(C1_lag(l,:,:))); title(['C1 ',num2str(lag),'-month lag'],'FontWeight','bold');
    xlabel('lon (degrees east)'); ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal');
    p=p+1;
 
    subplot(5,2,p); imagesc(lon,lon,squeeze(C2_lag(l,:,:))); title(['C2 ',num2str(lag),'-month lag'],'FontWeight','bold');
    xlabel('lon (degrees east)'); ylabel('lon (degrees east)'); colorbar; set(gca,'YDir','normal');
    p=p+1;
 
    l=l+1;
end
colormap(redblue);
 
% 6.
 
aa1=cat(2,pr_anomalies1./pr_timelonmean,olr_anomalies1./olr_timelonmean);
aa2=cat(2,pr_anomalies2./pr_timelonmean,olr_anomalies2./olr_timelonmean);
 
CC1=aa1' * aa1; CC2=aa2' * aa2;
 
figure('Color',[1 1 1]); subplot(1,2,1); imagesc(CC1); title('CC1'); colorbar; set(gca,'YDir','normal');
subplot(1,2,2); imagesc(CC2); title('CC2'); colorbar; set(gca,'YDir','normal'); colormap(redblue);