% InversionSettingFigure_v1

% Mike Wolovick, 12/13/2022

% This script makes the setting figure for the inversion paper.  

% Subplots:  
% a) Surfelev (colored and hillshaded)
% b) Bedelev (colored and hillshaded)
% c) Driving stress (in flow direction)
% d) Velocity observations

clear vars
tic
close all

%% Parameters:

% File names:
gridfile='/work/ollie/mwolovic/Inputs/AntarcticaCombinedGrids_v4.mat';
velfile='/work/ollie/mwolovic/Inputs/AntarcticaMergedVel_MEASURESv2withphase_plus_FRmosaic_velonly_900m.mat';
velsuffix='_merged';  
boundaryfile='/work/ollie/mwolovic/FRISModeling/Boundaries/BoundaryPickingTest_v3_50km_expanded.exp';
cmapfiles={'/home/ollie/mwolovic/GMA_default.mat';...
    '/home/ollie/mwolovic/GMA_combined.mat';...
    '/home/ollie/mwolovic/GMA_ocean.mat';...
    '/home/ollie/mwolovic/GMA_land.mat'};  % cell array of strings (can be commented)
figname='/home/ollie/mwolovic/InversionSettingFigure_v2.png'; 

% Ticks and limits and colorscales:
xlims=[-1750,1000];               % [1x2] km
xtick=500;                               % km
ylims=[-250,1750];                % [1x2] km
ytick=500;                               % km
surflims=[0,4100];                    % [1x2] m
surftick=1000;                           % m
surfcmap='GMA_land';            % valid colormap string, potentially including names of cmapfiles
bedlims=[-2000,2000];            % [1x2] m
bedtick=500;                           % m
bedcmap='GMA_ocean';          % valid colormap string, potentially including names of cmapfiles
stresslims=[-300,300];              % [1x2] kPa
stresstick=100;                             % kPa
stresscmap='FadedJet';             % valid colormap string, potentially including names of cmapfiles
vellims=[0.5,1500];                   % 1x2 m/yr (will be shown on log scale)
velcmap='GMA_default';            % valid colormap string, potentially including names of cmapfiles

% Lighting settings:
lightangle_surf=30;                    % deg up from horizon
lightangle_bed=60;                    % deg up from horizon
diffuselight=.5;                          % unitless
ambientlight=.5;                         % unitless
specularlight=0;                         % unitless

% Lat/lon graticule:
latinterval=5;                              % degrees
loninterval=45;                           % degrees
latlims=[-90,-60];                        % [1x2] degrees
lonlims=[-180,180];                    % [1x2] degrees
numgratpoints=1000;                % integer
gratcolor='k';                             % valid color string or RGB
gratlinewidth=0.5;                      % points

% Inset controls:
insetxlims=[-2700,2700];           % [1x2] km
insetylims=[-2700,2700];            % [1x2] km
insetxrange=700;                        % km
insetmaskcolor='k';                            % valid color string or RGB
insetmasklinewidth=1;                     % points
insetoutlinewidth=2;                      % points
insetinlinewidth=1;                         % points

% Other aesthetic settings:
boundarylinewidth=2;                % points
boundarylinecolor='r';                % valid color string or RGB
masklinewidth=1;                       % points
maskcolor='w';                            % valid color string or RGB
axesfontsize=14;                       % points
titlefontsize=15;                          % points
colorbarfontsize=14;                   % points

% Subplot controls:
vertaxesbuffer=.07;             % unitless
titlebuffer=.05;                     % unitless
horzaxesbuffer=.04;             % unitless
horznotextbuffer=.00;          % unitless
horzcolorbarbuffer=.11;        % unitless

% Output size and resolution:
pagewidth=18;                       % inches (height computed automatically)
resolution=300;                      % dpi

% Physical parameters:
rho_i=917;                               % kg/m^3
g=9.8;                                      % m/s^2

%% Preparation:

