clc
clear all
close all

% make selections
storm_name = 'Gulab';
save_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/' storm_name '/'];

start_date = '20210924';
% Grace: '20210813' 
% Florence: '20180830'; 
% Chaba: '20220628';
% Laura: '20200819'; 
% Ida: '20210820';
% Fred: '20210811';
% Iota: '20201113';  
% Eta :'20201101';   
 % Idai:'20190222';
stop_date = '20211004'; 
%Grace: '20210821';
% Florence: '20180918';
% Chaba: '20220705';
% Laura: '20200829'; 
%Ida : '20210905'; 
%Fred: '20210820'; 
% Iota: '20201119';
% Eta: '20201113';
% stop_date = '20190323'; Idai
lat_range = [14 35]; 
% Grace: [10 55]; 
%Florence: [10 55];
%Chaba: [5 30];
%Fred/Ida, Laura:  [10 55];
%Eta/Iota: [0 33];      % [southern_border northern_border];
%lat_range = [-24 -10];      % Idai
lon_range = [52 97]; 
% Grace : [-110 -55];
%Chaba: [100 130];
%Fred/Ida, Laura: [-110 -50];
% Eta/Iota: [-110 -70];  % [western_border eastern_border]
%lon_range = [25 52];  % Idai

plot_track = 'yes';
track_fmt = 'txt';

smap_stream = 'Vv5030';
ol_stream = 'OL5000';

% convert selections
smap_l4_path = ['/css/smapl4/public/L4_Products/L4_SM/' smap_stream '/aup/'];

start_year = str2double(start_date(1:4));
start_month = str2double(start_date(5:6));
start_day = str2double(start_date(7:8));
stop_year = str2double(stop_date(1:4));
stop_month = str2double(stop_date(5:6));
stop_day = str2double(stop_date(7:8));

% read track data

if (strcmp(plot_track,'yes') & strcmp(track_fmt,'mat'))
   track_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/track_files/' storm_name '_track.mat'];
   track_data = load(track_file);
   track_datenum = track_data.track_mat(:,1);
   track_lat = track_data.track_mat(:,6);
   track_lon = track_data.track_mat(:,7);
elseif (strcmp(plot_track,'yes') & strcmp(track_fmt,'txt'))
   track_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/track_files/' storm_name '_track.ascii'];
   track_data = dlmread(track_file);
   track_datenum = datenum(start_year, track_data(:,1), track_data(:,2), track_data(:,3), 0, 0);
   track_lat = track_data(:,5);
   track_lon = track_data(:,4);
end

t_vec = {'0000' '0300' '0600' '0900' '1200' '1500' '1800' '2100'};
month_length = [31 28 31 30 31 30 31 31 30 31 30 31];

mat_index = 0;
day_index = 0;

%v = VideoWriter(['/home/jkolassa/Projects/SMAP_TC/' storm_name '_SFMC_increments.avi']);
%v.FrameRate = 1;
%open(v);
%
%v_rz = VideoWriter(['/home/jkolassa/Projects/SMAP_TC/' storm_name '_RZMC_increments.avi']);
%v_rz.FrameRate = 1;
%open(v_rz);

