% Count number of particles and biomass inside a polygon (NAdriatic)
%
% 28.7.2023 - Martin

% close all
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inFold = '/home/martin/opendrift/PelagiaMFS09/';
subFold = 'C2k8_1-1/';
fPrefix = 'Pelagia_Med_fwd_';
% % Kind of N Adri
% PY = [43.9, 44.2, 45.9, 45.9, 43.9];
% PX = [11, 17, 17 , 10, 11];
% % Ancona-Zadar
PA(1).name = 'NAdri';
PA(1).PY = [43.50941394323307, 44.31755748599436, 45.9, 45.9, 43.50941394323307];
PA(1).PX = [13.1815090706803 , 15.879736836008636, 17 , 9, 13.1815090706803];
PA(1).Ecoregion_code = 1;
% % Pescara-Sibenik
% PY = [42.37832612003335, 43.883142056295846, 45.9, 45.9, 42.37832612003335];
% PX = [14.064512455701623, 16.056077443878248, 17 , 9, 14.064512455701623];

% % S Adri
SAdri = [43.70384736866087, 16.64251295690252; 41.33293832103632, 20.863039648129355;
40.56538290490766, 19.78780395897874; 42.71133262928798, 15.373439726329124;
43.70384736866087, 16.64251295690252];
PA(2).PY = SAdri(:,1);
PA(2).PX = SAdri(:,2);
PA(2).name = 'SAdri';
PA(2).Ecoregion_code = 3;
% Tunisia
Tun = [36.64003271392411, 8.706574837914962;
38.27800597251262, 8.718503398515272;
38.32949019580724, 11.623107904691125;
36.57300318947219, 11.485929457787543;
36.64003271392411, 8.706574837914962];
PA(3).PY = Tun(:,1);
PA(3).PX = Tun(:,2);
PA(3).name = 'Tunisia';
PA(3).Ecoregion_code = 5;

% Pelagia sightings index
sightFile = 'Jelly_NAdriatic_26.1.2024.xls';
sheetName = 'Pelagia';
sightFileS = 'Index_Pelagia_data_ALL_MV.csv';

saveIMG = 1;
IMGfold = 'IMG/Area_density_plot/';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mkdir(IMGfold);



%     % Plot area
%     EUlatlim = [30 46];
%     EUlonlim = [-7 36.5];
%     figure
%     m_proj('mercator' ,'lon',[EUlonlim(1) EUlonlim(2)], ...
%     'lat',[EUlatlim(1) EUlatlim(2)]);
%     m_grid('xtick',[EUlonlim(1):4:EUlonlim(2)],'ytick',[EUlatlim(1):2:EUlatlim(2)],'linestyle','none')
%     hold on
% %     m_plot(lons(:),lats(:),'.','Color',[0.5, 0.5, 0.5])
% %     m_plot(lons(:,end),lats(:,end),'*')
%     m_plot(PX,PY,'k')
%     m_plot(SaX,SaY,'b')
%     m_plot(TunX,TunY,'r')
%     hold off
%     m_gshhs_l('patch',[.7 .7 .7],'edgecolor','black')
%     title(inFiles(nf).name)
%     pause 


% Jelly sightings - Sonia's file
PS = readtable(sightFileS);
PS.dtime = datetime(PS.Year, PS.Month, 15);


inFiles = dir([inFold,subFold,fPrefix,'*.nc']);

for nf = 1:length(inFiles)
    inFile = [inFiles(nf).folder,'/',inFiles(nf).name];
    nc = netcdf(inFile);
    nctime = nc{'time'}(:);  % seconds since 1970-01-01 00:00:00
    biomass = nc{'biomass'}(:);  % Actually biomass of each particlez
    lons = nc{'lon'}(:);
    lats = nc{'lat'}(:);
    close(nc)

    dtime = datetime(1970,1,1) + seconds(nctime);
    % Count only for particles inside polygon marking eac area
    for p=1:length(PA)
        [totalMass, totalN] = count_area(lons, lats, biomass, PA(p).PX, PA(p).PY);   
        subTT = timetable(dtime, totalN, totalMass);
        if nf == 1
            PA(p).TT = subTT;
        else
            PA(p).TT = [PA(p).TT; subTT];
        end
    end
end

%% plot it

for p=1:length(PA)

% figure
% plot(PA(p).TT.dtime, PA(p).TT.totalN)
% title(['Number of particles', PA(p).name])


