%Buoy_Wind_scatter_CB.m
%
%Building a Data table that co-locates:
% - Buoy drift information from IABP
% - Winds (10m) from ERAi and ERA5
% - Geostrophic winds from ERAi 850hpa height
% - Ice Conditions: SIC (CDR-SIC); Thickness (CS2SMOS)
%
%CB - 01/03/2019
%   - 09/08/2019: U and V from geostrophic wind
%   - 11/11/2019: adding distance from the coast parameter
%   - 17/02/2020: rewriting to make a table format and adding ERA5 winds
%   - 20/06/2020: Finding PPFv4 Free drift estimates at buoy locations
%   - 08/07/2020: Adding PIOMAS thickness estimates
%% Common
%addpath('/aos/home/charles/Desktop/Matlab/') -> SIKU
addpath('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Code/Matlab/') % -> MACBOOK
FAITPaths_LOCAL_CB
Make_EASE_grid_CB

%IF needed (to-avoid re-computing everything when adding a new column):
load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/DATA_TABLE_DRIFT.mat');

%Path for buoy data on EASE:
% PPFcomponentEASEPath

%Era daily wind path:
ERAiDaily = [ERAipath 'Wind/Daily/EASE/'];
ERA5Daily = [ERA5path 'Wind/Daily/EASE/'];
ERAiGeo = [ERAipath 'GeostrophicWind/850hPaGeoWind/Daily/'];

%PPFv4_Wind_path
PPFv4Windp = '/Volumes/Yoda/FG2020/Data/NSIDCIceMotionV4/Raw/input_sources/wind/';

%PIOMAS daily path
PIOpath = '/Volumes/Yoda/FG2020/Data/PIOMAS/Daily/EASE/';

%Including a distance from the coast paramter, relevant for a free drift parameterization
% as described in Simizu et al. 2014
%Set up a two row matrix that lists the coordinates of all land points on
%the EASE grid:
load('bathyEASE');
[yland,xland]=find(bathyEASE>0);
LandCo = [xland,yland];

%% Variable names

VarNames = {
    'year';     ...
    'month';    ...
    'day';      ...
    'doy';      ...
    'EASE_x';   ...
    'EASE_y';   ...
    'u_buoy';   ...
    'v_buoy';   ...
    'id_buoy';  ...
    'u_ERAi';   ...
    'v_ERAi';   ...
    'theta_ERAi';...
    'ug_ERAi'; ...
    'vg_ERAi'; ...
    'thetag_ERAi';...
    'u_ERA5'; ...
    'v_ERA5'; ...
    'theta_ERA5';...
    'SIC'; ...
    'SIT'; ...
    'D2C'; ...
    'u_PPFv4Wind';...
    'v_PPFv4Wind';...
    'theta_PPFv4Wind';...
    'h_piomas'
    };

% When saving the table, make sure the variables are in the right order!!!

%% Co-locate data in space and time and make scatter plot

%One large matrix with daily buoy data for all years:
%[YEAR MONTH DAY DOY EASEx EASEy U V BUOYID]

% The BUOY_YYYY.mat files are produced by the routine
% ReadBuoyIceMotionPPFV4_CB.m in Matlab/PPF/

%Seed:
ALLBUOYMAT = NaN(1,9);

tic
for iyear = 1979:2019 % Cover satellite range

    load([PPFV4Proc 'BUOY_' int2str(iyear) '.mat']); %load buoy data

    ALLBUOYMAT = [ALLBUOYMAT; BUOY];
       
end
toc

ALLBUOYMAT(1,:)=[]; %-> Kill first row - just NaNs 

year = ALLBUOYMAT(:,1);
month = ALLBUOYMAT(:,2);
day = ALLBUOYMAT(:,3);
doy = ALLBUOYMAT(:,4);
EASE_x = ALLBUOYMAT(:,5);
EASE_y = ALLBUOYMAT(:,6);
u_buoy = ALLBUOYMAT(:,7);
v_buoy = ALLBUOYMAT(:,8);
id_buoy = ALLBUOYMAT(:,9);

week = ceil(doy./7);

