% Plot average biomass distribution through the area. New approach to
% provide median and standard deviation as well
%

clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inFold = '/home/martin/opendrift/PelagiaIBM/PelagiaMFS09/';
subFold = 'C2k8_1-1/';
Prefix = 'Pelagia_Med_fwd_';
ryears = 2000:2020;
initialMass = 10000;
gshhsFile = '/home/martin/Documents/MATLAB/m_map/private/gshhs_l.b';
EUlatlim = [30 46];
EUlonlim = [-7 36.5];
saveIMG = 1;
IMGfold = 'IMG/biomass_density_plot2/';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mkdir(IMGfold)

Noutputs = 0;
GridLat = EUlatlim(1):0.2:EUlatlim(2)+0.5;
GridLon = EUlonlim(1)-0.5:0.2:EUlonlim(2);
SumBiomass = zeros(length(ryears), length(GridLon), length(GridLat));
n = 1;

for ryear = ryears
    disp(['Reading ',num2str(ryear),' ...'])
    inFile = [inFold,subFold,Prefix,num2str(ryear),'.nc'];
    nc = netcdf(inFile);
    lons = nc{'lon'}(:);
    lats = nc{'lat'}(:);
    biomass = nc{'biomass'}(:);
    close(nc)
    
    % want to use just a subset of timesteps?
    LONS = lons(:,1:1:end);
    LATS = lats(:,1:1:end);
    MASS = biomass(:,1:1:end);
    
    % Limit to europe
    LONS = max(min(LONS(:),EUlonlim(2)),EUlonlim(1));
    LATS = max(min(LATS(:),EUlatlim(2)),EUlatlim(1));

    SumBiomassY = zeros(length(GridLon), length(GridLat));
    
    disp('Calculating biomass distribution ...')
    for nlat = 1:length(GridLat)
        Glat = GridLat(nlat);
        for nlon = 1:length(GridLon)
            Glon = GridLon(nlon);
            idx = LONS >= Glon-0.1 & LONS < Glon+0.1 & LATS >= Glat-0.1 & LATS < Glat+0.1;
            SumBiomassY(nlon,nlat) = sum(MASS(idx));
        end
    end
       
    NoutputsY = size(MASS,2);
    AvgBiomassY = SumBiomassY / NoutputsY / initialMass / (0.2*0.2); % divide by number of outputs and initial mass to get average biomass per square degree

       
    
    SumBiomass(n,:,:) = AvgBiomassY;
    Noutputs = Noutputs + NoutputsY;
    n=n+1;
end


disp('Statistics ...')
AvgBiomass = squeeze(mean(SumBiomass));
StdBiomasa = squeeze(std(SumBiomass));
MedBiomass = squeeze(median(SumBiomass));

% Save it
save('Biomass_distribution.mat','SumBiomass','AvgBiomass','StdBiomasa','MedBiomass');


