%%
% set(0,'DefaultFigureWindowStyle','docked')

srcgeotif = 'guam_mosaic_large.tif';
srcbathy = 'd:\_data\shoals\guam_100_b_depths-East1.xyz.tif';
srctransect = 'guam_transect_coordinates.mat';

[geotif,cmap,R,bbox] = geotiffread(srcgeotif);
geotif = flipdim(geotif,1);
[bathy,cmap2,R2,bbox2] = geotiffread(srcbathy);
bathy = flipdim(bathy,1);
load(srctransect)

% take care of nans:
bathy(bathy == -32767) = nan;
% vertical up:
bathy = -bathy;

%You must set the surface's FaceColor to texturemap to use ZData and CData
%of different dimensions
x0 = 1E5*2.58721; y0 = 1E6*1.479417;% corresponds to first sensor
x = [-100 1000]; y = [-1200 800];% grid range
xlm = [-500 1000]; ylm = [-1200 800];% plot range

% generate tiff image coordinates
idim = size(geotif);
ixf = bbox(1):(bbox(2)-bbox(1))/(idim(2)-1):bbox(2);
iyf = [bbox(3):(bbox(4)-bbox(3))/(idim(1)-1):bbox(4)];

% generate bathy image coordinates
bdim = size(bathy);
bxf = bbox2(1):(bbox2(2)-bbox2(1))/(bdim(2)-1):bbox2(2);
byf = [bbox2(3):(bbox2(4)-bbox2(3))/(bdim(1)-1):bbox2(4)];

% UTM coordinates to meters: units are the same, origin is different
ixf = ixf - x0; iyf = iyf -y0;
bxf = bxf - x0; byf = byf -y0;
transect = [transect(:,1)-x0 transect(:,2)-y0];

% find geotif pixels corresponding to grid:
[tmp,xc(1)] = min(abs(ixf-x(1))); [tmp,xc(2)] = min(abs(ixf-x(2))); 
[tmp,yc(1)] = min(abs(iyf-y(1))); [tmp,yc(2)] = min(abs(iyf-y(2)));
% cut geotif:
geotif = geotif(yc(1):yc(2),xc(1):xc(2),:);
idim = size(geotif);

% find bathy pixels corresponding to grid:
[tmp,xc(1)] = min(abs(bxf-x(1))); [tmp,xc(2)] = min(abs(bxf-x(2))); 
[tmp,yc(1)] = min(abs(byf-y(1))); [tmp,yc(2)] = min(abs(byf-y(2)));
% cut bathy:
bathy = bathy(yc(1):yc(2),xc(1):xc(2),:);
bdim = size(bathy);

% plotting grid for the cut data:
xg = linspace(x(1),x(2),bdim(2));
yg = linspace(y(1),y(2),bdim(1));

h.f1 = figure('Renderer','opengl','Menubar','figure','Toolbar','figure');
h.a1 = axes('Nextplot','add','XLim',x,'YLim',y,'DataAspectRatio',[1 1 1/5],...
    'OuterPosition',[-0.032 .1 .5 .9]);

h.tr = plot(transect(:,1),transect(:,2),'.m');
% text(transect(:,1),transect(:,2)+5,transect_names,'Color','w')
xlabel('East (m)'), ylabel('North (m)')
set(h.tr,'MarkerSize',20)

h.t1 = text(700,-1180,'T. Hilmer,UH','FontSize',7,'FontWeight','bold','Color','w');
h.t2 = title('');

% grid minor

h.morph = surface(xg,yg,double(bathy),geotif,...
    'FaceColor','Texturemap','EdgeColor','none','Parent',h.a1); 

h.push1 = uicontrol('Style','pushbutton','Units','normalized','Position',[.02 .07 .04 .04],...
    'Callback',{@push,h},'String','CLONE');
h.push2 = uicontrol('Style','pushbutton','Units','normalized','Position',[.12 .07 .04 .04],...
    'Callback',{@push,h,1},'String','SAVE');

% ALTER THE LEFT AXES HOWEVER YOU LIKE.  IT'S 3D, SO TRY ROTATING.  THEN
% PUSH THE BUTTON.