
clc; clear all; close all;

%% Load 9-km GPP data (NR and Ops)


fileList = {'../subsets/SMAP_L4_C_mdl_Vv3040_001_alaskaIds2018.h5',
            '../subsets/SMAP_L4_C_mdl_Vv3040_001_NRv41_alaskaIds2018.h5'};
datasetList = {'/GPP/gpp_mean','/SUBSET/year_doy','/SUBSET/site_name'};

dataList = {};
for i = 1:numel(fileList)
    for j = 1:numel(datasetList)
        dataList{i,j} = h5read(fileList{i},datasetList{j});
    end
end

%% Plot 9-km GPP data (NR and Ops)

dataOps = dataList{1}';
dataNR = dataList{2}';

dateOps = double(dataList{3}');
dateNR = double(dataList{4}');
tNR = (1:length(dateNR))';
startOps = find(dateNR(:,1)==dateOps(1,1) & dateNR(:,2)==dateOps(1,2),1);
endNR = find(dateOps(:,1)==dateNR(end,1) & dateOps(:,2)==dateNR(end,2),1);
tOps = (startOps:(startOps+length(dateOps)-1))';

[ yyyymmdd ] = convertDoy2MonthDay( [dateNR ; dateOps(endNR:end,:)] );

t = datenum(yyyymmdd(:,1),yyyymmdd(:,2),yyyymmdd(:,3));

siteName = dataList{5};

for i = 1:15
    figure
    hold on
    plot(t(tNR), dataNR(:,i),'-b')
    plot(t(tOps), dataOps(:,i),'-r')    
    title(siteName(i))
    ylim([0 10])
    ylabel('[g m-2 d-2]')
    xlim([t(1) t(end)])
    datetick('x','yyyy','keeplimits')
    set(gcf,'position',[130         475        1708         420])
end

%% Load 9-km x PFT GPP data 

pftList = {};
for i = 1:8
    pftList{i} = ['/GPP/gpp_pft' num2str(i) '_mean'];
end

dataPftList = {};
for i = 1:numel(fileList)
    for j = 1:numel(pftList)
        dataPftList{i,j} = h5read(fileList{i},pftList{j});
    end
end

%% Plot 9-km x PFT data for specified site

sitePlot = 'BON-SPR'

pftLabel = {'ENF','EBF','DNF','DBF','SHR','GRS','CCP','BCP'};

sIndx = find(strcmp(siteName,sitePlot),1);

figure
plt = []; pftIndx = []; pftCnt = 1;
for i = 1:8
    hold on
    dataPft = dataPftList{2,i}'; % use NR series (1= Ops series)
    dataPft = dataPft(:,sIndx);
    dataPft(dataPft<0) = NaN;
    if sum(~isnan(dataPft))>0
        plt(pftCnt) = plot(t(tNR),dataPft);
        pftIndx(pftCnt) = i;
        pftCnt = pftCnt + 1;
    end
end
ylabel('[g m-2 d-2]')
ylim([0 14])
xlim([t(1) t(end)])
datetick('x','yyyy','keeplimits')
title(sitePlot)
legend(plt,pftLabel(pftIndx));
set(gcf,'position',[130         475        1708         420])