for y = start_year:stop_year

    y_str = num2str(y);

    m_first = 1;
    m_last = 12;
    if (y==start_year)
        m_first = start_month;
    end

    if (y==stop_year)
        m_last = stop_month;
    end

    for m = m_first:m_last

        m_str = num2str(m,'%02i');

        day_first = 1;
        day_last = month_length(m);
        if (y==start_year)&(m==start_month)
           day_first = start_day;
        end
        if (y==stop_year)&(m==stop_month)
            day_last = stop_day;
        end

        for d = day_first:day_last
            day_index = day_index + 1;

            d_str = num2str(d,'%02i');
            t_index = 0;
            for t = 1:length(t_vec)

                mat_index = mat_index + 1;
                t_index = t_index + 1;

                smap_l4_file = [smap_l4_path 'Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_aup_' y_str m_str d_str 'T' t_vec{t} '00_' smap_stream '_001.h5'];

                if (mat_index==1)
                   lat_data = double(h5read(smap_l4_file,'/cell_lat'));
                   lon_data = double(h5read(smap_l4_file,'/cell_lon'));
                   lat_index = (lat_data(1,:)>=lat_range(1)) + (lat_data(1,:)<=lat_range(2));
                   lat_index = lat_index == 2;
                   lon_index = (lon_data(:,1)>=lon_range(1)) + (lon_data(:,1)<=lon_range(2));
                   lon_index = lon_index == 2;
                   domain_lats = lat_data(lon_index,lat_index);
                   domain_lons = lon_data(lon_index,lat_index);
                  % sm_diff_mat = NaN*ones(size(domain_lats,1),size(domain_lats,2),104);
                end
                sm_ana_data = double(h5read(smap_l4_file,'/Analysis_Data/sm_surface_analysis'));
                sm_fcst_data = double(h5read(smap_l4_file,'/Forecast_Data/sm_surface_forecast'));
                domain_sm_ana = sm_ana_data(lon_index,lat_index);
                domain_sm_fcst = sm_fcst_data(lon_index,lat_index);

                rz_ana_data = double(h5read(smap_l4_file,'/Analysis_Data/sm_rootzone_analysis'));  
                rz_fcst_data = double(h5read(smap_l4_file,'/Forecast_Data/sm_rootzone_forecast')); 
                domain_rz_ana = rz_ana_data(lon_index,lat_index);                                 
                domain_rz_fcst = rz_fcst_data(lon_index,lat_index);   

                domain_sm_diff(:,:,t_index) = domain_sm_ana - domain_sm_fcst;
                domain_rz_diff(:,:,t_index) = domain_rz_ana - domain_rz_fcst;
  
              %  sm_diff_mat(:,:,mat_index) = domain_sm_diff;
                test_vec(mat_index,1) = max(domain_sm_diff(:));
                test_vec(mat_index,2) = min(domain_sm_diff(:));
                end % t_vec loop

                % compute mean increment

                sm_daily_incr = nansum(domain_sm_diff,3);
                rz_daily_incr = nansum(domain_rz_diff,3); 
                % compute nearest track point
                if strcmp(plot_track,'yes')
                smap_datenum = datenum(y,m,d,0,0,0);            
                datenum_dist = sqrt((track_datenum-smap_datenum).^2);                                           
                date_index = find(datenum_dist==min(datenum_dist));  
                end

                % plot domain-wide surface SM increments

                load coastlines;
                [latcells, loncells] = polysplit(coastlat, coastlon);
                 latlim = lat_range;
                 lonlim = lon_range;

                 cmap_mat = cmocean('balance',20);
                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
                 cmap_mat = flipud(cmap_mat);
                 cmap = cmap_mat;
                 s1 = quantile(sm_daily_incr(:),0.99);
                 s2 = quantile(sm_daily_incr(:),0.01);
                 caxis_min = -1*max(abs(s1),abs(s2));
                 caxis_max = max(abs(s1),abs(s2));
                 caxis_min = -0.06;
                 caxis_max = 0.06;

                 figure
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
                 'FLatLimit',lat_range,'FLonLimit',lon_range);
                 framem off; gridm on;
                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
                 'MLabelParallel','south','FontSize',12)
                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
                  pcolorm(domain_lats,domain_lons,sm_daily_incr)
                  caxis([caxis_min caxis_max]);
                  set(gcf,'Colormap',cmap);
                  plotm(coastlat, coastlon,'k')
                  if strcmp(plot_track,'yes')
                  hold on 
                  plotm(track_lat(1:date_index),track_lon(1:date_index),'-m','LineWidth',2)   
                  end 
                  hold on
                  %geoshow(states,'facealpha',0);
                  set(gca,'XLim',lonlim);
                  set(gca,'YLim',latlim);
                  tightmap
                  title(['SFMC [m^{3} m^{-3}] increment ' y_str m_str d_str])
                  h = colorbar;
                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
                  set(gcf,'paperunits','centimeters')
                  set(gcf,'paperposition',[1 1 17 15])
                  saveas(gcf,[save_path storm_name '_L4_SFMC_increment_' y_str m_str d_str],'png')

                 caxis_min = -0.01;
                 caxis_max = 0.01;
  
                 figure                                                                           
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...          
                 'FLatLimit',lat_range,'FLonLimit',lon_range);                                    
                 framem off; gridm on;                                                            
                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...             
                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...                   
                 'MLabelParallel','south','FontSize',12)                                          
                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                            
                  pcolorm(domain_lats,domain_lons,rz_daily_incr)                                  
                  caxis([caxis_min caxis_max]);                                                   
                  set(gcf,'Colormap',cmap);                                                       
                  plotm(coastlat, coastlon,'k')                                                   
                  if strcmp(plot_track,'yes')
                  hold on
                  plotm(track_lat(1:date_index),track_lon(1:date_index),'-m','LineWidth',2)
                  end   
                  hold on                                                                         
                  %geoshow(states,'facealpha',0);                                                 
                  set(gca,'XLim',lonlim);                                                         
                  set(gca,'YLim',latlim);                                                         
                  tightmap                                                                        
                  title(['RZMC [m^{3} m^{-3}] increment ' y_str m_str d_str])                     
                  h = colorbar;                                                                   
                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);         
                  set(gcf,'paperunits','centimeters')                                             
                  set(gcf,'paperposition',[1 1 17 15])                                            
                  saveas(gcf,[save_path storm_name '_L4_RZMC_increment_' y_str m_str d_str],'png')
