%% Thesis: Patterns and trends of frontal activity in Australian marine hotspots
% Student: Kai Yang
% Student ID: 490063
% Data: AVHRR daytime SST data (6-day average)
% Edge detection: Adaptive Canny Edge detection
%% read SST data (quality level=4 & 5) and apply Adaptive Canny method
clear all; close all;
% read longitude and latitude
lon=ncread('**.nc', 'lon');
lat=ncread('**.nc', 'lat');
nx=length(lon);
ny=length(lat);
time=ncread('**.nc','time');
nt=length(time);
% read 6-day average SST data with quality level = 4 or 5
for r=6:6:nt
    data=ncread('**.nc', 'sea_surface_temperature', [1,1,r],[nx, ny, 1]);
    quality=ncread('**.nc', 'quality_level', [1,1,r],[nx, ny, 1]);
    for x=1:900;
        for y=1:1450;
            if quality(x,y)<4
                data(x,y)=nan;
            end
        end
    end
    % convert Kelvin into Celsius
    data=data-273.1;
    % South East Australia study region
    subAX=[146, 161 -43 -25];
    xidx=find(lon>=subAX(1) & lon<=subAX(2));
    yidx=find(lat>=subAX(3) & lat<=subAX(4));
    sub_data=data(xidx,yidx,:);
    % plot SST map and await to plot fronts on it
    figure(2);
    lonx=lon(xidx);
    laty=lat(yidx);
    imagesc(lonx, laty, sub_data');
    set(gca, 'ydir', 'normal');
    xlabel('lon');
    ylabel('lat');
    h=colorbar;
    ylabel(h, '\circC');
    hold on
    
    % Apply Adaptive Canny Method
    [~,threshOut] = edge(sub_data,'Canny');
    threshold = threshOut*2.0;
    BW1 = edge(sub_data,'Canny',threshold,2);
    % Plot detected front on SST map
    BW1=double(BW1);
    idx2=find(BW1==0);
    BW1(idx2)=NaN;
    colormap(jet);
    h2 = pcolor(lon(xidx),lat(yidx), BW1'), shading flat;
    title(filename);
    % save frontal maps and SST maps
    save(['F:\original data\EAST\1993\F_mat\A_F',num2str(r/6)],'lonx','laty', 'BW1'); % frontal maps
    save(['F:\original data\EAST\1993\SST_mat\A_SST',num2str(r/6)],'lonx','laty', 'sub_data'); % 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:58 % There are totally 58 SST maps and 58 frontal maps in 1993.
    % locate pixels with no SST data in SST maps
    filepath='F:\original data\EAST\1993\SST_mat\';
    load([filepath, 'A_SST', num2str(i), '.mat']);
    [r,c]=find(isnan(sub_data));
    % assign values to these pixels in frontal maps
    filepath='F:\original data\EAST\1993\F_mat\';
    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
    % no front detected --> 2
    [row,col]=find(isnan(BW1));
    for n=1:length(row)
        BW1(row(n),col(n))= 2;
    end
    % save
    save(['F:\original data\EAST\1993\Ffill_mat\A_1993_',num2str(i)], 'lonx','laty','BW1');
    % processed frontal maps
end
%% Calculate Annual probability of frontal encounter(APFE)
filepath='F:\original data\EAST\1993\Ffill_mat\'; % read processed frontal maps
% Calculate the number of times of a pixel was detected as a front for each
% pixel
FRONT=zeros(750,900);
for i=1:58 % There are totally 58 frontal maps in 1993.

    load([filepath,'A_1993_',num2str(i),'.mat']);
    front{i} = zeros(750, 900)
    for x=1:750;
        for y=1:900;
            if BW1(x,y)==1;
                front{i}(x,y)=1;
            end
        end
    end
    FRONT=FRONT+front{i};
end

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

% mask land area
[lon,lat] = meshgrid(lonx,laty);
land=landmask(lat,lon,'australia');
pfe(land)=NaN;

% plot PFE maps
imagesc(lonx,laty,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\EAST\1993\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=5; % There are 5 frontal maps in each month on average. 
    FRONT=zeros(750,900);
    m=1+(r-1)*5;
    M=m+4;
    if r==12
        M=M-2;
        q=3;
    else
        M=M;
        q=5;
    end
    
    for i=m:M
        
        load([filepath,'A_1993_',num2str(i),'.mat']);
        front{i}=zeros(750,900);
        for x=1:750;
            for y=1:900;
                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(lonx,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\EAST\1993\Ffill;_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./(750*900-count);
end
