clc
clear all
close all

sara_path = '/discover/nobackup/qzhang/LDAS_INPUTS/ZSCOREstats/CF0180x6C_GLOBAL/';
jana_path = '/gpfsm/dnb31/jkolassa/SMAP_TC_project/SMAP_innovations_run/SMAP_omf_C360_ens24/output/CF0360x6C_GLOBAL/stats/z_score_clim/';

int_precision   = 'int32';      % precision of fortran tag
float_precision = 'float32';    % precision of data in input file

grid_coord      = 'latlon_id'; %
data_product    = 'scaling';

N_out_fields = 10;

tile_coord_file = '/gpfsm/dnb31/jkolassa/SMAP_TC_project/SMAP_innovations_run/SMAP_omf_C360_ens24/output/CF0360x6C_GLOBAL/rc_out/SMAP_omf_C360_ens24.ldas_tilecoord.bin';
[tile_coord ] = read_tilecoord_GEOS( tile_coord_file, 0 );


% variables stored in the binary file
field_tags = {  'm_obs_H','s_obs_H','m_mod_H','s_mod_H','N_data_H',...
    'm_obs_V','s_obs_V','m_mod_V','s_mod_V','N_data_V'};

ylabel_tags = { '<Tb_H> [K]','s[Tb_H] [K]','<Tb_H> [K]','s[Tb_H] [K]','N\_data_H',...
    '<Tb_V> [K]','s[Tb_V] [K]','<Tb_V> [K]','s[Tb_V] [K]','N\_data_V'};

angles = [40];

for p = 1:73

    p_str = num2str(p,'%02i');

    sara_file = [sara_path 'SMAPc180e24_zscore_stats_2015_p19_2020_p18_hscale_0.00_W_9p_Nmin_20_D_p' p_str '.bin'];
    jana_file = [jana_path 'SMAPc360e24_zscore_stats_2015_p1_2021_p73_hscale_0.00_W_9p_Nmin_20_D_p' p_str '.bin'];

    ifp = fopen( sara_file, 'r', 'b' );
    [data_sara, inc_angle_sara, col_ind_sara, row_ind_sara, asc_flag_sara, Nmin_sara, O_start_time_sara, O_end_time_sara, N_grid_sara, tile_id_sara ] = ...
    read_SMOS_prep(sara_file, N_out_fields, grid_coord, data_product, ifp);

    ifp = fopen( jana_file, 'r', 'b' );
    [data_jana, inc_angle_jana, col_ind_jana, row_ind_jana, asc_flag_jana, Nmin_jana, O_start_time_jana, O_end_time_jana, N_grid_jana, tile_id_jana ] = ... 
    read_SMOS_prep(jana_file, N_out_fields, grid_coord, data_product, ifp);

    if (p == 1)
       m_obs_h_sara = NaN*ones(size(data_sara,2),73);
       m_mod_h_sara = NaN*ones(size(data_sara,2),73);
       m_obs_v_sara = NaN*ones(size(data_sara,2),73);
       m_mod_v_sara = NaN*ones(size(data_sara,2),73);

       m_obs_h_jana = NaN*ones(size(data_jana,2),73);
       m_mod_h_jana = NaN*ones(size(data_jana,2),73);
       m_obs_v_jana = NaN*ones(size(data_jana,2),73);
       m_mod_v_jana = NaN*ones(size(data_jana,2),73);
    end 
    m_obs_h_sara(:,p) = data_sara(1,:);
    m_mod_h_sara(:,p) = data_sara(3,:);
    m_obs_v_sara(:,p) = data_sara(6,:);
    m_mod_v_sara(:,p) = data_sara(8,:);

    m_obs_h_jana(:,p) = data_jana(1,:);
    m_mod_h_jana(:,p) = data_jana(3,:);
    m_obs_v_jana(:,p) = data_jana(6,:);
    m_mod_v_jana(:,p) = data_jana(8,:);
end

m_obs_h_sara(m_obs_h_sara==-9999) = NaN;
m_mod_h_sara(m_mod_h_sara==-9999) = NaN;
m_obs_v_sara(m_obs_v_sara==-9999) = NaN;
m_mod_v_sara(m_mod_v_sara==-9999) = NaN;

m_obs_h_jana(m_obs_h_jana==-9999) = NaN;
m_mod_h_jana(m_mod_h_jana==-9999) = NaN;
m_obs_v_jana(m_obs_v_jana==-9999) = NaN;
m_mod_v_jana(m_mod_v_jana==-9999) = NaN;

%test_lat = unique(tile_coord.com_lat);
%test_lon = unique(tile_coord.com_lon);
%
%[lon_mat,lat_mat] = meshgrid(test_lon,test_lat);
%
%mod_obs_h_jana_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),size(mod_obs_h_jana,2));
%
%for i = 1:length(lon_data)
%
%    lat_index = test_lat == lat_data(i);
%    lon_index = test_lon == lon_data(i);
%
%    mod_obs_h_jana_mat(lat_index,lon_index,:) = m_obs_h_jana(i,:);
%
%end

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


latlim = [-60 90];
lonlim = [-180 180];

for p = 1:73

 p_str = num2str(p,'%02i');

tmp = m_obs_h_jana(:,p);

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(tmp(:),0.05);
s2 = quantile(tmp(:),0.95);
caxis_min = min(mod_obs_h_jana(:));
caxis_max = max(mod_obs_h_jana(:));

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                                             
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
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])
scatterm(tile_coord.com_lat,tile_coord.com_lon,5,m_obs_h_jana(:,p),'s','filled')
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(['observation mean Tb_{H} [K] for pentad ' p_str])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/test_readzscore/m_obs_h_jana_' p_str],'png')

end