% Hardcoded subplot arrangement:
numvertplots=2;
numhorzplots=2;
% Compute subplot dimensions:
% Compute available space:
availablehorzspace=1-(numhorzplots-1)*horznotextbuffer-horzaxesbuffer-numhorzplots*horzcolorbarbuffer;
availablevertspace=1-numvertplots*titlebuffer-vertaxesbuffer;
% Compute boxes:
Boxes=cell(numvertplots,numhorzplots);
for ii=1:numvertplots
    for jj=1:numhorzplots
        Boxes{numvertplots-ii+1,jj}=[horzaxesbuffer+(jj-1)*(horznotextbuffer+horzcolorbarbuffer+availablehorzspace/numhorzplots),...
            vertaxesbuffer+(ii-1)*(titlebuffer+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end
% Modify boxes for colorbars:
for ii=1:numvertplots*numhorzplots
    Boxes{ii}=Boxes{ii}+[0,0,horzcolorbarbuffer,0];
end

% Compute page height:
pageheight=pagewidth*((ylims(2)-ylims(1))/(xlims(2)-xlims(1)))*(availablehorzspace/availablevertspace)*(numvertplots/numhorzplots);

% Compute inset box placement:
insetzoomratio=insetxrange/(insetxlims(2)-insetxlims(1));
insetbox=[xlims(1),ylims(2)-insetzoomratio*(insetylims(2)-insetylims(1)),xlims(1)+insetzoomratio*(insetxlims(2)-insetxlims(1)),ylims(2)];  % [left, bottom, right, top]

% Load colormaps:
if exist('cmapfiles','var')
    for ii=1:length(cmapfiles)
        % load the colormap:
        load(cmapfiles{ii},'thiscolormap')
        % rename the colormap:
        slashnum=strfind(cmapfiles{ii},'/');
        thiscmapname=cmapfiles{ii}(slashnum(end)+1:end-4);
        eval([thiscmapname,'=thiscolormap;'])
    end
end

% Load domain boundary:
DomainBoundary=expread(boundaryfile);

% Load grids:
load(gridfile,'X','Y','BedElev','SurfElev','IceThick','Mask','stereoprojection')

% Identify grid bounds within the figure:
ind1_grid=find(Y>=ylims(1)*1e3,1,'first');
ind2_grid=find(Y<=ylims(2)*1e3,1,'last');
ind3_grid=find(X>=xlims(1)*1e3,1,'first');
ind4_grid=find(X<=xlims(2)*1e3,1,'last');

% Rename coordinate variables:
X_grid=X;
Y_grid=Y;

% Reorder the mask:
ReorderedMask=Mask;  % 0 (ocean) and 2 (grounded ice) remain the same
ReorderedMask(Mask==3)=1; % shelf
ReorderedMask(Mask==1)=3; % rock
Mask=ReorderedMask;
clear ReorderedMask

% Evaluate grid size:
[ysize,xsize]=size(Mask);
dx=X_grid(2)-X_grid(1); % m

% Load velocity:
load(velfile,'X','Y',['U',velsuffix],['V',velsuffix])

% Rename variables:
X_vel=X;
Y_vel=Y;
eval(['U_vel=U',velsuffix,';'])
eval(['V_vel=V',velsuffix,';'])
eval(['clear X Y U',velsuffix,' V',velsuffix]) 

% Identify bounds on velocity grid:
ind1_vel=find(Y_vel>=ylims(1)*1e3-2*dx,1,'first');
ind2_vel=find(Y_vel<=ylims(2)*1e3+2*dx,1,'last');
ind3_vel=find(X_vel>=xlims(1)*1e3-2*dx,1,'first');
ind4_vel=find(X_vel<=xlims(2)*1e3+2*dx,1,'last');

% Interpolate velocity to main grid:
U=interp2(X_vel(ind3_vel:ind4_vel),Y_vel(ind1_vel:ind2_vel),U_vel(ind1_vel:ind2_vel,ind3_vel:ind4_vel),repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]),'linear');
V=interp2(X_vel(ind3_vel:ind4_vel),Y_vel(ind1_vel:ind2_vel),V_vel(ind1_vel:ind2_vel,ind3_vel:ind4_vel),repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]),'linear');

% Compute velopcity magnitude:
VelMag=sqrt(U.^2+V.^2);

% Compute driving stress:
[SurfGradX,SurfGradY]=gradient(SurfElev,dx);
DrivingStressX=-rho_i*g*IceThick.*SurfGradX;
DrivingStressY=-rho_i*g*IceThick.*SurfGradY;

% Resolve driving stress in the flow direction:
DrivingStressFlowDir=(DrivingStressX.*U+DrivingStressY.*V)./VelMag;
DrivingStressFlowDir(Mask==0|Mask==3|VelMag==0|isnan(VelMag))=0;

