%
% PURPOSE:  Convert Dave's .grd map to a .mat.  The .grd may be in either UTM or
%           Geographic coordinates.  Based on loadGrd.m
%
% AUTHORS: (1) Rob McEwen, 8 July 2021.
%
% INPUTS:  ncname - a string containing the name of the .grd file, including
%                   the .grd extension.
%
%

function plg = cvtGrid(mapName,projection)

%inPath   = '/Users/rob/Maps/';
%outPath  = '/Users/rob/Documents/MATLAB/';


% if( strncmp(computer,'PCWIN',5))
%     pathRoot = '//atlas';
%     [status hname] = system('hostname');
%     if( regexp(hname,'NORTH') )
%       pathRoot = 'C:\dorado\data'
%     end
% elseif( strncmp(computer,'MAC',3))
%     pathRoot = '/Volumes';
% else
%     fprintf('Unknown host %s\n',computer);
%     return
% end
%
% Single quotes specifiy the variable name.
if(exist('mapName') == 0 || isempty(mapName))
    fprintf('Please pass in the map name in quotes.'\n');
    return
end

if( exist('mapName') )
  pathName = mapName;
else
  fprintf('Error.  Can''t find %s.\n', mapName);
  return
end
pathName

% else
%   pathName = [inPath,mapName];
%   if(exist(pathName) == 0)
%     fprintf('Error.  Could not find %s either locally or in %s.\n', mapName, inPath);
%     return
%   end
% end

% description=ncreadatt(pathName,'/','source');
% [token remain]=strtok(description)
% [token remain]=strtok(remain)

% if(strcmp(token,'Geographic'))

if(strncmpi(projection,'Geographic',3))

  fprintf('Converting %s from Geographic .grd to UTM .mat %s.\n',pathName);

  tic

  [plg.lon plg.lat plg.z] = grdread2rob(pathName);
  plg.z = double(plg.z);

  toc

  fprintf('Finished loading %s.\n', pathName);

  %
  % Here I'm using the mapping convention; x is eastings, y northings.
  %

  %
  % Convert to UTM for plotting
  %
  tic
  fprintf('Converting to UTM.\n');
  len = length(plg.lon);
  indx = find( plg.lon > 180 );
  plg.lon(indx) = plg.lon(indx)-360;
  zone = utmzone(plg.lat(1), plg.lon(1));
  zoneEnd = utmzone(plg.lat(end), plg.lon(end));
  if( ~strcmp(zone,zoneEnd))
    fprintf('This data crosses UTM zone boundary. Using zone %s\n', ...
        zone);
  end
  mstruct = defaultm('utm');
  mstruct.zone = zone;
  mstruct  = defaultm(mstruct);
  %
  % This gets around mfwdtran's requirement that the x and y axes be
  % matrices of the same dimension of z.
  %
  [plg.glon,plg.glat]=meshgrid(plg.lon,plg.lat); 
  %
  % LL2utm is ryjo's code that runs mfwdtran.  There is also a ll2utm on the
  % File Exchange that transforms lat lon to utm using first principles
  % without dependence on the mapping toolbox.
  %
  [plg.eastings,plg.northings]=LL2utm(plg.glon,plg.glat);
  %rmfield(plg,'gx')
  %rmfield(plg,'gy')
  plg.glon = 0;
  plg.glat = 0;
  plg.eastings=plg.eastings(1,:); 
  plg.northings=plg.northings(:,1);
  toc
  fprintf('Finished UTM conversion.\n');
  %
elseif(strncmpi(projection,'UTM',3))

  fprintf('Converting UTM .grd to .mat %s.\n',pathName);

  tic

  [plg.eastings plg.northings plg.z] = grdread2rob(pathName);
  plg.z = double(plg.z);

  toc

else
  
  fprintf('Conversion from %s  Projection not yet implemented.\n');
  plg = 0;
  return
end  
% Extract the file name. Like baseName.
i=1;
remain = pathName;
while(~isempty(remain))
  [token, remain]=strtok(remain,'/')
  i = i+1;
  if( i > 50 )
    fprintf('Error.  Path to pathName is too long.\n')
    return
  end
end
token
fileName = replace(token,'.grd','.mat');
fileName
fprintf('Saving file to disk\n');
tic
if(strncmpi(projection,'Geographic',3))
  save(fileName,'plg','zone','zoneEnd','mstruct','pathName');
elseif(strncmpi(projection,'UTM',3))
  save(fileName,'plg','pathName');
end
toc
fprintf('Save completed.  End cvtGrd.\n');
return