%% Find Wind speed and sea ice concentration to ALLBUOYMAT

%Wind speed
u_ERAi = NaN(length(id_buoy),1);
v_ERAi = u_ERAi;
theta_ERAi = u_ERAi;
ug_ERAi = u_ERAi;
vg_ERAi = u_ERAi;
thetag_ERAi = u_ERAi;
u_ERA5 = u_ERAi;
v_ERA5 = u_ERAi;
theta_ERA5 = u_ERAi;
%Sea Ice Concentration
SIC = u_ERAi;
%Sea Ice Thickness
SIT = u_ERAi;
%Distance 2 Coast
D2C = u_ERAi;
%PPFv4 Free Drift
u_PPFv4Wind = u_ERAi;
v_PPFv4Wind = u_ERAi;
theta_PPFv4Wind = u_ERAi;
h_piomas = u_ERAi;

% LOAD DATA_TABLE and associate fields to these:

%%

tic
 f = waitbar(0,'Scanning IABP list');
for i = 1:length(id_buoy)
    
    %Identify time for each entry
    if month(i) < 10
       Mstring = ['0' int2str(month(i))];
    else
       Mstring = int2str(month(i));
    end
    
    if day(i) < 10
       Dstring = ['0' int2str(day(i))];
     else
       Dstring = int2str(day(i));
     end
   
    tstring = [ int2str(year(i)) Mstring Dstring];
    
    %find position of the buoy in the EASE grid:
    
    %PPFV4 : on 361x361 EASE Grid
     [r,c]=find(EASEx==round(EASE_x(i))...
            & EASEy==round(EASE_y(i)));
        
    if isempty(r)== 1
        continue %if the buoy don't have proper 
        %EASE grid coordinates, just skip and go to next datapoint    
    else
    
        erainame = [ERAiDaily 'ERAi_10m_UVwind_ease_' tstring '.mat'];
        if isfile(erainame)
        %load wind data for that day:
            load(erainame);
            %Get wind data at that position:       
            u_ERAi(i) = uwindEASE(r,c);
            v_ERAi(i) = vwindEASE(r,c);
        end
        
        eraigname = [ERAiGeo 'UV_Wind_Geostrophic_850hPa_' tstring '.mat'];
        if isfile(eraigname)
        %load geostrophic wind data for that day:
            load(eraigname);
            %Get wind data at that position:       
            ug_ERAi(i) = ug850(r,c);
            vg_ERAi(i) = vg850(r,c);
        end
        
        era5name = [ERA5Daily 'ERA5_10m_UV_wind_ease_' tstring '.mat'];
        if isfile(era5name)
        %load era5 wind data for that day:
            load(era5name);
            %Get wind data at that position:       
            u_ERA5(i) = uwindEASE(r,c);
            v_ERA5(i) = vwindEASE(r,c);
        end

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %load sea ice concentration data:
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        sicname = [CDRSICDholepath 'SSMI_CDR_SIC_' ...
                   tstring '.mat'];
        %Runs conditionally to the existence of data for that specific day
        %(Especially in the beginning of the time series when data is 
        %distributed every other day...)
        if isfile(sicname)
            load(sicname);
            SIC(i) = sicEASE(r,c);
        elseif ~isfile(sicname) % When data is only available every other day...
            bsic = [CDRSICDholepath 'SSMI_CDR_SIC_' datestr(datetime(year(i),month(i),day(i))-1,'yyyymmdd') '.mat'];% day before
            nsic = [CDRSICDholepath 'SSMI_CDR_SIC_' datestr(datetime(year(i),month(i),day(i))-1,'yyyymmdd') '.mat'];% day after
            if isfile(bsic) && isfile(nsic) % check both files exist
                losic = load(bsic);
                hisic = load(nsic);
                sicEASE = (losic.sicEASE + hisic.sicEASE)./2; %... take the average of day before and after 
                SIC(i) = sicEASE(r,c);    
            end
        end
        %Just making sure that there's no value above 100% SIC!
        if SIC(i)>1
            SIC(i)=1;
        end
        
        
        %Sea Ice Thickness
        %Go directly to CS2SMOS interpolated on the EASE grid (easier to associate
        %grid coordinates); They are distributed daily based on weekly averages. 
        cs2name = [CS2SMOShDPath 'cs2smos_EASE_' tstring '.mat'];
        if isfile(cs2name)==1
            load(cs2name);
            SIT(i) = cs2smos_h(r,c);
        end 
        
        %Distance2Coast
        %findind nearest coastal point and getting the distance:
        [~,dist] = dsearchn(LandCo,[EASE_x(i),EASE_y(i)]);
        D2C(i) = dist*25;
        
        %%%%%%%%%%%%%%%Angles%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Direction of the buoy drift [+90 degree of the Greenwich merridian]
        thetaB = atan2d(-v_buoy(i),u_buoy(i)); thetaB = thetaB + (thetaB<0)*360;% -> bring back to [-180 180] definition
        %Direction of the 10m wind ERAi
        thetaWerai = atan2d(-v_ERAi(i),u_ERAi(i)); thetaWerai = thetaWerai + (thetaWerai<0)*360;
        %Direction of the geostrophic wind ERAi
        thetaWgerai = atan2d(-vg_ERAi(i),ug_ERAi(i)); thetaWgerai = thetaWgerai + (thetaWgerai<0)*360;
        %Direction of the 10m wind ERA5
        thetaWera5 = atan2d(-v_ERA5(i),u_ERA5(i)); thetaWera5 = thetaWera5 + (thetaWera5<0)*360;
        
        dangleERAi = (thetaWerai-thetaB);
         if dangleERAi < -180
            dangleERAi = 360 + dangleERAi;
        elseif dangleERAi > 180
            dangleERAi = dangleERAi - 360;
         end
        
        dangleGEO = (thetaWgerai-thetaB);
        if dangleGEO < -180
            dangleGEO = 360 + dangleGEO;
        elseif dangleGEO > 180
            dangleGEO = dangleGEO - 360;
        end
        
        dangleERA5 = (thetaWera5-thetaB);
         if dangleERA5 < -180
         dangleERA5 = 360 + dangleERA5;
        elseif dangleERA5 > 180
            dangleERA5 = dangleERA5 - 360;
        end
        
        theta_ERAi(i) = dangleERAi;
        thetag_ERAi(i) = dangleGEO;
        theta_ERA5(i) = dangleERA5;
        
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end
    
    waitbar(i/length(id_buoy),f,['Scanning IABP List Part1: ' int2str(i) '/' int2str(length(id_buoy)) ', TOC=' int2str(toc) 's']);
    %pc = i/length(id_buoy)*100; 
    %disp([int2str(pc) '% done']);
    
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ADDING PPFv4 Free drift estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Path to data: PPFv4Windp