% Compute graticule:
[X_latpoints,Y_latpoints]=projfwd(stereoprojection,repmat(latinterval*(ceil(latlims(1)/latinterval):floor(latlims(2)/latinterval)),[numgratpoints,1]),...
    repmat(linspace(lonlims(1),lonlims(2),numgratpoints)',[1,length(ceil(latlims(1)/latinterval):floor(latlims(2)/latinterval))]));
[X_lonpoints,Y_lonpoints]=projfwd(stereoprojection,repmat(linspace(latlims(1),latlims(2),numgratpoints)',[1,length(ceil(lonlims(1)/loninterval):floor(lonlims(2)/loninterval))]),...
    repmat(loninterval*(ceil(lonlims(1)/loninterval):floor(lonlims(2)/loninterval)),[numgratpoints,1]));

% Generate contour structure for mask:
C_main=contourc(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid),[.5,1.5,2.5]);
C_inset=contourc(insetbox(1)+insetzoomratio*(X_grid/1000-insetxlims(1)),insetbox(2)+insetzoomratio*(Y_grid/1000-insetylims(1)),Mask,[.5,1.5]);

%% Make Figure:

% Make a figure:
figure(1)

% Make faded jet:
MakeFadedJet;


% Surface plot:
% Make subplot:
h1=subplot('Position',Boxes{1,1});
% Generate RGB:
% pre-allocate:
RGB=zeros(ind2_grid-ind1_grid+1,ind4_grid-ind3_grid+1,3);
% get this colormap:
eval(['thiscmap=colormap(',surfcmap,');'])
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(SurfElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid)-surflims(1))/(surflims(2)-surflims(1)))),'linear');
end
% set color outside of the ice:
RGB(repmat(Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==0,[1,1,3]))=1;
% Draw surface:
surf(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,SurfElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid),RGB,...
    'EdgeColor','none','FaceColor','interp','FaceLighting','gouraud','SpecularStrength',specularlight,'DiffuseStrength',diffuselight,'AmbientStrength',ambientlight)
