% DragLakeCandidatesFigure_v1 

% Lea-Sophie Hoeyns, 06/02/2024

% This script makes a figure of the basal drag and basal drag
% coefficient map with possible lake candidates plotted on top. 

% Colormaps citation: Zhaoxu Liu / slandarer (2024). 
% 200 colormap (https://www.mathworks.com/matlabcentral/fileexchange/120088-200-colormap), 
% MATLAB Central File Exchange. Retrieved February 13, 2024.

% a. Basal drag map with lake candidates; b. Zoom-in box for Thwaites
% Glacier with lakes; c. Zoom-in box for Mercer and Whillens Ice Stream
% with lakes; d. Zoom-in box for MacAyeal Ice Stream with lakes

clearvars
close all
tic

%% DATA PATHS 
issmscriptpath='/albedo/work/projects/landice/ISSM/ISSM_installations/ISSM_v4p22/trunk/';
functionpath='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions';
slancm='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions/slanCM';
addpath(genpath(issmscriptpath),functionpath,slancm)
lakepath='~/1_Models/WAIS_FrictionOptimization_v2/Data/shapefiles/WAIS_shp/SubglacialLakes_WAIS/';
lakefiles_exp=dir('~/1_Models/WAIS_FrictionOptimization_v2/Data/shapefiles/WAIS_shp/SubglacialLakes_WAIS/*mask.exp'); 
lakefiles_shp=dir('~/1_Models/WAIS_FrictionOptimization_v2/Data/shapefiles/WAIS_shp/SubglacialLakes_WAIS/*.shp');
% File names:
lcurvefile='/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/OldLcurveRange/Run_2/WAIS.dragInv_HO_m3_Nocp_1DTemp_Lcurve_19-Jan-2024_Lambda_0.1.mat';
%lcurvefile='/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/OldLcurveRange/Run_1-1/WAIS.dragInv_HO_m3_Ncuas_1DTemp_LambdaBest-0.5_22-Jan-2024.mat'; % cell array of strings, include trailing ".mat", first one is reference
figname_best='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/BestDragLakeCanidatesFigure_v1_final.png';
figname='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/DragLakeCanidatesFigure_v1_final.png';

%% PARAMETERS 
% Model descriptions:
%modeldescriptext={'Best basal drag map, m=3, Budd, N_{CUAS}'};  % cell array of strings
modeldescriptext={'Basal drag map, m=3, Budd, N_{op}'};  % cell array of strings

% Map limits: 
%xlims=[-2000,600];                     % km
%ylims=[-1400,600];                     % km
% Thwaites Glacier zoom-in 
xlims_thw=[-1520,-1290];                % km
ylims_thw=[-520,-350];                  % km
% Mercer and Whillans Ice Stream zoom-in
xlims_pi=[-570 -210];                   % km
ylims_pi=[-640 -420];                   % km
% MacAyeal Ice Stream zoom-in
xlims_sip=[-850 -550];                  % km
ylims_sip=[-950 -670];                  % km
% Axes limits:
xtick=500;                              % km
ytick=500;                              % km
% Caxis limits:
draglims=[0,200];%150];%300];           % [1x2] kPa
dragtick=50;                            % kPa
flipdragcmap=0;                         % logical

% Aesthetic settings:
% Colors:
backplanecolor=[.7,.7,.7];              % valid color string or RGB
lcolor='w';                             % valid color string
% Colormaps:
dragcmap=(slanCM('twilight_s'));        % valid colormap string
% Fontsizes:
axesfontsize=16;                        % points
titlefontsize=16;                       % points
titlefontsize2=16;                      % points
colorbarfontsize=16;                    % points
legendfontsize=16;                      % points
markersize=36;                          % points
% Linewidths:
axeslinewidth=2;                        % points
axeslinewidth2=1;                       % points
llinewidth=1.5;                         % points
llinewidth2=2;                          % points

% Material parameters:
secondsperyear=60*60*24*365.25;

% Output size and resolution:
TileSpacing='compact';%'tight';
Padding='compact';
resolution=300;                         % dpi

%% PREPARATION  

% Load basic model:
load(lcurvefile,'md','MyInversionParameters')

% Unpack inversion parameters:
unpack(MyInversionParameters)

% Compute mesh properties:
% Extract 2D mesh:
try
    x2d=md.mesh.x2d;
    y2d=md.mesh.y2d;
    elements2d=md.mesh.elements2d;
    oldnumlayers=md.mesh.numberoflayers;
catch
    x2d=md.mesh.x;
    y2d=md.mesh.y;
    elements2d=md.mesh.elements;
    oldnumlayers=1;
