Part_2_Ex_4.m

%% Posterior analysis through Principal Component Analysis
% Let us analyze the spatial and temporal patterns within the DA results through 
% a PCA.
%% Loading the results
% We will load the configuration file we just used to perform DA.

clc, clear all;
addpath(genpath('./Codes/')); dbstop if error;
load("Config_DA.mat");
%% 
% We can use the information stored in config to load the 2D files.

% Load 2D time series
ys = year(datetime(config.run.fromdate));
ye = year(datetime(config.run.todate));
fn = sprintf('%s/COMP_DA_%d_%d_ens.mat',config.summary.directory,ys,ye);
load(fn,'time','TWS_ens'); % This might take a few seconds
%% Perform PCA
% We will now perform the Principal Component Analysis of the spatially distributed 
% TWS time series. More specifically, we will analyze the four dominant Principal 
% Components. From the PCA, we will obtain 4 eigenvectors (spatial patterns) and 
% their corresponding Principal Components (PCs, temporal patterns).

[eigvector, eigvalue, PCAResult] = pcan2(TWS_ens','ReducedDim',4);
% Adjust results to put the dimensionality on the eigenvector instead of
% the PC
eigvector = std(PCAResult).*eigvector;
PCAResult = PCAResult./std(PCAResult);
%% Visualization of mode 1
% Let us now visualize the four dominant spatial and temporal patterns. First, 
% we will load the geographical information that we need for the scatter plot.

% Load latitude and longitude of grid points
% Latlon
fn = sprintf('%s/lonlat_50km_MDB.mat',fileparts(config.DAOptions.OtherData.DesignMatrix));
load(fn);
%% 
% We can now plot the spatial pattern corresponding to the first mode.

mode = 1;
f = figure; f.Position = f.Position + [-f.Position(3)/2 0 f.Position(3) 0]; % Make figure larger
subplot(1,2,1);
scatter(lon_MDB_50km,lat_MDB_50km,24,eigvector(:,mode),'filled');
cb = colorbar; title(sprintf('Spatial pattern - mode %d',mode));
xlabel('Lon'); ylabel('Lat'); 
ch = get(cb,'Title'); set(ch ,'String','TWS (mm)');
%% 
% And we can also represent the spatial pattern corresponding to it

subplot(1,2,2);
plot(time,PCAResult(:,mode));
title(sprintf('Temporal pattern - mode %d',mode));
%% 
% The eigenvector resulted from the PCA analysis can indicate us the percentage 
% of total signal variability that this mode represents.

percent = eigvalue(mode)/sum(eigvalue)*100;
text(datetime(2007,3,1),2,sprintf("Mode 1 variability percentage: %0.0f%%",percent));
print(gcf,'Output_figures/Ex4_mode1','-djpeg');
%% Visualization of modes 2, 3 and 4
% Can you now plot the spatial and temporal patterns for modes 2, 3 and 3?
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%

%% Reflection questions
%% 
% # What is the unit of the temporal pattern (Principal Component) of these 
% modes?
%% 
%% 
% # How is the PCA beneficial with respect to representing sub-basin averaged 
% time series?
%% 
% # What kind of scientific questions can the PCA help to answer?
%%