view([0,90])
hold on
% Compute max surface:
topz=max(max(SurfElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid)));
% Generate light objects:
light('Position',[cosd(lightangle_surf),0,sind(lightangle_surf)])
light('Position',[0,cosd(lightangle_surf),sind(lightangle_surf)])
% Plot graticule:
plot3(X_latpoints/1000,Y_latpoints/1000,(topz+1)*ones(size(X_latpoints)),'Color',gratcolor,'LineWidth',gratlinewidth)
plot3(X_lonpoints/1000,Y_lonpoints/1000,(topz+1)*ones(size(X_lonpoints)),'Color',gratcolor,'LineWidth',gratlinewidth)
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot3(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),(topz+2)*ones(size(C_main(1,startind+1:stopind))),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot3(DomainBoundary.x/1000,DomainBoundary.y/1000,(topz+3)*ones(size(DomainBoundary.x)),'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Make inset:
% make back plane:
patch('Vertices',[[insetbox(1);insetbox(3);insetbox(3);insetbox(1);insetbox(1)],...
    [insetbox(2);insetbox(2);insetbox(4);insetbox(4);insetbox(2)],(topz+4)*ones(5,1)],'Faces',1:5,'EdgeColor','k','FaceColor','w','LineWidth',insetoutlinewidth,'FaceLighting','flat','AmbientStrength',1)
% Plot graticule:
% for ii=1:size(X_latpoints,2)
%     ThesePoints=X_latpoints(:,ii)/1000>=insetxlims(1)&X_latpoints(:,ii)/1000<=insetxlims(2)&Y_latpoints(:,ii)/1000>=insetylims(1)&Y_latpoints(:,ii)/1000<=insetylims(2);
%     plot3(insetbox(1)+insetzoomratio*(X_latpoints(ThesePoints,ii)/1000-insetxlims(1)),insetbox(2)+insetzoomratio*(Y_latpoints(ThesePoints,ii)/1000-insetylims(1)),(topz+5)*ones(sum(ThesePoints),1),'Color',gratcolor,'LineWidth',gratlinewidth)
% end
% for ii=1:size(X_lonpoints,2)
%     ThesePoints=X_lonpoints(:,ii)/1000>=insetxlims(1)&X_lonpoints(:,ii)/1000<=insetxlims(2)&Y_lonpoints(:,ii)/1000>=insetylims(1)&Y_lonpoints(:,ii)/1000<=insetylims(2);
%     plot3(insetbox(1)+insetzoomratio*(X_lonpoints(ThesePoints,ii)/1000-insetxlims(1)),insetbox(2)+insetzoomratio*(Y_lonpoints(ThesePoints,ii)/1000-insetylims(1)),(topz+5)*ones(sum(ThesePoints),1),'Color',gratcolor,'LineWidth',gratlinewidth)
% end
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_inset(2,startind);
    plot3(C_inset(1,startind+1:stopind),C_inset(2,startind+1:stopind),(topz+6)*ones(size(C_inset(1,startind+1:stopind))),'Color',insetmaskcolor,'LineWidth',insetmasklinewidth)
    startind=stopind+1;
    if startind>size(C_inset,2)
        done=1;
    end
end
% Plot box representing main figure:
plot3(insetbox(1)+insetzoomratio*([xlims(1),xlims(2),xlims(2),xlims(1),xlims(1)]-insetxlims(1)),...
    insetbox(2)+insetzoomratio*([ylims(1),ylims(1),ylims(2),ylims(2),ylims(1)]-insetylims(1)),(topz+7)*ones(1,5),'k','LineWidth',insetinlinewidth)
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
ylabel('Northing (km)','FontSize',axesfontsize)
set(gca,'XTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale and colorbar:
colormap(gca,thiscmap)
caxis(surflims)
hc=colorbar;
set(hc,'YTick',surftick*[ceil(surflims(1)/surftick):1:floor(surflims(2)/surftick)])
set(get(hc,'Ylabel'),'String','m','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('a) Surface elevation','FontSize',titlefontsize)


% Bed plot:
% Make subplot:
h2=subplot('Position',Boxes{1,2});
% Generate RGB:
% pre-allocate:
RGB=zeros(ind2_grid-ind1_grid+1,ind4_grid-ind3_grid+1,3);
% get this colormap:
eval(['thiscmap=colormap(',bedcmap,');'])
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(BedElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid)-bedlims(1))/(bedlims(2)-bedlims(1)))),'linear');
end
% Draw surface:
surf(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,BedElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid),RGB,...
    'EdgeColor','none','FaceColor','interp','FaceLighting','gouraud','SpecularStrength',specularlight,'DiffuseStrength',diffuselight,'AmbientStrength',ambientlight)
