clear all;
close all;

format long
scrnsz = get(0,'ScreenSize');
fontsize = 16;
line_width_plot = 1.5;
set(0,'DefaultAxesFontSize',fontsize);
set(0,'DefaultFigureColor','White');
set(0,'DefaultFigurePosition',[100,100,0.6*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
%set(0,'DefaultFigurePosition',[100,100,0.65*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
%set(0,'DefaultFigurePosition',[100,100,0.65*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
set(0,'DefaultLineLineWidth',line_width_plot);

% Germany data (not part of original data set, downloaded from 
% https://gadm.org/download_country.html )
addpath('/Volumes/DATA3/FESSTHH2020data/gadm40_DEU_shp');

month_to_examine = '08';
num_days_in_month = 31;
missing_data_counts_pa = zeros(num_days_in_month,1);
missing_data_time_indices_pa = cell(num_days_in_month,1);
missing_data_datetimes_pa = cell(num_days_in_month,1);
missing_data_counts_ta = zeros(num_days_in_month,1);
missing_data_time_indices_ta = cell(num_days_in_month,1);
missing_data_datetimes_ta = cell(num_days_in_month,1);
all_missing_datetimes_pa = datetime(1970,1,1,0,0,0,0,'TimeZone','UTC');
all_missing_datetimes_ta = datetime(1970,1,1,0,0,0,0,'TimeZone','UTC');
for day_to_examine_debug = 1:num_days_in_month
    day_to_examine = num2str(day_to_examine_debug);
    while numel(day_to_examine) < 2
        day_to_examine = strcat('0',day_to_examine);
    end
    apollo_folder_name = strcat('/Volumes/DATA3/FESSTHH2020data/FESSTHH2020_',month_to_examine,'_APOLLO_L2/',day_to_examine,'/');
    wxt_folder_name = strcat('/Volumes/DATA3/FESSTHH2020data/FESSTHH2020_',month_to_examine,'_WXT_L2/',day_to_examine,'/');

    % For Figure 9 reproduction:
    stations_to_plot_apollo = {'018LMa','054OGa','082PGa','104PGa'};
    stations_to_plot_wxt = {};
    hour_to_start_plotting = 12;
    minute_to_start_plotting = 30;
    hour_to_end_plotting = 15;
    minute_to_end_plotting = 0;

    % For gridded data:
    hour_to_start_plotting_grid = 12;
    minute_to_start_plotting_grid = 30;
    hour_to_end_plotting_grid = 15;
    minute_to_end_plotting_grid = 0;
    plot_interval_minutes_grid = 5;
    hour_reference_grid = 12;
    minute_reference_grid = 30;

    darkblue = [1 0 137]./255;
    lightblue = [135 206 235]./255;
    redcolor = [220 20 60]./255;
    orange = [255 165 0]./255;
    set(groot,'DefaultAxesColorOrder',[redcolor;orange;lightblue;darkblue])



    apollo_vars = {'pa','ta'};
    num_apollo_vars = numel(apollo_vars);

    wxt_vars = {'hail','precip','hur','pa','ta','wdir','wspeed_max','wspeed'};
    num_wxt_vars = numel(wxt_vars);

    sample_nc_file_apollo = strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc');
    sample_nc_file_wxt = strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc');
    %ncdisp(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc'));
    %ncdisp(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_ta_v00_2020',month_to_examine,day_to_examine,'000000.nc'));

    %ncdisp(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc'));
    %ncdisp(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_ta_v00_2020',month_to_examine,day_to_examine,'000000.nc'));

    %start_datetime = datetime(1970,1,1,0,0,0,0,'TimeZone','UTCLeapSeconds');
    start_datetime = datetime(1970,1,1,0,0,0,0,'TimeZone','UTC');
    % Times
    timedata_apollo = ncread(sample_nc_file_apollo,'time');
    timedata_wxt = ncread(sample_nc_file_wxt,'time');
    % Start and end of time averaging intervals: (2 x num times)
    time_bnds_apollo = ncread(sample_nc_file_apollo,'time_bnds');
    time_bnds_wxt = ncread(sample_nc_file_wxt,'time_bnds');
    % Station identifier code:
    station_id_apollo = ncread(sample_nc_file_apollo,'station_id');
    station_id_wxt = ncread(sample_nc_file_wxt,'station_id');
    % Station name:
    station_name_apollo = ncread(sample_nc_file_apollo,'station_name');
    station_name_wxt = ncread(sample_nc_file_wxt,'station_name');
    % Latitude:
    lat_apollo = ncread(sample_nc_file_apollo,'lat');
    lat_wxt = ncread(sample_nc_file_wxt,'lat');
    % Longitude:
    lon_apollo = ncread(sample_nc_file_apollo,'lon');
    lon_wxt = ncread(sample_nc_file_wxt,'lon');
    % Number of stations:
    num_stations_apollo = numel(lat_apollo);
    num_stations_wxt = numel(lat_wxt);
    % Altitude of ground level relative to mean sea level:
    zsl_apollo = ncread(sample_nc_file_apollo,'zsl');
    zsl_wxt = ncread(sample_nc_file_wxt,'zsl');
    % Height of sensor above ground:
    zag_apollo = ncread(sample_nc_file_apollo,'zag');
    zag_wxt = ncread(sample_nc_file_wxt,'zag');
    % Local climate zone of station: (40 x num_stations)
    lcz_apollo = ncread(sample_nc_file_apollo,'lcz');
    lcz_wxt = ncread(sample_nc_file_wxt,'lcz');


    full_datetimes_apollo = start_datetime + seconds(timedata_apollo);
    full_datetimes_wxt = start_datetime + seconds(timedata_wxt);

    % Read in the meteorological data:

    for apollo_var_num = 1:num_apollo_vars
        apollo_var_name = apollo_vars{1,apollo_var_num};
        switch apollo_var_name
            case 'pa' % Pressure in Pa
                pa_apollo = ncread(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_',apollo_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),apollo_var_name);
            case 'ta' % Temperature in K
                ta_apollo = ncread(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_',apollo_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),apollo_var_name);
            otherwise
                error(strcat('Need to define what to do with the variable named:',32,apollo_var_name));
        end
    end

    % for wxt_var_num = 1:num_wxt_vars
    %     wxt_var_name = wxt_vars{1,wxt_var_num};
    %     switch wxt_var_name
    %         case 'hail' % Amount of hail: Number of hail hits within the measurement interval, m^-2
    %             hail_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'precip' % Accumulated amount of precipitation within measurement interval, kg/m^2
    %             precip_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'hur' % Relative humidity (fraction)
    %             hur_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'pa'
    %             pa_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'ta'
    %             ta_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'wdir' % Wind direction (average of 4-Hz data within measurement interval), degrees
    %             wdir_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'wspeed_max' % Wind gust: maximum 3-s average of 4-Hz data within measurement interval, m/s
    %             wspeed_max_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         case 'wspeed' % Wind speed: average of 4-Hz data within measurement interval, m/s
    %             wspeed_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
    %         otherwise
    %             error(strcat('Need to define what to do with the variable named:',32,wxt_var_name));
    %     end
    % end
    
    % Read in missing data for APOLLO station 49: 077PGa
    missing_data_pa_orig = isnan(pa_apollo(49,:));
    missing_data_pa_small_areas_removed = bwareaopen(missing_data_pa_orig,2);
    missing_data_pa_bw = and(missing_data_pa_orig,~missing_data_pa_small_areas_removed);
    
    missing_data_pa = find(missing_data_pa_bw);
    missing_data_time_indices_pa{day_to_examine_debug,1} = missing_data_pa;
    missing_data_counts_pa(day_to_examine_debug,1) = numel(missing_data_pa);
    missing_data_datetimes_pa{day_to_examine_debug,1} = full_datetimes_apollo(missing_data_pa);
    
    missing_data_ta_orig = isnan(ta_apollo(49,:));
    missing_data_ta_small_areas_removed = bwareaopen(missing_data_ta_orig,2);
    missing_data_ta_bw = and(missing_data_ta_orig,~missing_data_ta_small_areas_removed);
    
    missing_data_ta = find(missing_data_ta_bw);
    missing_data_time_indices_ta{day_to_examine_debug,1} = missing_data_ta;
    missing_data_counts_ta(day_to_examine_debug,1) = numel(missing_data_ta);
    missing_data_datetimes_ta{day_to_examine_debug,1} = full_datetimes_apollo(missing_data_ta);
    date_debug = strcat(month_to_examine,32,day_to_examine)
    if numel(missing_data_pa) < 25
        missing_data_pa = missing_data_pa
    end
    if numel(missing_data_ta) < 25
        missing_data_ta = missing_data_ta
    end
    if numel(missing_data_pa) > 0
        all_missing_datetimes_pa = [all_missing_datetimes_pa;full_datetimes_apollo(missing_data_pa,1)];
    end
    if numel(missing_data_ta) > 0
        all_missing_datetimes_ta = [all_missing_datetimes_ta;full_datetimes_apollo(missing_data_ta,1)];
    end
end
all_missing_datetimes_pa = all_missing_datetimes_pa(2:end);
all_missing_datetimes_ta = all_missing_datetimes_ta(2:end);

% window_plot = 30;
% center_index = 43203;
% indices_plot = (center_index-window_plot):(center_index+window_plot);
% figure;
% subplot(2,1,1)
% plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
% subplot(2,1,2)
% plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')
% 
% %window_plot = 50;
% center_index = 57603;
% indices_plot = (center_index-window_plot):(center_index+window_plot);
% figure;
% subplot(2,1,1)
% plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
% subplot(2,1,2)
% plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')
% 
% %window_plot = 50;
% center_index = 68403;
% indices_plot = (center_index-window_plot):(center_index+window_plot);
% figure;
% subplot(2,1,1)
% plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
% subplot(2,1,2)
% plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')

window_plot = 50;
center_index = 43203+3600;
indices_plot = (center_index-window_plot):(center_index+window_plot);
figure;
subplot(2,1,1)
plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
subplot(2,1,2)
plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')

window_plot = 50;
center_index = 43203+7200;
indices_plot = (center_index-window_plot):(center_index+window_plot);
figure;
subplot(2,1,1)
plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
subplot(2,1,2)
plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')

window_plot = 50;
center_index = 43203+10800;
indices_plot = (center_index-window_plot):(center_index+window_plot);
figure;
subplot(2,1,1)
plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
subplot(2,1,2)
plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')

window_plot = 50;
center_index = 57603+3600;
indices_plot = (center_index-window_plot):(center_index+window_plot);
figure;
subplot(2,1,1)
plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
subplot(2,1,2)
plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')

window_plot = 50;
center_index = 57603+7200;
indices_plot = (center_index-window_plot):(center_index+window_plot);
figure;
subplot(2,1,1)
plot(full_datetimes_apollo(indices_plot),ta_apollo(49,indices_plot),'o-')
subplot(2,1,2)
plot(full_datetimes_apollo(indices_plot),pa_apollo(49,indices_plot),'o-')