%
%                  frame = getframe(gcf);                                                          
%                  writeVideo(v_rz,frame);
        end

    end
end

%close(v);
%close(v_rz);
%!sum_sm_diff = nansum(sm_diff_mat,3);
%!avg_sm_diff = nanmean(sm_diff_mat,3);
%!
%!cmap_mat = cmocean('balance',20);                                                               
%!                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                
%!                 cmap_mat = flipud(cmap_mat);
%!                 cmap = cmap_mat;
%!                 s1 = quantile(sum_sm_diff(:),0.99);                                                          
%!                 s2 = quantile(sum_sm_diff(:),0.01);                                                          
%!                 caxis_min = -1*max(abs(s1),abs(s2));                                                            
%!                 caxis_max = max(abs(s1),abs(s2));   
%!
%!figure                                                                                          
%!                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...                         
%!                 'FLatLimit',lat_range,'FLonLimit',lon_range);                                                   
%!                 framem off; gridm on;
%!                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...                            
%!                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...                                  
%!                 'MLabelParallel','south','FontSize',12)                                                         
%!                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                           
%!                  pcolorm(domain_lats,domain_lons,sum_sm_diff)                                                
%!                  caxis([caxis_min caxis_max]);                                                                  
%!                  set(gcf,'Colormap',cmap);                                                                      
%!                  plotm(coastlat, coastlon,'k')                                                                  
%!                  hold on                                                                                        
%!                  %geoshow(states,'facealpha',0);                                                                
%!                  set(gca,'XLim',lonlim);                                                                        
%!                  set(gca,'YLim',latlim);                                                                        
%!                  tightmap                                                                                       
%!                  title(['total SFMC [m^{3} m^{-3}] increment '])                           
%!                  h = colorbar;                                                                                  
%!                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                        
%!                  set(gcf,'paperunits','centimeters')                                                            
%!                  set(gcf,'paperposition',[1 1 17 15])                                                           
%!                  saveas(gcf,['/home/jkolassa/Projects/SMAP_TC/' storm_name '_SFMC_increment_total'],'png')
%!
%!cmap_mat = cmocean('balance',20);
%!                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
%!                 cmap_mat = flipud(cmap_mat);
%!                 cmap = cmap_mat;
%!                 s1 = quantile(avg_sm_diff(:),0.99);
%!                 s2 = quantile(avg_sm_diff(:),0.01);
%!                 caxis_min = -1*max(abs(s1),abs(s2));
%!                 caxis_max = max(abs(s1),abs(s2));
%!
%!figure
%!                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...                         
%!                 'FLatLimit',lat_range,'FLonLimit',lon_range);
%!                 framem off; gridm on;
%!                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...                            
%!                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...                                  
%!                 'MLabelParallel','south','FontSize',12)
%!                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
%!                  pcolorm(domain_lats,domain_lons,avg_sm_diff)
%!                  caxis([caxis_min caxis_max]);
%!                  set(gcf,'Colormap',cmap);
%!                  plotm(coastlat, coastlon,'k')
%!                  hold on
%!                  %geoshow(states,'facealpha',0);                                                                
%!                  set(gca,'XLim',lonlim);
%!                  set(gca,'YLim',latlim);
%!                  tightmap
%!                  title(['average SFMC [m^{3} m^{-3}] increment '])
%!                  h = colorbar;
%!                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
%!                  set(gcf,'paperunits','centimeters')
%!                  set(gcf,'paperposition',[1 1 17 15])
%!                  saveas(gcf,['/home/jkolassa/Projects/SMAP_TC/' storm_name '_SFMC_increment_avg'],'png')