view([0,90])
hold on
% Compute max surface:
topz=max(max(BedElev(ind1_grid:ind2_grid,ind3_grid:ind4_grid)));
% Generate light objects:
light('Position',[cosd(lightangle_bed),0,sind(lightangle_bed)])
light('Position',[0,cosd(lightangle_bed),sind(lightangle_bed)])
% Plot graticule:
plot3(X_latpoints/1000,Y_latpoints/1000,(topz+1)*ones(size(X_latpoints)),'Color',gratcolor,'LineWidth',gratlinewidth)
plot3(X_lonpoints/1000,Y_lonpoints/1000,(topz+1)*ones(size(X_lonpoints)),'Color',gratcolor,'LineWidth',gratlinewidth)
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot3(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),(topz+2)*ones(size(C_main(1,startind+1:stopind))),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot3(DomainBoundary.x/1000,DomainBoundary.y/1000,(topz+3)*ones(size(DomainBoundary.x)),'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
set(gca,'YTickLabel',[])
set(gca,'XTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale and colorbar:
colormap(gca,thiscmap)
caxis(bedlims)
hc=colorbar;
set(hc,'YTick',bedtick*[ceil(bedlims(1)/bedtick):1:floor(bedlims(2)/bedtick)])
set(get(hc,'Ylabel'),'String','m','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('b) Bed elevation','FontSize',titlefontsize)


% Velocity plot:
% Make subplot:
h3=subplot('Position',Boxes{2,1});
% Generate RGB:
% pre-allocate:
RGB=zeros(ind2_grid-ind1_grid+1,ind4_grid-ind3_grid+1,3);
% get this colormap:
eval(['thiscmap=colormap(',velcmap,');'])
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(log(VelMag(ind1_grid:ind2_grid,ind3_grid:ind4_grid))-log(vellims(1)))/(log(vellims(2))-log(vellims(1))))),'linear');
end
% set areas outside the ice or no data to gray:
RGB(repmat(Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==0|Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==3|...
    isnan(VelMag(ind1_grid:ind2_grid,ind3_grid:ind4_grid))|VelMag(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==0,[1,1,3]))=.75;
% Draw image:
image(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,RGB)
hold on
% Plot graticule:
plot(X_latpoints/1000,Y_latpoints/1000,'Color',gratcolor,'LineWidth',gratlinewidth)
plot(X_lonpoints/1000,Y_lonpoints/1000,'Color',gratcolor,'LineWidth',gratlinewidth)
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
set(gca,'ydir','normal')
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
ylabel('Northing (km)','FontSize',axesfontsize)
xlabel('Easting (km)','FontSize',axesfontsize)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale:
colormap(gca,thiscmap)
caxis(log10(vellims))
% Compute first base and tick for color bar:
thisbase=10^floor(log10(vellims(1)));
thistick=thisbase*ceil(vellims(1)/thisbase);
% Create all ticks and labels for color bar:
velticks=zeros(100,1);
velticklabels=cell(100,1);
for ii=1:100
    % Record this tick:
    velticks(ii)=log10(thistick);
    % Check if we hit the next power of ten:
    if abs(thistick-10*thisbase)/thisbase<1e-12 || (ii==1 && abs(thistick-thisbase)/thisbase<1e-12)
        % Generate label:
        velticklabels{ii}=['10^{',num2str(round(log10(thistick))),'}'];
        % Advance base:
        thisbase=thistick;
    end
    % Advance tick:
    thistick=thistick+thisbase;
    % Break from loop and truncate:
    if thistick>(1+1e-12)*vellims(2)
        velticks=velticks(1:ii);
        velticklabels=velticklabels(1:ii);
        break
    end
end
% Make colorbar:
hc=colorbar;
set(hc,'YTick',velticks)
set(hc,'YTickLabel',velticklabels)
set(get(hc,'Ylabel'),'String','m/a','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('c) Observed velocity','FontSize',titlefontsize)


% Driving stress plot:
% Make subplot:
h4=subplot('Position',Boxes{2,2});
% Generate RGB:
% pre-allocate:
RGB=zeros(ind2_grid-ind1_grid+1,ind4_grid-ind3_grid+1,3);
% get this colormap:
eval(['thiscmap=colormap(',stresscmap,');'])
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(DrivingStressFlowDir(ind1_grid:ind2_grid,ind3_grid:ind4_grid)/1e3-stresslims(1))/(stresslims(2)-stresslims(1)))),'linear');
end
% set areas outside the ice or no data to gray:
RGB(repmat(Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==0|Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==3|...
    isnan(VelMag(ind1_grid:ind2_grid,ind3_grid:ind4_grid))|VelMag(ind1_grid:ind2_grid,ind3_grid:ind4_grid)==0,[1,1,3]))=.75;
% Draw image:
image(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,RGB)
hold on
% Plot graticule:
plot(X_latpoints/1000,Y_latpoints/1000,'Color',gratcolor,'LineWidth',gratlinewidth)
plot(X_lonpoints/1000,Y_lonpoints/1000,'Color',gratcolor,'LineWidth',gratlinewidth)
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
set(gca,'ydir','normal')
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
set(gca,'YTickLabel',[])
xlabel('Easting (km)','FontSize',axesfontsize)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale:
colormap(gca,thiscmap)
caxis(stresslims)
% Make colorbar:
hc=colorbar;
set(hc,'YTick',stresstick*[ceil(stresslims(1)/stresstick):1:floor(stresslims(2)/stresstick)])
set(get(hc,'Ylabel'),'String','kPa','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('d) Driving stress in the flow direction','FontSize',titlefontsize)


% Reset previous axes colormaps:
eval(['colormap(h1,',surfcmap,')'])
eval(['colormap(h2,',bedcmap,')'])
eval(['colormap(h3,',velcmap,')'])

% Save figure:
set(gcf,'PaperSize',[pagewidth,pageheight])
set(gcf,'PaperPosition',[0,0,pagewidth,pageheight])
print('-dpng',figname,['-r',num2str(resolution)])

% Final display:
disp('Done!')
toc