%
% PURPOSE:  Cut a large .grd map up into equally-sized tiles.  Each tile
% is also a .grd.  Then use grdOctreeMaker to convert the .grd tiles to
% .bo's. 
%
% First, overlay Dave's EMARK_vents map with grid
% Rob M., 12 July 2018
% Rob M., 29 Mar 2023
%
% To plot the data:
%
% cd /Users/rob/logs/2023.083.00
% close all;clear all; cd ../2023.083.00;missionFileName='EMARK_13v3.cfg.pp';...
% loc='EMARK';loadData;deci=20;vertEx=1;loadGrd;review;plotSurf;set(gca,'zlim',[-4532 -3000])
%
% To plot the sim:
%
% cd /Volumes/qnx_partition/midAtlantic/dorado-auv-qnx/auv/altex/onboard/logs
% close all;clear all;cd ../2023.094.08;deci=5;loc='portSlope';vertEx=1;...
% missionFileName='PortTrapezoid.cfg.pp';loadData;loadGrd;review;plotSurf
%
% Path name to map:
  parentMapName = '~/Maps/EMARK_vents_2023/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd';
%
%  Use: 
%
%    gmt grdinfo EMARK_vents_2023/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd
%
%  to find the boundaries of the .grd.
%
% Now define the tile boundaries.  Let them be 1.5kmx1.5km, with center
% separations of 1 km, which creates non-overlapped regions of
% 1x1km. So there are 6x6/1.0^2 = 36 tiles.
%
% Start with arrays of the lower left corners of the NON-OVERLAPPED regions:
%
% For ~/Maps/EMARK_vents_2023/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd
  tileSize = 1000;
  edgeSize = 250;
  ec = [ 499000:tileSize: 504000]; % Easting corners
  nc = [2593000:tileSize:2599000]; % Northing corners

if(~exist('position'))
  position = get(0,'screensize');
  npposition(1:2) = position(1:2) + round([position(3)*.02 position(4)*.04]);
  npposition(3:4) = round([position(3)*.95 position(4)*.90]);
  position = npposition;
end

  lstr={};
  if exist( 'plgFigPtr' ) ~= 1
    plgFigPtr= figure;
    set( plgFigPtr, 'Name', 'Axial Grid', 'DeleteFcn', 'clear plgFigPtr' );
    set(gcf,'position',position)
    if(exist('fontsize')) 
      MbTrnFontsize=fontsize;
    else
      MbTrnFontsize=20;
    end
    MbTrnFontweight='normal';
    set(gca,'fontsize',MbTrnFontsize,'fontweigh',MbTrnFontweight)
    set(gcf,'defaultaxesfontsize',MbTrnFontsize,'defaulttextfontsize', ...
	    MbTrnFontsize)
  else
    figure(plgFigPtr);
  end

% parentMapName  
  if( exist('plg') == 0 )
    if( exist(parentMapName) == 2 )
      [plg.eastings plg.northings plg.z]=grdread2rob(parentMapName);
    else
      fprintf('Map %s doesn''t exist.  Exit.\n', parentMapName );
    end
  end
  
%
% Compute the tile boundaries using ec and nc defined above.  
  llc = [];
  for( i=1:length(nc) )
    llc = [llc;[ec' nc(i)*ones(length(ec),1)]];
  end
%
% Find the centers:
  nTiles = length(llc);
  tc =  llc + ones(nTiles,1)*[1 1]*tileSize/2;
  
%
% Find the boundaries:
  outerEdge = tileSize/2 + edgeSize;
  tb = [tc(:,1)-outerEdge tc(:,1)+outerEdge tc(:,2)-outerEdge tc(:,2)+outerEdge];
  
  if( exist('tiles.csv')  ~= 2 && exist(parentMapName) == 2 )
    %
    % Create .csv file:
    tilesFid = fopen('tiles.csv','w');
    fprintf( tilesFid, 'TileName , Easting , Northing , %d\n', nTiles );
    %
    % Build grdcut argument strings and loop through:
    %
    % gmt6 has a bug.  The grdcut -N argument doesn't work anymore.
    %
    % gcutmap = sprintf('/usr/local/bin/gmt grdcut  ');
    %
    % Revert to gmt5, built locally:

    gcutmap = sprintf('/Users/rob/gmt5/install/bin/grdcut ');
    for( i=1:nTiles )
      tileName    = sprintf('MidAtlantic%d', i);
      grdTileName = strcat(tileName,'.grd');
      boTileName  = strcat(tileName,'.bo');
      if( exist(grdTileName) == 2 )
	%
	% If it exists, don't recompute it.
	fprintf('Tile %s already exists; don''t recreate.\n',grdTileName);
	%continue;
        return;
      end
      gcutNameArg = sprintf(' -G%s ', grdTileName);
      gcutBounds  = sprintf(' -R%d/%d/%d/%d -N ',tb(i,1),tb(i,2),...
			    tb(i,3),tb(i,4));
      %
      % grd-octree-maker takes easting limits in reverse order than grdcut:
      gcutBounds2 = sprintf(' -R %d %d %d %d',tb(i,2),tb(i,1),...
			    tb(i,3),tb(i,4));
      gcutCmd = [gcutmap parentMapName gcutNameArg gcutBounds];
      fprintf('%s', gcutCmd);
      fprintf(tilesFid, '%s , %d , %d\n', boTileName, tc(i,1), tc(i,2)); 
      tic;
      system(gcutCmd);
      fprintf('Created tile %s.\n', grdTileName);
      toc;
      %
      % Now create the .bo:
      boCmd = ['~/libtrnav/bin/grd-octree-maker ' tileName gcutBounds2];
      boCmd = ['~/libtrnav/bin/grd-octree-maker ' tileName ];  %It figures out the bounds.
      tic
      system(boCmd);
      fprintf('Created tile %s.\n', boTileName);
      toc;
      fprintf('\n\n');
    end
  else
    fprintf('tiles.csv exists, skip tile generation loop.\n');
  end
  
  Levels = [-105:1:-50];
  tic
  fprintf('Plotting Map\n');
  deci = 10;
  [cc hh] = contourf(downsample(plg.eastings,deci),downsample(plg.northings,deci),downsample(downsample(plg.z,deci)',deci)');
  clabel(cc,hh);
  axis equal
  set(gca,'ytick',[nc nc(end)+1500]);
  set(gca,'xtick',[ec ec(end)+1500]);
  toc
  grid on
  hold on

  if( exist( 'missionFileName') );
    %dataWP=corridor(['~/missions/',missionFileName]);
    dataWP=corridor(missionFileName);
  end
  
  linewidth=1;
  if( exist( 'dataWP' ) == 1 & min(size(dataWP)) > 1 )
    plot( dataWP(:,3), dataWP(:,2), 'r', 'linewidth', linewidth);
    lstr(length(lstr)+1,1) = {'Trackline'};
    plot( dataWP(:,3), dataWP(:,2), 'gd','linewidth', linewidth);
    lstr(length(lstr)+1,1) = {'Waypoints'};
    sizeDataWP = size(dataWP);
  end
  
  for( i=1:nTiles )
    locStr = sprintf('%d',i);
    tHandle = text(tc(i,1),tc(i,2),locStr);
    set(tHandle,'color','m');
    set(tHandle,'fontsize',12);
  end
