initializing app sources
This commit is contained in:
parent
be8e826d76
commit
8be5cca7f9
14
README.md
14
README.md
@ -1,16 +1,6 @@
|
||||
--- TO BE REMOVED ---
|
||||
# MIDSTREAM: Model Testing and Forecasting — Official Application Repository
|
||||
|
||||
This is a template for an official application repository. After creating a new repository, fill in the placeholder fields and remove this text.
|
||||
|
||||
✅ Follow [Section 1: Migrating an Application to Git repository](https://docs.cyfronet.pl/x/W4MdDg) when migrating an existing application from `apps_to_migrate`.
|
||||
|
||||
✅ Follow [Section 2: Managing Code Changes in Official Application Repositories](https://docs.cyfronet.pl/x/W4MdDg) to develop or apply updates.
|
||||
|
||||
--- END TO BE REMOVED ---
|
||||
|
||||
# <app_name> — Official Application Repository
|
||||
|
||||
This repository contains the source code and configuration files for the `<app_name>` application used within the [EPISODES Platform](https://EpisodesPlatform.eu/).
|
||||
This repository contains the source code and configuration files for the `MIDSTREAM: Model Testing and Forecasting` application used within the [EPISODES Platform](https://EpisodesPlatform.eu/).
|
||||
|
||||
📦 To test or modify this application in the EPISODES Platform environment, follow the guide:
|
||||
https://docs.cyfronet.pl/display/ISDOC/Editing+application+codes+guide
|
||||
|
||||
66
src/calcrate_exponential.m
Normal file
66
src/calcrate_exponential.m
Normal file
@ -0,0 +1,66 @@
|
||||
function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_exponential(...
|
||||
Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ...
|
||||
Dt_fcst)
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Resample the IR data for considering the Dt for forecasts:
|
||||
|
||||
Ttmp = TEST_IR_time(1);
|
||||
j=1;
|
||||
while Ttmp + Dt_fcst <= TEST_IR_time(end)
|
||||
ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time<Ttmp+Dt_fcst));
|
||||
ind_seism = find(and(TEST_Tev_i>=Ttmp, TEST_Tev_i<Ttmp+Dt_fcst));
|
||||
TEST_IR_time_resampled(j,1) = median(TEST_IR_time(ind_ir_sel));
|
||||
if numel(ind_ir_sel)>0
|
||||
TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel));
|
||||
else
|
||||
TEST_IR_data_resampled(j,1) = 0;
|
||||
end
|
||||
|
||||
TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism));
|
||||
TEST_seism_Nev(j,1) = length(ind_seism);
|
||||
j=j+1;
|
||||
Ttmp = Ttmp + Dt_fcst;
|
||||
clear ind_ir_sel ind_seism
|
||||
end
|
||||
|
||||
%Find IR values = 0 and set them as NaN:
|
||||
clear ind
|
||||
ind = find(TEST_IR_data_resampled<1e-6);
|
||||
TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled);
|
||||
|
||||
[Nsmpls, Npars] = size(Parameter_Samples);
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Calculate seism. rtaes from IR data using the Power Law model
|
||||
clear SeismRate
|
||||
for i=1:Nsmpls
|
||||
%seism rate, Nev/day, rescaled to the forecasting window Dt_fcst:
|
||||
%SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst
|
||||
%indtmp = isnan(SeismRate(:,i));
|
||||
%SeismRate(indtmp,i) = 0;
|
||||
%clear indtmp
|
||||
%CumSeis(:,i) = cumsum(SeismRate(:,i));
|
||||
MU(:,i) = Parameter_Samples(i,1).*(10.^((TEST_IR_data_resampled).*(Parameter_Samples(i,2)))); % days
|
||||
indtmp = isnan(MU(:,i));
|
||||
MU(indtmp,i) = median(MU(~isnan(MU(:,i)),i));
|
||||
%SeismRate(:,i) = Dt_fcst.*(1./MU(:,i));
|
||||
for j=1:numel(TEST_IR_data_resampled)
|
||||
SeismRate(j,i) = Dt_fcst.*(1/median(exprnd(MU(j,i), 100,1)));
|
||||
end
|
||||
clear indtmp mu_tmp
|
||||
CumSeis(:,i) = cumsum(SeismRate(:,i));
|
||||
|
||||
end
|
||||
|
||||
for i=1:numel(TEST_IR_time_resampled)
|
||||
Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)];
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
71
src/calcrate_powerlaw.m
Normal file
71
src/calcrate_powerlaw.m
Normal file
@ -0,0 +1,71 @@
|
||||
function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_powerlaw(...
|
||||
Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ...
|
||||
Dt_fcst)
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Resample the IR data for considering the Dt for forecasts:
|
||||
|
||||
Ttmp = TEST_IR_time(1);
|
||||
j=1;
|
||||
while Ttmp + Dt_fcst <= TEST_IR_time(end)
|
||||
ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time<Ttmp+Dt_fcst));
|
||||
ind_seism = find(and(TEST_Tev_i>=Ttmp, TEST_Tev_i<Ttmp+Dt_fcst));
|
||||
TEST_IR_time_resampled(j,1) = median(TEST_IR_time(ind_ir_sel));
|
||||
if numel(ind_ir_sel)>0
|
||||
TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel));
|
||||
else
|
||||
TEST_IR_data_resampled(j,1) = 0;
|
||||
end
|
||||
|
||||
TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism));
|
||||
TEST_seism_Nev(j,1) = length(ind_seism);
|
||||
j=j+1;
|
||||
Ttmp = Ttmp + Dt_fcst;
|
||||
clear ind_ir_sel ind_seism
|
||||
end
|
||||
|
||||
%Find IR values = 0 and set them as NaN:
|
||||
clear ind
|
||||
ind = find(TEST_IR_data_resampled<1e-6);
|
||||
TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled);
|
||||
|
||||
[Nsmpls, Npars] = size(Parameter_Samples);
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Calculate seism. rtaes from IR data using the Power Law model
|
||||
clear SeismRate
|
||||
% for i=1:Nsmpls
|
||||
% %seism rate, Nev/day, rescaled to the forecasting window Dt_fcst:
|
||||
% SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst
|
||||
% indtmp = isnan(SeismRate(:,i));
|
||||
% SeismRate(indtmp,i) = 0;
|
||||
% clear indtmp
|
||||
% CumSeis(:,i) = cumsum(SeismRate(:,i));
|
||||
% end
|
||||
|
||||
for i=1:Nsmpls
|
||||
%seism rate, Nev/day, rescaled to the forecasting window Dt_fcst:
|
||||
%SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst
|
||||
%MU(:,i) = 10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)); % days
|
||||
MU(:,i) = Parameter_Samples(i,1).*((TEST_IR_data_resampled).^(Parameter_Samples(i,2))); % days
|
||||
indtmp = isnan(MU(:,i));
|
||||
MU(indtmp,i) = median(MU(~isnan(MU(:,i)),i));
|
||||
%SeismRate(:,i) = Dt_fcst.*(1./MU(:,i));
|
||||
for j=1:numel(TEST_IR_data_resampled)
|
||||
SeismRate(j,i) = Dt_fcst.*(1/median(exprnd(MU(j,i), 100,1)));
|
||||
end
|
||||
clear indtmp mu_tmp
|
||||
CumSeis(:,i) = cumsum(SeismRate(:,i));
|
||||
end
|
||||
|
||||
for i=1:numel(TEST_IR_time_resampled)
|
||||
Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)];
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
55
src/calcrate_stationary.m
Normal file
55
src/calcrate_stationary.m
Normal file
@ -0,0 +1,55 @@
|
||||
function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_stationary(...
|
||||
Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ...
|
||||
Dt_fcst)
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Resample the IR data for considering the Dt for forecasts:
|
||||
|
||||
Ttmp = TEST_IR_time(1);
|
||||
j=1;
|
||||
while Ttmp + Dt_fcst <= TEST_IR_time(end)
|
||||
ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time<Ttmp+Dt_fcst));
|
||||
ind_seism = find(and(TEST_Tev_i>=Ttmp, TEST_Tev_i<Ttmp+Dt_fcst));
|
||||
TEST_IR_time_resampled(j,1) = median(TEST_IR_time(ind_ir_sel));
|
||||
if numel(ind_ir_sel)>0
|
||||
TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel));
|
||||
else
|
||||
TEST_IR_data_resampled(j,1) = 0;
|
||||
end
|
||||
TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism));
|
||||
TEST_seism_Nev(j,1) = length(ind_seism);
|
||||
j=j+1;
|
||||
Ttmp = Ttmp + Dt_fcst;
|
||||
clear ind_ir_sel ind_seism
|
||||
end
|
||||
|
||||
%Find IR values = 0 and set them as NaN:
|
||||
clear ind
|
||||
ind = find(TEST_IR_data_resampled<1e-6);
|
||||
TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled);
|
||||
[Nsmpls, Npars] = size(Parameter_Samples);
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Calculate seism. rtaes from IR data using the Power Law model
|
||||
clear SeismRate
|
||||
for i=1:Nsmpls
|
||||
%seism rate, Nev/day, rescaled to the forecasting window Dt_fcst:
|
||||
SeismRate(:,i) = Dt_fcst.* ( (1./(Parameter_Samples(i,1).*(ones(numel(TEST_IR_data_resampled),1)))));
|
||||
indtmp = isnan(SeismRate(:,i));
|
||||
SeismRate(indtmp,i) = median(SeismRate(~isnan(SeismRate(:,i)),i));
|
||||
%SeismRate(indtmp,i) = 0;
|
||||
clear indtmp
|
||||
CumSeis(:,i) = cumsum(SeismRate(:,i));
|
||||
end
|
||||
|
||||
for i=1:numel(TEST_IR_time_resampled)
|
||||
Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)];
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
244
src/midstream_model_testing.m
Normal file
244
src/midstream_model_testing.m
Normal file
@ -0,0 +1,244 @@
|
||||
function[WarningMessages, Exit_ID] = midstream_model_testing(Training_Model_Samples, ...
|
||||
Testing_dataset_for_testing, model, Model_Selection, ir_units, ...
|
||||
iet_units, ir_threshold, iet_threshold, Dt_fcst, FcastTime_window, ...
|
||||
Path_output_results, HighPrctile, LowPrctile, ScaleNumEventForecast, ...
|
||||
Script_Mode)
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Load data
|
||||
disp(['Reading model training parameter samples (model: ' model.folder{Model_Selection} ')'])
|
||||
Parameter_Samples = csvread(Training_Model_Samples, 1);
|
||||
|
||||
disp(['Reading input data saved for testing: ' Testing_dataset_for_testing])
|
||||
%
|
||||
if Script_Mode.ID == 1
|
||||
WarningMessages = {'Using the tool for testing the forecasting performance with the testing dataset',...
|
||||
'List of warning messages:'};
|
||||
disp(cell2mat(WarningMessages(1)))
|
||||
load(Testing_dataset_for_testing);
|
||||
else
|
||||
WarningMessages = {'Using the tool for forecasting seismicity using a custom timeseries of Injection rate data', ...
|
||||
'Error:', ...
|
||||
'This functionality is not yet available! Program Exit without results'};
|
||||
disp(cell2mat(WarningMessages(1)))
|
||||
return
|
||||
end
|
||||
% It loads:
|
||||
% CatMagnitude 1x2 4 char
|
||||
% DelayTime 1x1 8 double
|
||||
% Mc 1x1 8 double
|
||||
% TEST_IETij 93x1 744 double
|
||||
% TEST_IR_data_raw 37077x1 296616 double
|
||||
% TEST_IR_time_raw 37077x1 296616 double
|
||||
% TEST_IRij 93x1 744 double
|
||||
% TEST_Magev_i 94x1 752 double
|
||||
% TEST_Tev_i 94x1 752 double
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% check if events in the test dataset are above Mc:
|
||||
%WarningMessages = {'List of warning messages:'};
|
||||
|
||||
ind_mc = find(TEST_Magev_i<Mc);
|
||||
if numel(ind_mc)>0
|
||||
WarningMessages(numel(WarningMessages)+1) = {'ERROR: events below Mc have been found. Inconsistency in input data. Check'};
|
||||
end
|
||||
clear ind_mc
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% select data above min. IR & IET defined thresholds
|
||||
ind = find(and(TEST_IRij>ir_threshold, TEST_IETij>iet_threshold));
|
||||
%Keep source data in these variables (useful for plots)
|
||||
IR = TEST_IRij(ind);
|
||||
IET = TEST_IETij(ind);
|
||||
clear TEST_IRij TEST_IETij
|
||||
TEST_IRij = IR;
|
||||
TEST_IETij = IET;
|
||||
clear IR IET
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% % Check if the sampling rate of IR data and observed seimsicity rates are
|
||||
% coherent with the Dt time window selected for Seism. rate forecasts.
|
||||
% Eventually write a warning message
|
||||
|
||||
Mean_iet = median(TEST_IETij);
|
||||
Mean_dt_ir = median(diff(TEST_IR_time_raw));
|
||||
|
||||
disp(['Dt - seism. rate forecast (' num2str(Dt_fcst) ' days; ' num2str(Dt_fcst*24) ' hours; ' num2str(Dt_fcst*24*60) ' seconds)'])
|
||||
disp(['Median IET seismicity for testing (' num2str(Mean_iet) ' days; ' num2str(Mean_iet*24) ' hours; ' num2str(Mean_iet*24*60) ' seconds)'])
|
||||
disp(['Median sampling rate of IR data (' num2str(Mean_dt_ir) ' days; ' num2str(Mean_dt_ir*24) ' hours; ' num2str(Mean_dt_ir*24*60) ' seconds)'])
|
||||
|
||||
if Dt_fcst < Mean_iet
|
||||
WarningMessages(numel(WarningMessages)+1) = {'Warning: Dt for forecast is lower than the median IET. Increase Dt.'};
|
||||
elseif Dt_fcst < Mean_dt_ir
|
||||
WarningMessages(numel(WarningMessages)+1) = {'Warning: Dt for forecast is lower than the median Sampling rate if IR data. Increase Dt.'};
|
||||
end
|
||||
|
||||
%start forecast from the time of the firat recorded event in the testing
|
||||
%period
|
||||
To = min(TEST_Tev_i);
|
||||
ind_ir = find(TEST_IR_time_raw >= To);
|
||||
|
||||
TEST_IR_time = TEST_IR_time_raw(ind_ir);
|
||||
TEST_IR_data = TEST_IR_data_raw(ind_ir);
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Seism rate calculation for differnet models:
|
||||
|
||||
|
||||
if Model_Selection== 1 %stationary
|
||||
disp('Testing the stationary model...')
|
||||
[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, ...
|
||||
CumSeis, Nsmpls] = calcrate_stationary(Parameter_Samples, TEST_IR_time, ...
|
||||
TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst);
|
||||
|
||||
elseif Model_Selection== 2 % Power Law
|
||||
disp('Testing the Power-law model...')
|
||||
[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, ...
|
||||
CumSeis, Nsmpls] = calcrate_powerlaw(Parameter_Samples, TEST_IR_time, ...
|
||||
TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst);
|
||||
|
||||
elseif Model_Selection== 3 % Exponential
|
||||
disp('Testing the Exponential model...')
|
||||
[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ...
|
||||
TEST_IR_data_resampled, Forecast_cumulated_Nev, ...
|
||||
CumSeis, Nsmpls] = calcrate_exponential(Parameter_Samples, TEST_IR_time, ...
|
||||
TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst);
|
||||
end
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Plot results, save output data and plots
|
||||
|
||||
%%
|
||||
%Create folder
|
||||
path2SaveOutData = [Path_output_results model.folder{Model_Selection} '/Forecast_Nev_per_' cell2mat(FcastTime_window)];
|
||||
try
|
||||
mkdir(path2SaveOutData)
|
||||
catch
|
||||
disp('Forder for output data already exists')
|
||||
end
|
||||
|
||||
%%
|
||||
figure(1)
|
||||
set(1, 'Position', [943 264 964 670])
|
||||
|
||||
% Times normalized respect to the injection start:
|
||||
To=TEST_IR_time(1);
|
||||
|
||||
for i=1:Nsmpls
|
||||
if i==1
|
||||
plt1 = plot(TEST_IR_time_resampled - To, 1+CumSeis(:,i), '.', 'Color', [0.8 0.8 0.8]);
|
||||
else
|
||||
plot(TEST_IR_time_resampled - To, 1+CumSeis(:,i), '.', 'Color', [0.8 0.8 0.8]);
|
||||
end
|
||||
hold on
|
||||
%pause
|
||||
end
|
||||
|
||||
plt2 = plot(TEST_seism_time - To, cumsum(TEST_seism_Nev), '.r', 'LineWidth', 1.5);
|
||||
plt3 = plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,2), '.-', 'Color', 'k', 'LineWidth', 1.5);
|
||||
plt4 = plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,1), '.-', 'Color', 'b', 'LineWidth', 1.);
|
||||
plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,3), '.-', 'Color', 'b', 'LineWidth', 1.);
|
||||
hold off
|
||||
grid on
|
||||
legend([plt2 plt3 plt4 plt1], 'Cumulated number of events (observed in the testing period)', ...
|
||||
'Cumulated n. events (median forecast)', ...
|
||||
['Cumulated n. events (' num2str(LowPrctile) 'th - ' num2str(HighPrctile) 'th prctiles)'], ...
|
||||
'Cumulated number of events (samples of model parameters)', 'Location', 'southeast');
|
||||
try
|
||||
set(gca, 'YScale', ScaleNumEventForecast)
|
||||
catch
|
||||
%default: log plot
|
||||
set(gca, 'YScale', 'log')
|
||||
end
|
||||
|
||||
xlim([min(TEST_IR_time_resampled)-To max(TEST_IR_time_resampled)-To])
|
||||
xlabel(['Time (days) after ' datestr(To)])
|
||||
ylabel('Cumulated number of events')
|
||||
title({'Comparison between observed and forecast seimsicity', ...
|
||||
['(cumulated number of events considering event rates / ' cell2mat(FcastTime_window) ')']})
|
||||
|
||||
|
||||
%%
|
||||
figure(2)
|
||||
set(2, 'Position', [1321 107 588 846])
|
||||
subplot(3,1,1)
|
||||
plot(TEST_Tev_i-To, TEST_Magev_i, 'o', 'MarkerSize',6 , 'MarkerEdgeColor','k', 'MarkerFaceColor', [0.8 0.8 0.8])
|
||||
hold on
|
||||
plot([min(TEST_Tev_i)-To max(TEST_Tev_i)-To], [Mc Mc], '--r', 'LineWidth', 1.5)
|
||||
xlim([min(TEST_Tev_i)-To max(TEST_Tev_i)-To])
|
||||
grid on
|
||||
xlabel(['Time (days) after ' datestr(To)])
|
||||
ylabel('Magnitude')
|
||||
legend('Seismic event', 'Mc', 'Location', 'northeast')
|
||||
CumNum_observed_events = cumsum(ones(length(TEST_Tev_i),1));
|
||||
% disp(['Time injection start: ' datestr(min(TEST_IR_time_raw))])
|
||||
% disp(['Time first event: ' datestr(min(TEST_Tev_i))])
|
||||
% Tau = min(TEST_Tev_i) -min(TEST_IR_time_raw);
|
||||
% disp(['\tau (Tev-Tinj) = ' num2str(Tau) ' days'])
|
||||
|
||||
subplot(3,1,2)
|
||||
colororder({'k','k'})
|
||||
yyaxis left
|
||||
plot(TEST_Tev_i-To, CumNum_observed_events, 'k', 'LineWidth', 1.5)
|
||||
ylabel({'Cumulated number of recorded', ['events (above Mc= ' num2str(Mc) ')']})
|
||||
grid on
|
||||
yyaxis right
|
||||
plot(TEST_IR_time_raw(ind_ir)-To, TEST_IR_data_raw(ind_ir), 'r', 'LineWidth', 1.5)
|
||||
xlim([min(TEST_Tev_i)-To max(TEST_Tev_i)-To])
|
||||
ylabel(['Injection rate, ' ir_units])
|
||||
xlabel(['Time (days) after ' datestr(To)])
|
||||
legend('cumulated number of events', 'injection rate', 'Location', 'northwest')
|
||||
title('Original data')
|
||||
|
||||
subplot(3,1,3)
|
||||
colororder({'k','k'})
|
||||
yyaxis left
|
||||
plot(TEST_seism_time, cumsum(TEST_seism_Nev), 'k', 'LineWidth', 1.5)
|
||||
xlabel('time')
|
||||
ylabel({'Cumulated number of recorded', ['events (above Mc= ' num2str(Mc) ')']})
|
||||
grid on
|
||||
yyaxis right
|
||||
plot(TEST_IR_time_resampled, TEST_IR_data_resampled, 'r', 'LineWidth', 1.5)
|
||||
ylabel(['Injection rate, ' ir_units])
|
||||
legend('cumulated number of events', 'injection rate')
|
||||
title(['Data averaged in time windows of 1 ' cell2mat(FcastTime_window)])
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% save plot
|
||||
fout1 = [path2SaveOutData '/' model.folder{Model_Selection} '_SeismicityRateForecast_comparison_Nev_per_' cell2mat(FcastTime_window) '.jpg'];
|
||||
saveas(1,fout1,'jpg')
|
||||
fout2 = [path2SaveOutData '/' model.folder{Model_Selection} '_ProcessedData_Nev_per_' cell2mat(FcastTime_window) '.jpg'];
|
||||
saveas(2,fout2,'jpg')
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% Save output parameter samples
|
||||
NomOutSolution_model = [path2SaveOutData '/' model.folder{Model_Selection} '_Seismicity_rate_forecast_Nev_per_' cell2mat(FcastTime_window) '.txt'];
|
||||
cHeader = {['year'] ['month'] ['day'] ['hour'] ['minute'] ['second'] ...
|
||||
['lenght_window(days)'] ['ObservedSeismicity(cumulated)'] ...
|
||||
['SeismicityForecast(cumulated)-' num2str(LowPrctile) 'prctile'] ...
|
||||
['SeismicityForecast(cumulated)-median'] ...
|
||||
['SeismicityForecast(cumulated)-' num2str(HighPrctile) 'prctile']};
|
||||
commaHeader = [cHeader;repmat({','},1,numel(cHeader))];
|
||||
commaHeader = commaHeader(:)';
|
||||
textHeader = cell2mat(commaHeader); %cHeader in text with commas
|
||||
%write header to file
|
||||
fid = fopen(NomOutSolution_model,'w');
|
||||
fprintf(fid,'%s\n',textHeader);
|
||||
fclose(fid);
|
||||
%write data to end of file
|
||||
%Data2Write = [datevec(TEST_IR_time_resampled) repmat(Dt_fcst,numel(TEST_IR_time_resampled),1) cumsum(TEST_seism_Nev) Forecast_cumulated_Nev];
|
||||
dlmwrite(NomOutSolution_model,[datevec(TEST_IR_time_resampled) repmat(Dt_fcst,numel(TEST_IR_time_resampled),1) cumsum(TEST_seism_Nev) Forecast_cumulated_Nev],'-append');
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user