%% Plot - using matlab mapping tool - average
hFig=figure;
clf
set(hFig, 'Position', [10 10 1000 600]) 
set(hFig,'Color','w')
h = worldmap(EUlatlim,EUlonlim);
setm(h, 'FFaceColor', 'w');   % white background
setm(h,'MapProjection','mercator')
setm(h,'MLineLocation',5)   % interval of lat lon lines (1 deg)
setm(h,'PLineLocation',5)
setm(h,'MLabelLocation',5)  % interval of lat lon labels (1 deg)
setm(h,'PLabelLocation',5)
setm(h,'FontWeight','bold') 
setm(h,'Grid','off')
pcolorm(GridLat-0.1,GridLon-0.1,AvgBiomass')
shading flat
cb = colorbar('east'); 
set(cb,'position',[.9 .245 .03 .615])

% Mediterranean Sea
S = gshhs(gshhsFile, EUlatlim, EUlonlim);
geoshow([S.Lat], [S.Lon])
levels = [S.Level];
land = (levels == 1);
lake = (levels == 2);
island = (levels == 3);
geoshow(S(land),  'FaceColor', [0.5 0.5 0.5])  

setm(h, 'FFaceColor', [239/1000 149/1000 659/1000]);   % blue background
setm(h, 'FEdgeColor', 'k');   % black frame 

% framem('FlineWidth',5,'FEdgeColor','black')

% title('Average biomass','interpreter','none')

if saveIMG
    % save figure
    fOut = [IMGfold,'Pelagia_density_avg_',strip(subFold,'/')];
    export_fig(fOut, '-png', '-pdf');
end

%% Plot - using matlab mapping tool - median
hFig=figure;
clf
set(hFig, 'Position', [10 10 1000 600]) 
set(hFig,'Color','w')
h = worldmap(EUlatlim,EUlonlim);
setm(h, 'FFaceColor', 'w');   % white background
setm(h,'MapProjection','mercator')
setm(h,'MLineLocation',5)   % interval of lat lon lines (1 deg)
setm(h,'PLineLocation',5)
setm(h,'MLabelLocation',5)  % interval of lat lon labels (1 deg)
setm(h,'PLabelLocation',5)
setm(h,'FontWeight','bold') 
setm(h,'Grid','off')
pcolorm(GridLat-0.1,GridLon-0.1,MedBiomass')
% caxis([0 35])
shading flat
cb = colorbar('east'); 
set(cb,'position',[.9 .245 .03 .615])

% Mediterranean Sea
S = gshhs(gshhsFile, EUlatlim, EUlonlim);
geoshow([S.Lat], [S.Lon])
levels = [S.Level];
land = (levels == 1);
lake = (levels == 2);
island = (levels == 3);
geoshow(S(land),  'FaceColor', [0.5 0.5 0.5])  

setm(h, 'FFaceColor', [239/1000 149/1000 659/1000]);   % blue background
setm(h, 'FEdgeColor', 'k');   % black frame 

% framem('FlineWidth',5,'FEdgeColor','black')

% title('Average biomass','interpreter','none')

if saveIMG
    % save figure
    fOut = [IMGfold,'Pelagia_density_median_',strip(subFold,'/')];
    export_fig(fOut, '-png', '-pdf');
end

%% Plot - using matlab mapping tool - standard deviation
hFig=figure;
clf
set(hFig, 'Position', [10 10 1000 600]) 
set(hFig,'Color','w')
h = worldmap(EUlatlim,EUlonlim);
setm(h, 'FFaceColor', 'w');   % white background
setm(h,'MapProjection','mercator')
setm(h,'MLineLocation',5)   % interval of lat lon lines (1 deg)
setm(h,'PLineLocation',5)
setm(h,'MLabelLocation',5)  % interval of lat lon labels (1 deg)
setm(h,'PLabelLocation',5)
setm(h,'FontWeight','bold') 
setm(h,'Grid','off')
pcolorm(GridLat-0.1,GridLon-0.1,StdBiomasa')
% caxis([0 35])
shading flat
cb = colorbar('east'); 
set(cb,'position',[.9 .245 .03 .615])

% Mediterranean Sea
S = gshhs(gshhsFile, EUlatlim, EUlonlim);
geoshow([S.Lat], [S.Lon])
levels = [S.Level];
land = (levels == 1);
lake = (levels == 2);
island = (levels == 3);
geoshow(S(land),  'FaceColor', [0.5 0.5 0.5])  

setm(h, 'FFaceColor', [239/1000 149/1000 659/1000]);   % blue background
setm(h, 'FEdgeColor', 'k');   % black frame 

% framem('FlineWidth',5,'FEdgeColor','black')

% title('Average biomass','interpreter','none')

if saveIMG
    % save figure
    fOut = [IMGfold,'Pelagia_density_std_',strip(subFold,'/')];
    export_fig(fOut, '-png', '-pdf');
end






%% Log Plot - using matlab mapping tool - average
hFig=figure;
clf
set(hFig, 'Position', [10 10 1000 600]) 
set(hFig,'Color','w')
h = worldmap(EUlatlim,EUlonlim);
% setm(h, 'FFaceColor', 'w');   % white background
setm(h,'MapProjection','mercator')
setm(h,'MLineLocation',5)   % interval of lat lon lines (1 deg)
setm(h,'PLineLocation',5)
setm(h,'MLabelLocation',5)  % interval of lat lon labels (1 deg)
setm(h,'PLabelLocation',5)
setm(h,'FontWeight','bold') 
setm(h,'Grid','off')
pcolorm(GridLat+0.5,GridLon+1,AvgBiomass')
pcolorm(GridLat-0.1,GridLon-0.1,AvgBiomass')
shading flat
cb = colorbar('east'); 
set(cb,'position',[.9 .245 .03 .615])
set(gca,'ColorScale','log')

% Mediterranean Sea
S = gshhs(gshhsFile, EUlatlim, EUlonlim);
geoshow([S.Lat], [S.Lon])
levels = [S.Level];
land = (levels == 1);
lake = (levels == 2);
island = (levels == 3);
geoshow(S(land),  'FaceColor', [0.5 0.5 0.5])  

setm(h, 'FFaceColor', [239/1000 149/1000 659/1000]);   % blue background
setm(h, 'FEdgeColor', 'k');   % black frame 

% framem('FlineWidth',5,'FEdgeColor','black')
% title('Average biomass (logarithmic scale)','interpreter','none');

if saveIMG
    % save figure
    fOut = [IMGfold,'Pelagia_density_log10_',strip(subFold,'/')];
    export_fig(fOut, '-png');
end

%% Log Plot - using matlab mapping tool - median
hFig=figure;
clf
set(hFig, 'Position', [10 10 1000 600]) 
set(hFig,'Color','w')
h = worldmap(EUlatlim,EUlonlim);
% setm(h, 'FFaceColor', 'w');   % white background
setm(h,'MapProjection','mercator')
setm(h,'MLineLocation',5)   % interval of lat lon lines (1 deg)
setm(h,'PLineLocation',5)
setm(h,'MLabelLocation',5)  % interval of lat lon labels (1 deg)
setm(h,'PLabelLocation',5)
setm(h,'FontWeight','bold') 
setm(h,'Grid','off')
pcolorm(GridLat+0.5,GridLon+1,MedBiomass')
pcolorm(GridLat-0.1,GridLon-0.1,MedBiomass')
shading flat
cb = colorbar('east'); 
set(cb,'position',[.9 .245 .03 .615])
set(gca,'ColorScale','log')

% Mediterranean Sea
S = gshhs(gshhsFile, EUlatlim, EUlonlim);
geoshow([S.Lat], [S.Lon])
levels = [S.Level];
land = (levels == 1);
lake = (levels == 2);
island = (levels == 3);
geoshow(S(land),  'FaceColor', [0.5 0.5 0.5])  

setm(h, 'FFaceColor', [239/1000 149/1000 659/1000]);   % blue background
setm(h, 'FEdgeColor', 'k');   % black frame 

% framem('FlineWidth',5,'FEdgeColor','black')
% title('Average biomass (logarithmic scale)','interpreter','none');

if saveIMG
    % save figure
    fOut = [IMGfold,'Pelagia_density_median_log10_',strip(subFold,'/')];
    export_fig(fOut, '-png');
end