end
numelements2d=size(elements2d,1);
numvertices2d=length(x2d);
% Compute edge lengths:
a=sqrt((x2d(elements2d(:,1))-x2d(elements2d(:,2))).^2+(y2d(elements2d(:,1))-y2d(elements2d(:,2))).^2);
b=sqrt((x2d(elements2d(:,2))-x2d(elements2d(:,3))).^2+(y2d(elements2d(:,2))-y2d(elements2d(:,3))).^2);
c=sqrt((x2d(elements2d(:,3))-x2d(elements2d(:,1))).^2+(y2d(elements2d(:,3))-y2d(elements2d(:,1))).^2);
% Compute area:
s=(a+b+c)/2;
Area=sqrt(s.*(s-a).*(s-b).*(s-c));
% Compute hydraulic diameter:
HydraulicDiameter_elements=2*Area./s;
% Compute element centers:
CentX=mean(x2d(elements2d),2);
CentY=mean(y2d(elements2d),2);
% Find mesh limits:
xlims=[min(x2d),max(x2d)];
ylims=[min(y2d),max(y2d)];
% Extract masks on 2D mesh and truncate:
OceanLevelSet_vertices=max(-1,min(1,md.mask.ocean_levelset(1:numvertices2d)));
IceLevelSet_vertices=max(-1,min(1,md.mask.ice_levelset(1:numvertices2d)));
% Compute ocean and ice area fractions:
OceanAreaFraction_elements=(1-mean(OceanLevelSet_vertices(elements2d),2))/2;
IceAreaFraction_elements=(1-mean(IceLevelSet_vertices(elements2d),2))/2;

% Generate the active mask:
ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0&md.mesh.vertexonboundary(1:numvertices2d)==0;
ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);

% Extract variables we need:
if numlayers>1 && ~strcmp(flowequation,'SSA')
    ModelVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1);
    ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
    %ObsVel=md.inversion.vel_obs(md.mesh.vertexonsurface==1);
    ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
    if ~strcmp(nsource,'none')
        EffectivePressure_vertices=md.friction.effective_pressure(md.mesh.vertexonbase==1);
    end
else
    ModelVel=md.results.StressbalanceSolution.Vel;
    ModelBasalVel=ModelVel;
    ObsVel=md.inversion.vel_obs;
    ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient;
    if ~strcmp(nsource,'none')
        EffectivePressure_vertices=md.friction.effective_pressure;
    end
end
 
% Compute drag:
Drag=(ModelDragCoeff.^2).*((ModelBasalVel/secondsperyear).^(1/m));
if ~strcmp(nsource,'none')
    Drag=Drag.*(EffectivePressure_vertices.^r);
end
% Set drag to zero under floating shelf and in ice-free areas:
Drag(OceanLevelSet_vertices<0|IceLevelSet_vertices>0)=0;

%% MAKE FIGURE Best Drag Map  
FigH = figure('Position', get(0, 'Screensize'));
tiledlayout(2,2,'TileSpacing',TileSpacing,'Padding',Padding); 

%% MAKE zoom-in SUBFIGURE Thwaites
nexttile(2,[1,1])
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
% Overlay possible lake candidates: 
hold on
for i=1:size(lakefiles_shp,1)
    S = shaperead([lakepath,lakefiles_shp(i).name]);%,'UseGeoCoords',true);
    hold on 
    [rows,columns]=size(struct2table(S));
    for k=1:rows
        hold on 
        sz=10;%8;%12;
        scatter(S(k).X/1000,S(k).Y/1000,markersize,'o','MarkerEdgeColor',lcolor,'LineWidth',llinewidth2)
    end
end 
hold on 
for j=1:size(lakefiles_exp,1)
     exp_orig=expread([lakepath,lakefiles_exp(j).name]);
     hold on
     l=plot(exp_orig.x/1000,exp_orig.y/1000,'-','Color',lcolor,'linewidth',llinewidth2);
end
% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
% Set up axes:
xlim(xlims_thw)
ylim(ylims_thw)
% xlabel('x (km)','FontSize',axesfontsize)
% ylabel('y (km)','FontSize',axesfontsize)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
% Make title:
title('\rm (b) Zoom-in Thwaites Glacier','FontSize',titlefontsize2)  

%% MAKE zoom-in SUBFIGURE Mercer and Whillans Ice Stream
nexttile(3,[1,1])
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
% Overlay possible lake candidates: 
hold on
for i=1:size(lakefiles_shp,1)
    S = shaperead([lakepath,lakefiles_shp(i).name]);%,'UseGeoCoords',true);
    hold on 
    [rows,columns]=size(struct2table(S));
    for k=1:rows
        hold on 
        sz=10;%8;%12;
        scatter(S(k).X/1000,S(k).Y/1000,markersize,'o','MarkerEdgeColor',lcolor,'LineWidth',llinewidth2)
    end
