clc
clear all
close all

file_1 = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2021/M12/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.202112.nc4';
file_2 = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/202001_to_present/GEOSldas_e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2021/M12/GEOSldas_e0004s_transientCO2_05.tavg24_1d_lnd_Nt.monthly.202112.nc4';

lai_1 = double(ncread(file_1,'CNTLAI'));
lai_2 = double(ncread(file_2,'CNTLAI'));

scuns_1 = double(ncread(file_1,'SCUNS'));
scuns_2 = double(ncread(file_2,'SCUNS'));

lat_1 = double(ncread(file_1,'lat'));
lon_1 = double(ncread(file_1,'lon'));

lat_2 = double(ncread(file_2,'lat'));
lon_2 = double(ncread(file_2,'lon'));

test_lat = unique(lat_1);
test_lon = unique(lon_1);

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

lai_1_mat = NaN*ones(size(lon_mat));
lai_2_mat = NaN*ones(size(lon_mat));

scuns_1_mat = NaN*ones(size(lon_mat));
scuns_2_mat = NaN*ones(size(lon_mat));

for i = 1:length(lon_1)

    lat_index = test_lat == lat_1(i);
    lon_index = test_lon == lon_1(i);

    lai_1_mat(lat_index,lon_index) = lai_1(i);
    lai_2_mat(lat_index,lon_index) = lai_2(i);

end

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(lai_1);
caxis_max = max(lai_1);
clims = [caxis_min caxis_max];

figure
imagesc(lai_1_mat,clims)
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['LAI_1'],'png')

figure
imagesc(lai_2_mat,clims)
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['LAI_2'],'png')


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


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

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(lai_1);
caxis_max = max(lai_1);

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])
pcolorm(lat_mat,lon_mat,lai_1_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('LAI [-]')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['LAI_1'],'png')

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])
pcolorm(lat_mat,lon_mat,lai_2_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('LAI [-]')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['LAI_2'],'png')