%PPFv4 Free Drift
%u_PPFv4Wind = u_ERAi;
%v_PPFv4Wind = u_ERAi;
%theta_PPFv4Wind = u_ERAi;

%DEaling with EASE50km to EASE25km conversion:
    ease25 = 1:361;
    [x25,y25] = meshgrid(ease25,ease25);
    ease50 = 1:181; 
    ease50to25= 2.*ease50-1;
    [x50,y50] = meshgrid(ease50to25,ease50to25);
    xin = x50(:);
    yin = y50(:);
    
    load ArcticOceanMask4IceEdge.mat

 tic
 f = waitbar(0,'Scanning IABP list -> PPF Wind dataset');
% Loop on years
for year_i = 1979:2019

    %Find & load the netcdf PPFv4Wind file for that year:
    filename = [PPFv4Windp 'im_from_wind_nh_50km_' int2str(year_i) '0101_' int2str(year_i) '1231_v4.1.nc'];
    u_ppffd = ncread(filename, 'u');
    v_ppffd = ncread(filename, 'v');
    
    % Loop on day of year
    
    for doy_i = 1:365
        
    % Find elements corresponding to year + day-of-year
        target = find(year==year_i & doy==doy_i);
        if ~isempty (target)
        
            %Find EASE coordinates of buoy locations for the given year/doy
                % Need to add +1 since we work on a 1:361grid system
            x_target = EASE_x(target)+1;
            y_target = EASE_y(target)+1;

            %Note the flip left-right + transpose operation to turn map (0,0)
            %corner to top-left
            u_work = fliplr(u_ppffd(:,:,doy_i))'; u_work = u_work(:);
            v_work = fliplr(v_ppffd(:,:,doy_i))'; v_work = v_work(:);
            nonan_ind = find(~isnan(u_work));
                % In the eventuality that a field is missing: eg.: 2015_337
                if isempty(nonan_ind)
                    u_work = fliplr((u_ppffd(:,:,doy_i-1)+u_ppffd(:,:,doy_i+1))./2)'; u_work = u_work(:);
                    v_work = fliplr((v_ppffd(:,:,doy_i-1)+v_ppffd(:,:,doy_i+1))./2)'; v_work = v_work(:);
                    nonan_ind = find(~isnan(u_work));
                end
                
            %Patch holes in the gridded data on 50km x 50km
            u50k = griddata(xin(nonan_ind),yin(nonan_ind),u_work(nonan_ind),x50,y50,'cubic');
            v50k = griddata(xin(nonan_ind),yin(nonan_ind),v_work(nonan_ind),x50,y50,'cubic');

            %Regrid to a 25kmx25km
            u25k = interp2(u50k);
            v25k = interp2(v50k);

            % Define xease
            % Define yease 

            % Find data @ buoy location
            uppf = interp2(x25,y25,u25k,x_target,y_target);
            vppf = interp2(x25,y25,v25k,x_target,y_target);

            %Store in PPFv4_free drift uv columns:
            u_PPFv4Wind(target,1) = uppf;
            v_PPFv4Wind(target,1) = -vppf; % Minus sign on the v to get it in the correct direction...

            waitbar(doy_i/365,f,['Scanning IABP List Part2: ' int2str(year_i) '/' int2str(doy_i) ', TOC=' int2str(toc) 's']);
        end
    end
