-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathreconstruction_skill_calc.m
More file actions
executable file
·239 lines (211 loc) · 9.61 KB
/
Copy pathreconstruction_skill_calc.m
File metadata and controls
executable file
·239 lines (211 loc) · 9.61 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
% This script will calculate the range of the error in the calibration windows
% for varying station numbers, using staggered windows. It also outputs a
% mat file containing alot of statistics produced for num_stns and calibration window.
%
% Step 1: Run it with default setup (number of stations = 2:70)
% Step 2: Change NUM_STNS to 1 below in the loop, you also need to remove
% the transpose command here: squeeze(all_stn_precip(n,:,CAL_WDW(c,:)))');
%
% This should take roughly 18 hours to run in its entirety (all window sizes, all regions - run on a desktop i5.
%
% Ryan Batehup, 2014
% Edited for calculation of SAM reconstruction using CPS only
% by Willem Huiskamp, 2015
clear
%% Setup
tic;
load DataFiles/model_output.mat
load site_range.mat
for windowsize = [31]; % The running window in years
for region = [1 2 3 4 5 6]
%% Loading proxies
DIR_NAME = ['Proxies/NoResample/',num2str(windowsize),'yrWindow'];
% script fails at 1, it is noted in the script what to change if you only want one site
NUM_YRS = 500; NUM_TRIALS = 1000;
% Calibration windows set to being 10 overlapping windows over 500 years
NUM_CAL_WDW = 10; clear CAL_WDW;
overlap = ceil(-(NUM_YRS-NUM_CAL_WDW*windowsize)/9.0);
for c=0:9
CAL_WDW(c+1,:) = (1+c*(windowsize-overlap)):((c*(windowsize-overlap))+windowsize);
end
%% Beginning the Loop - extracting the pseudoproxies from the precip file
for c= 1:size(CAL_WDW,1)
numstnstocompare = 2:precip_min(floor(windowsize/30),region);
if max(numstnstocompare) > 70
numstnstocompare = 2:70;
end
for NUM_STNS = numstnstocompare
all_stn_precip=zeros(NUM_TRIALS,NUM_STNS,NUM_YRS);
load([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/',num2str(NUM_STNS),'stns_1000prox.mat']);
% Select region for analysis: 1 = Global (SH really), 2 = South America only, 3 = Australia + New Zealand, 4 = SH without Antarctica
% 5 = Antarctica only, 6 = South Africa only
if region == 1
lat_prox = stn_lat;
lon_prox = stn_lon;
elseif region == 2
lat_prox = stn_lat_SA;
lon_prox = stn_lon_SA;
elseif region == 3
lat_prox = stn_lat_AuNz;
lon_prox = stn_lon_AuNz;
elseif region == 4
lat_prox = stn_lat_noAA;
lon_prox = stn_lon_noAA;
elseif region == 5
lat_prox = stn_lat_AAo;
lon_prox = stn_lon_AAo;
elseif region == 6
lat_prox = stn_lat_SoA;
lon_prox = stn_lon_SoA;
end
stn_precip = nan(NUM_TRIALS,NUM_STNS,NUM_YRS,'single');
for m=1:NUM_TRIALS
for n=1:NUM_STNS
stn_precip(m,n,:) = single(precip_detr(:,lat_prox(m,n),lon_prox(m,n)));
end
end
% Standardise the proxies
stn_precip_mn=mean(stn_precip,3);
stn_precip_std=std(stn_precip,0,3);
for n=1:NUM_TRIALS
for m=1:NUM_STNS
all_stn_precip(n,m,:) = single((stn_precip(n,m,:)-stn_precip_mn(n,m))./(stn_precip_std(n,m)));
end
end
clear lat_prox lon_prox stn_precip stn_precip_mn stn_precip_std
%% Esper et al 2005 CPS Method
for n=1:NUM_TRIALS
corr_matrix = corr(SAM(CAL_WDW(c,:))*ones(1,NUM_STNS), squeeze(all_stn_precip(n,:,CAL_WDW(c,:)))'); % If you want to calculate for one site, change NUM_STNS manually and run this last bit once more.
stn_CPS(n,:) = corr_matrix(1,:)*squeeze(all_stn_precip(n,:,:));
end
% Normalising (it already has mean ~0)
for n=1:NUM_TRIALS
stn_CPS(n,:) = single(squeeze(stn_CPS(n,:))./std(squeeze(stn_CPS(n,:))'));
end
% Skill Evaluation
for n=1:NUM_TRIALS
stn_corr_CPS(n) = single(corr(squeeze(stn_CPS(n,:)'),SAM));
stn_rmse_CPS(n) = single(sqrt(mean((SAM'-squeeze(stn_CPS(n,:))).^2)));
end
all_stn_CPS(NUM_STNS,:,:) = stn_CPS;
all_stn_corr_CPS(NUM_STNS,:) = stn_corr_CPS;
all_stn_rmse_CPS(NUM_STNS,:) = stn_rmse_CPS;
end
numstnstocompare = 2:sat_min(floor(windowsize/30),region);
if max(numstnstocompare) > 70
numstnstocompare = 2:70;
end
for NUM_STNS = numstnstocompare
all_stn_sat=zeros(NUM_TRIALS,NUM_STNS,NUM_YRS);
load([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/',num2str(NUM_STNS),'stns_1000prox.mat']);
% Select region for analysis: 1 = Global (SH really), 2 = South America only, 3 = Australia + New Zealand, 4 = SH without Antarctica
% 5 = Antarctica only, 6 = South Africa only
if region == 1
lat_prox = sat_lat;
lon_prox = sat_lon;
elseif region == 2
lat_prox = sat_lat_SA;
lon_prox = sat_lon_SA;
elseif region == 3
lat_prox = sat_lat_AuNz;
lon_prox = sat_lon_AuNz;
elseif region == 4
lat_prox = sat_lat_noAA;
lon_prox = sat_lon_noAA;
elseif region == 5
lat_prox = sat_lat_AAo;
lon_prox = sat_lon_AAo;
elseif region == 6
lat_prox = sat_lat_SoA;
lon_prox = sat_lon_SoA;
end
stn_sat = nan(NUM_TRIALS,NUM_STNS,NUM_YRS,'single');
for m=1:NUM_TRIALS
for n=1:NUM_STNS
stn_sat(m,n,:) = single(sat_detr(:,lat_prox(m,n),lon_prox(m,n)));
end
end
% Standardise the proxies
stn_sat_mn=mean(stn_sat,3);
stn_sat_std=std(stn_sat,0,3);
for n=1:NUM_TRIALS
for m=1:NUM_STNS
all_stn_sat(n,m,:) = single((stn_sat(n,m,:)-stn_sat_mn(n,m))./(stn_sat_std(n,m)));
end
end
clear lat_prox lon_prox stn_sat stn_sat_mn stn_sat_std
%% Esper et al 2005 CPS Method
for n=1:NUM_TRIALS
corr_matrix_sat = corr(SAM(CAL_WDW(c,:))*ones(1,NUM_STNS), squeeze(all_stn_sat(n,:,CAL_WDW(c,:)))');
sat_CPS(n,:) = corr_matrix_sat(1,:)*squeeze(all_stn_sat(n,:,:));
end
% Normalising (it already has mean ~0)
for n=1:NUM_TRIALS
sat_CPS(n,:) = single(squeeze(sat_CPS(n,:))./std(squeeze(sat_CPS(n,:))'));
end
% Skill Evaluation
for n=1:NUM_TRIALS
sat_corr_CPS(n) = single(corr(squeeze(sat_CPS(n,:)'),SAM));
sat_rmse_CPS(n) = single(sqrt(mean((SAM'-squeeze(sat_CPS(n,:))).^2)));
end
all_sat_CPS(NUM_STNS,:,:) = sat_CPS;
all_sat_corr_CPS(NUM_STNS,:) = sat_corr_CPS;
all_sat_rmse_CPS(NUM_STNS,:) = sat_rmse_CPS;
end
% Save each region to file
if region == 1
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS','all_stn_corr_CPS','all_stn_rmse_CPS','all_sat_CPS','all_sat_corr_CPS','all_sat_rmse_CPS');
elseif region == 2
all_stn_CPS_SA = all_stn_CPS;
all_stn_corr_CPS_SA = all_stn_corr_CPS;
all_stn_rmse_CPS_SA = all_stn_rmse_CPS;
all_sat_CPS_SA = all_sat_CPS;
all_sat_corr_CPS_SA = all_sat_corr_CPS;
all_sat_rmse_CPS_SA = all_sat_rmse_CPS;
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS_SA','all_stn_corr_CPS_SA','all_stn_rmse_CPS_SA','all_sat_CPS_SA','all_sat_corr_CPS_SA','all_sat_rmse_CPS_SA','-append');
elseif region == 3
all_stn_CPS_AuNz = all_stn_CPS;
all_stn_corr_CPS_AuNz = all_stn_corr_CPS;
all_stn_rmse_CPS_AuNz = all_stn_rmse_CPS;
all_sat_CPS_AuNz = all_sat_CPS;
all_sat_corr_CPS_AuNz = all_sat_corr_CPS;
all_sat_rmse_CPS_AuNz = all_sat_rmse_CPS;
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS_AuNz','all_stn_corr_CPS_AuNz','all_stn_rmse_CPS_AuNz','all_sat_CPS_AuNz','all_sat_corr_CPS_AuNz','all_sat_rmse_CPS_AuNz','-append');
elseif region == 4
all_stn_CPS_AA = all_stn_CPS;
all_stn_corr_CPS_AA = all_stn_corr_CPS;
all_stn_rmse_CPS_AA = all_stn_rmse_CPS;
all_sat_CPS_AA = all_sat_CPS;
all_sat_corr_CPS_AA = all_sat_corr_CPS;
all_sat_rmse_CPS_AA = all_sat_rmse_CPS;
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS_AA','all_stn_corr_CPS_AA','all_stn_rmse_CPS_AA','all_sat_CPS_AA','all_sat_corr_CPS_AA','all_sat_rmse_CPS_AA','-append');
elseif region == 5
all_stn_CPS_AAo = all_stn_CPS;
all_stn_corr_CPS_AAo = all_stn_corr_CPS;
all_stn_rmse_CPS_AAo = all_stn_rmse_CPS;
all_sat_CPS_AAo = all_sat_CPS;
all_sat_corr_CPS_AAo = all_sat_corr_CPS;
all_sat_rmse_CPS_AAo = all_sat_rmse_CPS;
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS_AAo','all_stn_corr_CPS_AAo','all_stn_rmse_CPS_AAo','all_sat_CPS_AAo','all_sat_corr_CPS_AAo','all_sat_rmse_CPS_AAo','-append');
elseif region == 6
all_stn_CPS_SoA = all_stn_CPS;
all_stn_corr_CPS_SoA = all_stn_corr_CPS;
all_stn_rmse_CPS_SoA = all_stn_rmse_CPS;
all_sat_CPS_SoA = all_sat_CPS;
all_sat_corr_CPS_SoA = all_sat_corr_CPS;
all_sat_rmse_CPS_SoA = all_sat_rmse_CPS;
save([DIR_NAME,'/CalWdw',num2str(CAL_WDW(c,1)),'_',num2str(CAL_WDW(c,end)),'/tonsofstats.mat'],...
'all_stn_CPS_SoA','all_stn_corr_CPS_SoA','all_stn_rmse_CPS_SoA','all_sat_CPS_SoA','all_sat_corr_CPS_SoA','all_sat_rmse_CPS_SoA','-append');
end
clear all_stn_CPS all_stn_corr_CPS all_stn_rmse_CPS all_sat_CPS all_sat_corr_CPS all_sat_rmse_CPS
c
end
region
end
windowsize
end