-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathWRF_ROMS_OBS_data_analysis.m
137 lines (115 loc) · 5.8 KB
/
WRF_ROMS_OBS_data_analysis.m
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
function WRF_ROMS_OBS_data_analysis(strdate,roms_matdir,wrf_matdir)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function WRF_ROMS_OBS_data_analysis(strdate,roms_matdir,wrf_matdir)
%
% Computes basic BIAS/RMSE/CORR statistics of WRF vs observations and ROMS
% vs observations and generates LaTeX table to be included in the Rissaga
% sections of the Rissaga report.
%
% Input arguments:
% strdate: YYYYMMDD
% roms_matdir: directory with ROMS .mat files for that date (generated
% by plot_sealevel_ROMS_BRIFS_OBS.m)
% wrf_matdir: directory with WRF .mat files for that date (generated by
% plot_pressure_WRF_BRIFS_OBS.m)
%
% Author: Baptiste Mourre - SOCIB
% Matjaz Licer - NIB MBS
% Date of creation: Jun-2015
% Last modification: 3-May-2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IFF data exists, perform statistics for CIUTADELLA station:
% strdate='20150730';
% roms_matdir = ['/home/mlicer/BRIFSverif/plots/ROMS/'];
% wrf_matdir = ['/home/mlicer/BRIFSverif/plots/WRF/'];
ROMSmatfile = [roms_matdir 'ROMS_OBS_stations_' strdate '.mat'];
WRFmatfile = [wrf_matdir 'pressure_WRF-OBS_' strdate '.mat'];
load(ROMSmatfile)
load(WRFmatfile)
%% spectral analysis of ROMS SSH and observed SSH:
fft_ROMS_ciutadella = fft_h(P2(:,2));
fft_OBS_ciutadella = fft_h(currentProfilersHF.ciutadella.WTR_PRE);
close all
figure(1);clf
plot(fft_OBS_ciutadella.period,fft_OBS_ciutadella.power,'r',fft_ROMS_ciutadella.period,fft_ROMS_ciutadella.power,'b');
legend('OBS','ROMS','location','northeast','orientation','horizontal')
set(gca,'fontsize',18)
set(gca,'GridLineStyle','--')
title(['SSH signal FFT at CIUTADELLA: ' strdate])
xlabel('FFT period [minutes]')
ylabel('FFT power density [m^2]')
grid on
box on
xlim([0 40])
epsname = [dirname_plots 'ROMS_OBS_SSH_FFT_' strdate '.eps']
print(epsname,'-depsc','-r300')
close all
figure(1);clf
semilogy(fft_OBS_ciutadella.period,fft_OBS_ciutadella.power,'r',fft_ROMS_ciutadella.period,fft_ROMS_ciutadella.power,'b');
legend('OBS','ROMS','location','southwest','orientation','horizontal')
set(gca,'fontsize',18)
set(gca,'GridLineStyle','--')
title(['LOG SSH signal FFT at CIUTADELLA: ' strdate])
xlabel('FFT period [minutes]')
ylabel('Log FFT power density')
grid on
box on
xlim([0 40])
epsname = [dirname_plots 'ROMS_OBS_SSH_FFTlog_' strdate '.eps']
print(epsname,'-depsc','-r300')
ratioOfSSHExtremes = max(P2(:,2))/max(currentProfilersHF.ciutadella.WTR_PRE);
%% WRF ANALYSIS:
% WRF P_stations:
% lon(10)=3.8;lat(10)=39.98; % Ciutadella
% lon(11)=4.4;lat(11)=39.84; % Maó
% lon(12)=2.7;lat(12)=41.1; % Close to Cabo Begur
% lon(13)=3.0;lat(13)=39.4; % Sarapita
% lon(14)=1.3;lat(14)=39.0; % Sant Antoni
% lon(15)=3.1;lat(15)=39.9; % Pollensa
% lon(16)=4.3;lat(16)=39.9; % Lamola
% lon(17)=3.2;lat(17)=39.7; % Colonia St Pere
% lon(18)=2.4;lat(18)=39.5; % Andratx
% CIUTADELLA:
if ~isempty(barometers.ciutadella.AIR_PRE) && any(barometers.ciutadella.AIR_PRE)
WRF_stat_ciutadella = basicStatistics(time,P_stations(:,10),barometers.ciutadella.time,barometers.ciutadella.AIR_PRE)
end
% SANTANTONI:
if ~isempty(barometers.santantoni.AIR_PRE) && any(barometers.santantoni.AIR_PRE)
WRF_stat_santantoni = basicStatistics(time,P_stations(:,14),barometers.santantoni.time,barometers.santantoni.AIR_PRE)
end
%% ROMS ANALYSIS:
if ~isempty(currentProfilersHF.ciutadella.WTR_PRE) && any(currentProfilersHF.ciutadella.WTR_PRE)
ROMS_stat_ciutadella = basicStatistics(time_child,removeROMSLowFrequencies(P2(:,2)),currentProfilersHF.ciutadella.time,currentProfilersHF.ciutadella.WTR_PRE)
end
%% Correlations between WRF PRESSURE ERROR IN SANTANTONI and ROMS SEALEVEL ERROR IN CIUTADELLA:
% figure(1);clf;hold on
% plot(WRF_stat_santantoni.observationTime,WRF_stat_santantoni.modelObservationDifference,'b',...
% ROMS_stat_ciutadella.observationTime,ROMS_stat_ciutadella.modelObservationDifference,'r')
% interpolate both errors to common timeseries:
wrf_error_at_roms_error_times_santantoni = interp1(WRF_stat_santantoni.observationTime,...
WRF_stat_santantoni.modelObservationDifference,ROMS_stat_ciutadella.observationTime);
% compute correlations:
WRF_ROMS_RHO_santantoni = corr(wrf_error_at_roms_error_times_santantoni,ROMS_stat_ciutadella.modelObservationDifference)
%% Correlations between WRF PRESSURE ERROR IN CIUTADELLA and ROMS SEALEVEL ERROR IN CIUTADELLA:
% figure(2);clf;hold on
% plot(WRF_stat_ciutadella.observationTime,WRF_stat_ciutadella.modelObservationDifference,'b',...
% ROMS_stat_ciutadella.observationTime,ROMS_stat_ciutadella.modelObservationDifference,'r')
% interpolate both errors to common timeseries:
wrf_error_at_roms_error_times_ciutadella = interp1(WRF_stat_ciutadella.observationTime,...
WRF_stat_ciutadella.modelObservationDifference,ROMS_stat_ciutadella.observationTime);
% compute correlations:
WRF_ROMS_RHO_ciutadella = corr(wrf_error_at_roms_error_times_ciutadella,ROMS_stat_ciutadella.modelObservationDifference)
% generate LaTeX table with statistics:
fid = fopen(['/home/mlicer/latex/roms_wrf_statistics_' strdate '.tex'],'w')
fprintf(fid,'\\begin{minipage}{\\linewidth}\n\\centering\n')
fprintf(fid,'\\captionof{table}{Model skill on %s.\\newline Ciutadella SSH ratio $SSH_{max}^{ROMS}/SSH_{max}^{OBS}$: %8.4f}\n',strdate,ratioOfSSHExtremes)
fprintf(fid,'\\label{tab::stat_table_%s}\n',strdate)
fprintf(fid,'\\begin{tabular}{lcc}\n')
fprintf(fid,'Model: Field & BIAS & RMSE\\\\\\hline\n')
fprintf(fid,'WRF: MSLP Sant Antoni [hPa] & %8.4f & %8.4f\\\\\n',WRF_stat_santantoni.bias,WRF_stat_santantoni.rmse)
fprintf(fid,'WRF: MSLP Ciutadella [hPa] & %8.4f & %8.4f\\\\\n',WRF_stat_ciutadella.bias,WRF_stat_ciutadella.rmse)
fprintf(fid,'ROMS: SSH Ciutadella [m] & %8.4f & %8.4f\\\\\\hline\n',ROMS_stat_ciutadella.bias,ROMS_stat_ciutadella.rmse)
fprintf(fid,'\\end{tabular}\n')
fprintf(fid,'\\end{minipage}\n')
fclose(fid)