end

%Calculate angle difference between Buoy and PPFv4FreeDrift
thetaB = atan2d(-v_buoy,u_buoy); thetaB = thetaB + (thetaB<0).*360;% -> bring back to [-180 180] definition
thetaPPFv4FD = atan2d(-v_PPFv4Wind,u_PPFv4Wind); thetaPPFv4FD = thetaPPFv4FD + (thetaPPFv4FD<0).*360;

danglePPFv4FD = (thetaPPFv4FD-thetaB);
less180 = find(danglePPFv4FD < -180);
danglePPFv4FD(less180) = 360 + danglePPFv4FD(less180);
more180 = find(danglePPFv4FD > 180);
danglePPFv4FD(more180) = danglePPFv4FD(more180) - 360;

theta_PPFv4Wind = danglePPFv4FD;


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ADDING PIOMAS thicknesses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
f = waitbar(0,'Scanning IABP list -> PIOMAS ice thickness');
% Loop on years
for year_i = 1979:2019

    %Find & load the EASE grid PIOMAS thickness mat-files:
    filename = [PIOpath 'PIOMAS_heff_EASEGrid_' int2str(year_i) '.mat'];
    load(filename);
    
    for doy_i = 1:365 % Loop on day of year
        
        h = heffEASE(:,:,doy_i); % pick the daily h field corresponding to doy
        
        target = find(year==year_i & doy==doy_i); % find whic h entries correspond to that day & year  
        
        if ~isempty (target)
            
            h_temp = NaN(numel(target),1); % temporary thickness storage
            
            for i_entry = 1:numel(target) % loop on corresponding buoy entries
            
                [r,c]=find(EASEx==round(EASE_x(target(i_entry)))...
                & EASEy==round(EASE_y(target(i_entry))));
        
                h_temp(i_entry) = h(r,c);
            end
            
            %Feed h_temp to h_piomas DATA_TABLE column:
            h_piomas(target) = h_temp;
            
         waitbar(doy_i/365,f,['Scanning IABP List Part3: ' int2str(year_i) '/' int2str(doy_i) ', TOC=' int2str(toc) 's']);
        end
    end
end

 %% Assemble data table and save output
 
