%% Thesis: Patterns and trends of frontal activity in Australian marine hotspots
% Student: Kai Yang
% Student ID: 490063
% Data: MODIS daytime SST data (8-day average)
% Edge detection: Adaptive Canny Edge detection
%% read SST data and apply Adaptive Canny method
close all; clear;
datadir='F:\original_data\MODIS\MODIS_data\2003\';
filelist=dir([datadir,'*.nc']);
k=length(filelist);
for s=1:k
    filename=[datadir,filelist(s).name];
    ncid=netcdf.open(filename,'NC_NOWRITE');
    sst  = ncread(filename,'sst'); %뗍흙긴좆precip
    quality  = ncread(filename,'qual_sst'); %뗍흙긴좆time
    LonData  = ncread(filename,'lon'); %뗍흙긴좆lon
    LatData  = ncread(filename,'lat'); %뗍흙긴좆lat
    netcdf.close(ncid);   % 밑균匡숭
    % South East study region
    subAX=[146, 161 -43 -25];
    xidx=find(LonData>=subAX(1) & LonData<=subAX(2));
    yidx=find(LatData>=subAX(3) & LatData<=subAX(4));
    sub_sst=sst(xidx,yidx,:);
    % plot SST map and await to plot fronts on it
    figure(s);
    imagesc(LonData(xidx), LatData(yidx), sub_sst');
    set(gca, 'ydir', 'normal');
    xlabel('lon');
    ylabel('lat');
    h=colorbar;
    ylabel(h, '\circC');
    hold on
    % Applying Adaptive Canny Method
    [~,threshOut] = edge(sub_sst,'Canny');
    threshold = threshOut*2.0;
    BW1 = edge(sub_sst,'Canny',threshold);
    % plot fronts on SST maps
    BW1=double(BW1);
    idx2=find(BW1==0);
    BW1(idx2)=NaN;
    colormap(jet);
    h2 = pcolor(LonData(xidx),LatData(yidx), BW1'), shading flat;
    sub_lonx=LonData(xidx);
    sub_laty=LatData(yidx);
    save(['F:\original data\MODIS\2003\F_mat\A_F',num2str(s)], 'BW1','sub_lonx','sub_laty'); % save frontal maps
    save(['F:\original data\MODIS\2003\SST_mat\A_SST',num2str(s)], 'sub_sst','sub_lonx','sub_laty'); %save SST maps
end;
%% Further processing frontal maps

% Before processing:
% 1 = front
% NaN = no front detected or no front due to missing SST data

% After processing:
% 1 = front
% 0 = no front due to missing SST data
% 2 = no front detected
for i=1:46 % There are 46 frontal maps and SST maps in 2003.
    filepath='F:\original data\MODIS\2003\SST_mat\';
    % locate pixels with no SST data in SST maps
    load([filepath, 'A_SST', num2str(i), '.mat']);
    [r,c]=find(isnan(sub_sst));
    
    filepath='F:\original data\MODIS\2003\F_mat\';
    % assign values to these pixels in frontal maps
    load([filepath,'A_F',num2str(i),'.mat']);
    % no front due to missing SST data --> 0
    for m =1:length(r)
        BW1(r(m),c(m))= 0;
    end
    [row,col]=find(isnan(BW1));
    % no front detected --> 2
    for n=1:length(row)
        BW1(row(n),col(n))= 2;
    end
    save(['F:\original data\MODIS\1993\Ffill_mat\A1993_',num2str(i)], 'BW1','sub_lonx','sub_laty');
    % save processed frontal maps
end
%% Calculate Annual probability of frontal encounter(APFE)
filepath='F:\original data\MODIS\2003\Ffill_mat\'; % read processed frontal maps
% Calculate the number of times of a pixel was detected as a front for each
% pixel
FRONT=zeros(360,432);
for i=1:46 % There are totally 46 frontal maps in 2003.
    load([filepath,'A2003_',num2str(i),'.mat']);
    front{i}=zeros(360,432);
    for x=1:360;
        for y=1:432;
            if BW1(x,y)==1;
                front{i}(x,y)=1;
            end
        end
    end
    FRONT=FRONT+front{i};
end

FRONT=FRONT';
pfe=FRONT./(46-FRONT); % Calculate PFE for each pixel

% mask land area
[lon,lat] = meshgrid(sub_lonx,sub_laty);
land=landmask(lat,lon,'australia');
pfe(land)=NaN;
% plot PFE maps
imagesc(sub_lonx, sub_laty, flipud(rot90(pfe')));
set(gca,'YDir', 'normal');

% calculate annual PFE value for the whole study region
meanpfe = mean(pfe(:),'omitnan');

%% Calculate Monthly PFE (MPFE)
filepath='F:\original data\MODIS\West\2004\Ffill_mat\'; % read processed frontal maps
% Calculate the number of times of a pixel was detected as a front for each
% pixel
for r=1:12 % 12 months
    q=4; % There are 4 frontal maps in each month on average
    FRONT=zeros(360,432);
    m=1+(r-1)*4;
    M=m+3;
    for i=m:M
        load([filepath,'A2003_',num2str(i),'.mat']);
        front{i}=zeros(360,432);
        for x=1:360;
            for y=1:432;
                if BW1(x,y)==1;
                    front{i}(x,y)=1;
                end
            end
        end
        FRONT=FRONT+front{i};
    end
    
    FRONT=FRONT';
    pfe=FRONT./(q-FRONT); % calculate PFE value for each pixel
    
    % mask land area
    [lon,lat] = meshgrid(sub_lonx,sub_laty);
    land=landmask(lat,lon,'australia');
    pfe(land)=NaN;
    count=sum(land(:));% count the number of land pixels
    
    % calculate PFE value for the whole study region
    pfe=pfe(~isinf(pfe)) ;
    meanpfe(r)=mean(pfe(:),'omitnan');
end
%% Calculate frontal density (FD)
datadir='F:\original data\MODIS\2003\F_mat\'; % read processed frontal maps
filelist=dir([datadir,'A*.mat']); 

k=length(filelist);

for s=1:k
  filename=[filelist(s).name];
   load([datadir, filename]);
  e=find(BW1==1);
  Fpixelsum=length(e); 
  FD(s) = Fpixelsum./(360*432-count);
end
