%
% 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. 
%
% Rob M., 12 July 2018
% Rob M., 29 Mar 2023
%
% This requires that the .grd map and a corresponding mission (.cfg or .cfg.pp)
% be in the local directory. Typing "ls *[gp]" should show them.
%
% Example:  Dave's EMARK_vents map
%
%  To find the boundaries of the .grd use: 
%
%  gmt grdinfo EMARK_vents_2023/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd
%  
% Then, choose tile size.  The non-overlapped regions are tileSize x tileSize
% meters. The whole tile size is then (tileSize+edgeSize) x (tileSize+edgeSize).
%
% Start with arrays of the lower left corners of the NON-OVERLAPPED regions:
%
  if(1)
    %
    % Standard Boulder-Crossing test:
    %
    % On the Matlab cmd line: 
    %
    % close all;clear all;missionFileName='20140529_PortugueseLedge_CrossRoute_1v3.cfg';tileMaker
    %
    parentMapName = 'PortugueseLedge20080424TopoUTM_NoNan.grd';
    tileNameStr =   'PortugueseLedgeTile';
    tileSize = 500;
    edgeSize = 250;
    ec = [ 594000:tileSize: 595000]; % Easting corners
    nc = [4061000:tileSize:4062000]; % Northing corners
    %
    % For plotting:
    deci = 2;
    Levels = [-105:1:-50];
    ytickVal = [4061000 4061500 4062000 4062500];
    xlimVal=[ 594000  595500];
    ylimVal=[4061000 4062500];
  else
    % On the Matlab cmd line: 
    %
    % close all;clear all;missionFileName='EMARK_13v3.cfg.pp';tileMaker
    %
    parentMapName = '~/Maps/EMARK_vents_2023_Test/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd';
    tileNameStr = 'MidAtlantic';
    tileSize = 1000;
    edgeSize = 250;
    ec = [ 499000:tileSize: 504000]; % Easting corners
    nc = [2593000:tileSize:2598000]; % Northing corners
    %
    % For plotting:
    deci = 10;
    xtickVal = [ec ec(end)+1500];
    ytickVal = [nc nc(end)+1500];
  end

  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', '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 );
      return
    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('%s%d', tileNameStr, 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
  
  tic
  fprintf('Plotting Map\n');
  if(exist('Levels'))
    [cc hh] = ...
        contourf(downsample(plg.eastings,deci),downsample(plg.northings,deci),downsample(downsample(plg.z,deci)',deci)', Levels);
  else
    [cc hh] = ...
        contourf(downsample(plg.eastings,deci),downsample(plg.northings,deci),downsample(downsample(plg.z,deci)',deci)');
  end
  clabel(cc,hh);
  axis equal
  if(exist('xtickVal'))
    set(gca,'xtick',xtickVal);
  end
  if(exist('ytickVal'))
    set(gca,'ytick',ytickVal);
  end
  if(exist('xlimVal'))
    set(gca,'xlim',xlimVal);
  end
  if(exist('ylimVal'))
    set(gca,'ylim',ylimVal);
  end
  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

  fprintf('There will be %d tiles.\n', nTiles);
  
  