%     If starting from previous DATA_DRIFT_TABLE:
            % load('/Users/charlesbrunette/Dropbox/FG2020/Analysis/Freedrift/Data/DATA_TABLE_DRIFT_save.mat');
            % year = DATA_TABLE_DRIFT.year;
            % month = DATA_TABLE_DRIFT.month;
            % day = DATA_TABLE_DRIFT.day;
            % doy = DATA_TABLE_DRIFT.doy;
            % EASE_x = DATA_TABLE_DRIFT.EASE_x;
            % EASE_y = DATA_TABLE_DRIFT.EASE_y;
            % u_buoy = DATA_TABLE_DRIFT.u_buoy;
            % v_buoy = DATA_TABLE_DRIFT.v_buoy;
            % id_buoy = DATA_TABLE_DRIFT.id_buoy;
            % u_ERAi = DATA_TABLE_DRIFT.u_ERAi;
            % v_ERAi = DATA_TABLE_DRIFT.v_ERAi;
            % theta_ERAi = DATA_TABLE_DRIFT.theta_ERAi;
            % ug_ERAi = DATA_TABLE_DRIFT.ug_ERAi;
            % vg_ERAi = DATA_TABLE_DRIFT.vg_ERAi;
            % thetag_ERAi = DATA_TABLE_DRIFT.thetag_ERAi;
            % u_ERA5 = DATA_TABLE_DRIFT.u_ERA5;
            % v_ERA5 = DATA_TABLE_DRIFT.v_ERA5;
            % theta_ERA5 = DATA_TABLE_DRIFT.theta_ERA5;
            % SIC = DATA_TABLE_DRIFT.SIC;
            % SIT = DATA_TABLE_DRIFT.SIT;
            % D2C = DATA_TABLE_DRIFT.D2C;
            % u_PPFv4Wind = DATA_TABLE_DRIFT.u_PPFv4Wind;
            % v_PPFv4Wind = DATA_TABLE_DRIFT.v_PPFv4Wind;
            % theta_PPFv4Wind = DATA_TABLE_DRIFT.theta_PPFv4Wind;
            % h_piomas = DATA_TABLE_DRIFT.h_piomas;

%%

DATA_TABLE_DRIFT = table( ...
    year,     ...
    month,    ...
    day,      ...
    doy,      ...
    EASE_x,   ...
    EASE_y,   ...
    u_buoy,   ...
    v_buoy,   ...
    id_buoy,  ...
    u_ERAi,   ...
    v_ERAi,   ...
    theta_ERAi,...
    ug_ERAi, ...
    vg_ERAi, ...
    thetag_ERAi, ...
    u_ERA5, ...
    v_ERA5, ...
    theta_ERA5,...
    SIC, ...
    SIT, ...
    D2C, ...
    u_PPFv4Wind, ...
    v_PPFv4Wind, ...
    theta_PPFv4Wind, ...
    h_piomas, ...
    'VariableNames',...
    VarNames);

save([FDpath_d 'DATA_TABLE_DRIFT_v2021.mat'], 'DATA_TABLE_DRIFT');

%Also provide a netcdf export! -> Need to work on it.
 
     
%% Calculate the drift angle relative to the wind speed
% 
% Appendix - Test for correct calculation of the angle direction
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Test with a 45 degree vector and a 30 degree vector
% %Should return a 15 degree difference
% 
% % u1 = 0; v1=-(-1); %Wind test
% % u2 = 1; v2=-(-0.5774); %buoy test
% % 
% % %theta1 = atand(v1/u1);
% % %theta2 = atand(v2/u2);
% % 
% % % matlab function atan2d returns angles in the range [-180,180] so we need
% % % to do a correction for negative angles to reteive angles on a [1 360]
% % % range..
% % 
% % theta1 = atan2d(v1,u1); theta1 = theta1 + (theta1<0)*360;
% % theta2 = atan2d(v2,u2); theta2 = theta2 + (theta2<0)*360;
% % 
% % %We need to take the inverse of the angle since v is defined downward
% % %positive and 0,0 is located at the top-left corner of the grid.
% % 
% % dang = (theta1-theta2);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%