end 
hold on 
for j=1:size(lakefiles_exp,1)
     exp_orig=expread([lakepath,lakefiles_exp(j).name]);
     hold on
     l=plot(exp_orig.x/1000,exp_orig.y/1000,'-','Color',lcolor,'linewidth',llinewidth2);
end
% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
% Set up axes:
xlim(xlims_pi)
ylim(ylims_pi)
% xlabel('x (km)','FontSize',axesfontsize)
% ylabel('y (km)','FontSize',axesfontsize)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTick',[])
set(gca,'YTick',[])
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
box on 
% Make title:
title('\rm (c) Zoom-in Mercer and Whillans Ice Streams','FontSize',titlefontsize2)  

%% MAKE zoom-in SUBFIGURE MacAyeal Ice Stream 
nexttile(4,[1,1]) 
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
% Overlay possible lake candidates: 
hold on
for i=1:size(lakefiles_shp,1)
    S = shaperead([lakepath,lakefiles_shp(i).name]);%,'UseGeoCoords',true);
    hold on 
    [rows,columns]=size(struct2table(S));
    for k=1:rows
        hold on 
        sz=10;%8;%12;
        scatter(S(k).X/1000,S(k).Y/1000,markersize,'o','MarkerEdgeColor',lcolor,'LineWidth',llinewidth2)
    end
end 
hold on 
for j=1:size(lakefiles_exp,1)
     exp_orig=expread([lakepath,lakefiles_exp(j).name]);
     hold on
     l=plot(exp_orig.x/1000,exp_orig.y/1000,'-','Color',lcolor,'linewidth',llinewidth2);
end
% Set up color scale and colorbar:
colormap(gca,dragcmap) 
caxis(draglims)
% Set up axes:
xlim(xlims_sip)
ylim(ylims_sip)
% xlabel('x (km)','FontSize',axesfontsize)
% ylabel('y (km)','FontSize',axesfontsize)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
% Make title:
title('\rm (d) Zoom-in MacAyeal Ice Stream','FontSize',titlefontsize2)

%% PLOT best drag map SUBFIGURE 
nexttile
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
hold on 
rectangle('Position',[xlims_thw(1) ylims_thw(1) xlims_thw(2)-xlims_thw(1) ylims_thw(2)-ylims_thw(1)],'LineWidth',axeslinewidth)
text(-1410,-300,['(',char(96+2),')'],'color','w','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_pi(1) ylims_pi(1) xlims_pi(2)-xlims_pi(1) ylims_pi(2)-ylims_pi(1)],'LineWidth',axeslinewidth)
text(-540,-580,['(',char(96+3),')'],'color','w','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_sip(1) ylims_sip(1) xlims_sip(2)-xlims_sip(1) ylims_sip(2)-ylims_sip(1)],'LineWidth',axeslinewidth)
text(-710,-750,['(',char(96+4),')'],'color','w','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')

% Overlay possible lake candidates: 
hold on
for i=1:size(lakefiles_shp,1)
    S = shaperead([lakepath,lakefiles_shp(i).name]);%,'UseGeoCoords',true);
    hold on 
    [rows,columns]=size(struct2table(S));
    for k=1:rows
        hold on 
        sz=12;%8;%12;
        scatter(S(k).X/1000,S(k).Y/1000,sz,'o','MarkerEdgeColor',lcolor,'LineWidth',llinewidth)
    end
end 
hold on 
for j=1:size(lakefiles_exp,1)
     exp_orig=expread([lakepath,lakefiles_exp(j).name]);
     hold on
     l=plot(exp_orig.x/1000,exp_orig.y/1000,'-','Color',lcolor,'linewidth',llinewidth);
end

% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
hc=colorbar;
set(get(hc,'Ylabel'),'String','kPa','FontSize',colorbarfontsize)
%set(hc,'YTick',dragtick*[ceil(draglims(1)/dragtick):1:floor(draglims(2)/dragtick)])
set(hc,'Ticks',[0 50 100 150 200],'TickLabels',{'0','50','100','150','\geq 200'})
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize,'LineWidth',axeslinewidth)
% Set up axes:
xlim(xlims/1000)
ylim(ylims/1000)
set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
xlabel('x (km)','FontSize',axesfontsize)
ylabel('y (km)','FontSize',axesfontsize)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'FontSize',axesfontsize,'LineWidth',axeslinewidth)
box on 
% Make legend: 
legend([l],{'Subglacial lakes'});
legend('Location','northeast','FontSize',legendfontsize,'Color',backplanecolor); 
% Make title:
title(['\rm',modeldescriptext],'FontSize',titlefontsize)

%% SAVE FIGURE
saveas(FigH,figname,'png');
exportgraphics(FigH,figname,'resolution',resolution)
%saveas(FigH,figname_best,'png');
%exportgraphics(FigH,figname_best,'resolution',resolution)

disp('Done!')
toc