hFig = figure;
maxVal = max(PA(p).TT.totalMass/10000);
set(hFig, 'Position', [10 10 700 300]); % fixed figure dimensions
set(hFig,'Color','w')
clf
c = [79.2/100 15.4/100 11.7/100]; % dot (jelly index) color
JellyIndDates = PS.dtime(PS.Ecoregion_code == PA(p).Ecoregion_code);
JellyInd2 = PS.index(PS.Ecoregion_code == PA(p).Ecoregion_code);
JellyInd2(JellyInd2==0) = NaN;
hold on
% plot model values
plot(PA(p).TT.dtime, PA(p).TT.totalMass/10000)
% plot observations
line([PA(p).TT.dtime(1) PA(p).TT.dtime(end)],[-maxVal/7 -maxVal/7],'Color','k')
scatter(JellyIndDates, JellyInd2*0-maxVal/7, JellyInd2*40,c,'filled')
% legend
xPos = datetime(year(PA(p).TT.dtime(end))-2,8,1);
xTextPos = datetime(year(PA(p).TT.dtime(end))-2,12,1);
yPos = maxVal * 0.95; 
scatter(xPos,yPos,40,c,'filled')
text(xTextPos,yPos,'1')
scatter(xPos,yPos*0.9,40*2,c,'filled')
text(xTextPos,yPos*0.9,'2')
scatter(xPos,yPos*0.8,40*3,c,'filled')
text(xTextPos,yPos*0.8,'3')
% end legend
hold off
xlim([PA(p).TT.dtime(1) PA(p).TT.dtime(end)])
ylim([-maxVal/7*1.5 maxVal*1.1])
yticks(linspace(0,yPos,5)) 
% title([PA(p).name,' total biomass'])

if saveIMG
    % save figure
    fOut = [IMGfold,PA(p).name,'_biomass_',strip(subFold,'/')];
    export_fig(fOut, '-png','-pdf');
end


% % LOGARITHMIC SCALE
% 
% hFig = figure;
% maxVal = max(PA(p).TT.totalMass/10000);
% set(hFig, 'Position', [10 10 700 300]); % fixed figure dimensions
% set(hFig,'Color','w')
% clf
% c = [79.2/100 15.4/100 11.7/100]; % dot (jelly index) color
% JellyIndDates = PS.dtime(PS.Ecoregion_code == PA(p).Ecoregion_code);
% JellyInd2 = PS.index(PS.Ecoregion_code == PA(p).Ecoregion_code);
% JellyInd2(JellyInd2==0) = NaN;
% hold on
% % plot model values
% plot(PA(p).TT.dtime, PA(p).TT.totalMass/10000)
% % plot observations
% line([PA(p).TT.dtime(1) PA(p).TT.dtime(end)],[-maxVal/7 -maxVal/7],'Color','k')
% scatter(JellyIndDates, JellyInd2*0-maxVal/7, JellyInd2*40,c,'filled')
% set(gca, 'YScale', 'log')
% % legend
% xPos = datetime(year(PA(p).TT.dtime(end))-2,8,1);
% xTextPos = datetime(year(PA(p).TT.dtime(end))-2,12,1);
% yPos = maxVal * 0.95; 
% scatter(xPos,yPos,40,c,'filled')
% text(xTextPos,yPos,'1')
% scatter(xPos,yPos*0.9,40*2,c,'filled')
% text(xTextPos,yPos*0.9,'2')
% scatter(xPos,yPos*0.8,40*3,c,'filled')
% text(xTextPos,yPos*0.8,'3')
% % end legend
% hold off
% xlim([PA(p).TT.dtime(1) PA(p).TT.dtime(end)])
% ylim([-maxVal/7*1.5 maxVal*1.1])
% yticks(linspace(0,yPos,5)) 
% % title([PA(p).name,' total biomass'])
% 
% if saveIMG
%     % save figure
%     fOut = [IMGfold,PA(p).name,'_biomass_',strip(subFold,'/'),'_log'];
%     export_fig(fOut, '-png','-pdf');
% end

end




%%
function [totalMass, totalN] = count_area(lons, lats, biomass, PX, PY)
% Number of particles and total biomass inside polygon
inArea = inpolygon(lons(:,:),lats(:,:),PX',PY');
mass = biomass.*inArea;
totalMass = sum(mass)'; % total biomass
totalN = sum(inArea)'; % number of particles in the simulation - doesn't change in a year
end
