clc
clear all
close all

% make selections
storm_name = 'Esther';
start_date = '20200213';
stop_date = '20200301';
lat_range = [-40 0];      % [southern_border northern_border]
lon_range = [110 160];  % [western_border eastern_border]
period_length = 'one_month';

smap_stream = 'Vv4030';
ol_stream = 'OL4001';

clim_file = ['/home/jkolassa/Projects/SMAP_TC/' storm_name '_clim_' period_length '.mat'];
clim_data = load(clim_file);

% convert selections
smap_l4_path = ['/css/smapl4/public/L4_Products/L4_SM/' smap_stream '/gph/'];
smap_ol_path = ['/css/smapl4/public/Model-Only/SoilMoisture/' ol_stream '/gph/'];

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

track_file = ['/home/jkolassa/Projects/SMAP_TC/' 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);

t_vec = {'0130' '0430' '0730' '1030' '1330' '1630' '1930' '2230'};
month_length = [31 28 31 30 31 30 31 31 30 31 30 31];

mat_index = 0;

v_sm = VideoWriter(['/home/jkolassa/Projects/SMAP_TC/' storm_name '_L4_SFMC_anom']);
v_sm.FrameRate = 1;
open(v_sm);

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

            d_str = num2str(d,'%02i');

            for t = 1:length(t_vec)

                mat_index = mat_index + 1;

                smap_l4_file = [smap_l4_path 'Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_gph_' y_str m_str d_str 'T' t_vec{t} '00_' smap_stream '_001.h5'];
                smap_ol_file = [smap_ol_path 'Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_gph_' y_str m_str d_str 'T' t_vec{t} '00_' ol_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);
                   domain_lai_mat = NaN*ones(size(domain_lats,1),size(domain_lats,2),104);
                  % mwrtm_param_file = ['/css/smapl4/public/L4_Products/L4_SM/' smap_stream '/rc_out/Y' y_str '/M11' ...
                  % '/SPL4SM_' smap_stream '.ldas_mwRTMparam.' y_str '1101_0000z.nc4'];
                  % mwrtm_lewt = ncread(mwrtm_param_file,'MWRTM_LEWT');
                end
                sm_l4_data = double(h5read(smap_l4_file,'/Geophysical_Data/sm_surface'));
                sm_ol_data = double(h5read(smap_ol_file,'/Geophysical_Data/sm_surface'));
                lai_l4_data = double(h5read(smap_ol_file,'/Geophysical_Data/leaf_area_index'));
                domain_lai_l4 = lai_l4_data(lon_index,lat_index);
                domain_sm_l4 = sm_l4_data(lon_index,lat_index)-clim_data.clim_sm_l4;
                domain_sm_ol = sm_ol_data(lon_index,lat_index)-clim_data.clim_sm_ol;

                domain_lai_mat(:,:,mat_index) = domain_lai_l4;

                % find nearest track point

                smap_datenum = datenum(y,m,d,str2double(t_vec{t}(1:2)),str2double(t_vec{t}(3:4)),0);
                datenum_dist = sqrt((track_datenum-smap_datenum).^2);
                date_index = find(datenum_dist==min(datenum_dist));
% 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 = flipud(cmap_mat);
                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
                 cmap = cmap_mat;
                 s1 = quantile(domain_sm_l4(:),0.99);
                 s2 = quantile(domain_sm_l4(:),0.01);
                 caxis_min = -1*max(abs(s1),abs(s2));
                 caxis_max = max(abs(s1),abs(s2));
                 caxis_min = -0.1;
                 caxis_max = 0.1

                 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,domain_sm_l4)                                               
                  caxis([caxis_min caxis_max]);                                                                 
                  set(gcf,'Colormap',cmap);                                                                     
                  plotm(coastlat, coastlon,'k')                                                                 
                  hold on                        
                  plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-r','LineWidth',2)
                  hold on                                                               
                  %geoshow(states,'facealpha',0);                                                               
                  set(gca,'XLim',lonlim);                                                                       
                  set(gca,'YLim',latlim);                                                                       
                  tightmap                                                                                      
                  title(['L4 SFMC [m^{3} m^{-3}] ' y_str m_str d_str t_vec{t}])                 
                  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 '_L4_SFMC_anom_' y_str m_str d_str t_vec{t}],'png')

                  frame = getframe(gcf);                                                                        
                  writeVideo(v_sm,frame);


            end

        end

    end
end

%close(v_diff);
close(v_sm);

tile_coord_file = ['/nfs3m/css/curated01/smapl4/data/public/L4_Products/L4_SM/Vv5014/rc_out/SPL4SM_Vv5014.ldas_tilecoord.bin'];
[tile_coord ] = read_tilecoord_le( tile_coord_file, 0 )
avg_lai = nanmean(domain_lai_mat,3);

%for i = 1:size(domain_lai_mat,1)
%    if mod(i,50)==0
%       i
%    end
%    for j = 1:size(domain_lai_mat,2)
%
%        dist_vec = sqrt((tile_coord.com_lat - domain_lats(i,j)).^2 + (tile_coord.com_lon-domain_lons(i,j)).^2);
%        point_index = find(dist_vec==min(dist_vec));
%        mwrtm_lewt_mat(i,j) = mwrtm_lewt(point_index(1));
%    end
%end
%
%avg_VWC = avg_lai.*mwrtm_lewt_mat;
%nan_index = avg_lai == -9999;
%avg_VWC(nan_index) = NaN;
%
%cmap_mat = cmocean('speed',20);
%cmap = cmap_mat;
%                 s1 = quantile(avg_VWC(:),0.99);
%                 s2 = quantile(avg_VWC(:),0.01);
%                 caxis_min = -1*max(abs(s1),abs(s2));
%                 caxis_max = max(abs(s1),abs(s2));
%                 caxis_min = 0;
%                 caxis_max = 10;
%
%                 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_VWC)
%                  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 vegetation water content [kg m^{-2}] '])
%                  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 '_VWC'],'png')
