clc
clear all
close all

pft2_file = '/discover/nobackup/jkolassa/restarts/GEOSldas_CNCLM45_NT/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT.catchcnclm51_internal_rst.19800101_0000';
pft4_file = '/discover/nobackup/jkolassa/restarts/GEOSldas_CNCLM45_NT/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT.catchcnclm51_internal_rst.19800101_0000.4PFT';
clm45_file = '/discover/nobackup/jkolassa/restarts/GEOSldas_CNCLM45_NT/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT.catchcnclm45_internal_rst.19800101_0000';
latlon_file = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN45/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN45.tavg24_1d_lnd_Nt.monthly.198001.nc4';

lat_vec = ncread(latlon_file,'lat');
lon_vec = ncread(latlon_file,'lon');

tile_id = ncread(pft2_file,'TILE_ID');

%cnpft_2 = ncread(pft2_file,'CNPFT');
%cnpft_4 = ncread(pft4_file,'CNPFT');
%cnpft_clm45 = ncread(clm45_file,'CNPFT');
%
ityp_2 = ncread(pft2_file,'ITY');
%ityp_4 = ncread(pft4_file,'ITY');
%
fveg_2 = ncread(pft2_file,'FVG');
%fveg_4 = ncread(pft4_file,'FVG');
%fveg_clm45 = ncread(clm45_file,'FVG');
%
%psnsunm_2 = ncread(pft2_file,'PSNSUNM');
%psnsunm_4 = ncread(pft4_file,'PSNSUNM');
%
%nc = 83057;
%var_index = 69;
%
%esai_2 = [];
%esai_4 = [];
%esai_clm45 = [];
%
%for nz = 1:3
%    for nv = 1:2
%
%        new_index = (var_index-1)*6 + (nv-1)*3 + nz;
%        nv_old = (nv-1)*2 + [1 2];
%        old_index = (var_index-1)*12 + (nv_old-1)*3 + nz;
%
%        esai_2 = [esai_2 cnpft_2(nc,new_index)];
%        esai_4 = [esai_4 cnpft_4(nc,old_index)];
%        esai_clm45 = [esai_clm45 cnpft_clm45(nc,old_index)];
%  end 
%end

location_vec = [ 2776 ...
        3359 ...
        3440 ...
        3930 ...
        6705 ...
        7511 ...
        7516 ...
        7869 ...
       10767 ...
       11160 ...
       11638 ...
       12707 ...
       12825 ...
       13062 ...
       13300 ...
       13301 ...
       14713 ...
       15439 ...
       16345 ...
       16461 ...
       17394 ...
       18476 ...
       19730 ...
       20440 ...
       22040 ...
       23341 ];

table_data = readtable('/discover/nobackup/jkolassa/towers_to_tile_num.csv');
global_locations = table_data.tile_id;
global_locations = unique(global_locations);


pft3_vec = find(ityp_2(:,1)==3 & fveg_2(:,1)>=0.7);
pft5_vec = find(ityp_2(:,1)==5 & fveg_2(:,1)>=0.7);
pft9_vec = find(ityp_2(:,1)==9 & fveg_2(:,1)>=0.7);

am_point = [0 -70];
am_dist = sqrt((lat_vec-am_point(1)).^2 + (lon_vec-am_point(2)).^2);
am_location = find(am_dist==min(am_dist));

con_point = [0 25];
con_dist = sqrt((lat_vec-con_point(1)).^2 + (lon_vec-con_point(2)).^2);
con_location = find(con_dist==min(con_dist));


global_locations_extended = [global_locations; pft3_vec(1); ...
                             pft5_vec(1,:); pft9_vec(1,:); ...
                             am_location(1,:); con_location(1,:)];

data_vec = ones(size(global_locations_extended));


% define plot domain
load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);
lat_range = [-60 90];
lon_range = [-180 180];
latlim = [-60 90];
lonlim = [-180 180];

    figure
    ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
    framem
    gridm on;
    setm(ax,'MLineLocation',60,'PLineLocation',60,'MeridianLabel','on',...
         'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',60,...
          'MLabelParallel','south','FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    scatterm(lat_vec(global_locations_extended),lon_vec(global_locations_extended),30,data_vec,"filled")
    caxis([0 1]);
    colormap(jet(2))
    plotm(coastlat, coastlon,'k')
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['Locations'])
    h = colorbar;
    set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,['test'],'png')

