clc
clear all
close all

data_path = '/discover/nobackup/jkolassa/Eunjee_test/';
zth_file = [data_path 'zth_file.txt'];
zenav_file = [data_path 'zenav_file.txt'];

tile_coord_file = [data_path 'IDS_fcst_CTRL_M09_2007_ens1_zth/output/SMAP_EASEv2_M09_GLOBAL/rc_out/IDS_fcst_CTRL_M09_2007_ens1_zth.ldas_tilecoord.bin'];
tile_coord = read_tilecoord_le(tile_coord_file);

zth_data = dlmread(zth_file);
%zth_data = zth_data(:);
zenav_data = dlmread(zenav_file);
%zenav_data = zenav_data(:);

zth_data_new = NaN*ones(size(zth_data,1)*size(zth_data,2),1);
zenav_data_new = NaN*ones(size(zth_data,1)*size(zth_data,2),1);

for i = 1:size(zth_data,1)
     
     if mod(i,10000)==0
        i
     end
     start_index = (i-1)*size(zth_data,2) + 1;
     stop_index = i*size(zth_data,2);
     zth_data_new(start_index:stop_index) = zth_data(i,:)';
     zenav_data_new(start_index:stop_index) = zenav_data(i,:)';
end

time_steps = 2*24*60/7.5;
zth_mat = NaN*ones(tile_coord.N_tile,time_steps);
zenav_mat = NaN*ones(tile_coord.N_tile,time_steps);

for i = 1:time_steps

    start_index = (i-1)*tile_coord.N_tile + 1;
    stop_index = i*tile_coord.N_tile;

    
    zth_mat(:,i) = zth_data_new(start_index:stop_index);
    zenav_mat(:,i) = zenav_data_new(start_index:stop_index);
end

test_lat = unique(tile_coord.com_lat);
test_lon = unique(tile_coord.com_lon);

[lon_mat,lat_mat] = meshgrid(test_lon,test_lat);

zth_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),time_steps);
zenav_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),time_steps);

for i = 1:length(tile_coord.com_lon)

    lat_index = test_lat == tile_coord.com_lat(i);
    lon_index = test_lon == tile_coord.com_lon(i);

    zth_plot_mat(lat_index,lon_index,:) = zth_mat(i,:);
    zenav_plot_mat(lat_index,lon_index,:) = zenav_mat(i,:);
end

v_class = VideoWriter(['/home/jkolassa/Projects/Eunjee_test/class_animation']);
v_class.FrameRate = 1;
open(v_class);

v_tmpreal = VideoWriter(['/home/jkolassa/Projects/Eunjee_test/tmpreal_animation']);
v_tmpreal.FrameRate = 1;
open(v_tmpreal);

for i = 1:time_steps

    class_mat = NaN*ones(size(zth_plot_mat,1),size(zth_plot_mat,2));

    zero_index = zth_plot_mat(:,:,i) > 0;
    class_mat(zero_index) = 1;

    mid_index = (zth_plot_mat(:,:,i)<=0) + (zenav_plot_mat(:,:,i)<=0);
    mid_index = mid_index == 2;
    class_mat(mid_index) = 2;

    else_index = isnan(class_mat);
    class_mat(else_index) = 3;

    tmpreal_mat = zth_plot_mat(:,:,i)./zenav_plot_mat(:,:,i);
    tmpreal_mat(isinf(tmpreal_mat)) = NaN;

    load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);


latlim = [40 60];
lonlim = [40 80];

cmap_mat = cmocean('thermal',3);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = 0.5;
caxis_max = 3.5;

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,...
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off; gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,class_mat)
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('classification')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])

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

cmap_mat = cmocean('tempo',10);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = 0; %min(tmpreal_mat(:));
caxis_max = 1; %max(tmpreal_mat(:));

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,...
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off; gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,tmpreal_mat)
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('classification')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])

frame = getframe(gcf);
writeVideo(v_tmpreal,frame);
end

close(v_class);
close(v_tmpreal)
