function plotBgrid_obj(d)
% plotBgrid_obj(d)
% where d is a structure array with fields
% A         [lat x lon] grid of data; missing or 'no transformation'
%               values are now NaN
% lat       latitudes (decimal degrees) corresponding to 
%               rows of A
% lon       longitudes (decimal degrees) corresponding to 
%               columns of A
% fbase     file base with naming convention
%               old_datum.new_datum.region

addpath Q:\VERTCON3\VERTCON3toolbox\cptcmap
t1 = 'VERTCON v3.0';

% plot title
trans = cell2struct(split(d.fbase,'.'), {'oldDat','newDat','region'});
if isfield(d, 'grid')
tstr = sprintf('%s %s minus %s OHT(%d sec) %s-%s %s-f%03d',...
    t1, trans.newDat, trans.oldDat, d.grid, trans.region, d.zone,...
    d.pltType, d.filter);
else
tstr = sprintf('%s %s minus %s OHT %s-%s %s-f%03d',...
    t1, trans.newDat, trans.oldDat, trans.region, d.zone,...
    d.pltType, d.filter);
end

[X,Y] = meshgrid(d.lon,d.lat);
figure;
s = surf(X,Y,d.A);
title(tstr);
axis equal
view(0,90)
s.EdgeColor = 'none';
% xlabel('E longitude (deg)')
% ylabel('N latitude (deg)')
c = colorbar('southoutside');
c.Label.String = 'meters';

% GMT 'rainbow' colormap, with 8 discrete colors
% and tick marks/labels at each of the color transitions
cptcmap('gmt_rainbow','ncol',8);
cl = caxis;
tk.m = 100;
dcl = diff(cl);
%if dcl<0.5; tk.m = 100; end
tk.inc = ceil(tk.m*(dcl/8))/tk.m;
tk.start = floor(cl(1)/tk.inc)*tk.inc;
%tk.start = floor(cl(1)/tk.m)/tk.m;
caxis([tk.start tk.start+tk.inc*8]);
c.XTick = tk.start:tk.inc*2:tk.start+tk.inc*8;
% a = pbaspect; % number of tick labels determined by plot
%         % aspect ratio
% if a(1)<a(2)
%     c.XTick = tk.start:tk.inc*2:tk.start+tk.inc*8;
% else
%     c.XTick = tk.start:tk.inc:tk.start+tk.inc*8; 
% end

% dcl = diff(cl);
% tk.inc = ceil(10*(dcl/8))/10;
% tk.start = floor(cl(1)*10)/10;
% caxis([tk.start tk.start+tk.inc*8]);
% c.XTick = tk.start:tk.inc:tk.start+tk.inc*8; 

% set axes background to light gray, for printing
set(gca,'color',[.9 .9 .9])
set(gcf,'color',ones(3,1))
set(gcf, 'InvertHardcopy', 'off')

% set axes limits to boundaries of full grid
axis(d.gridLim);
disp([d.gridLim(1:2) d.gridLim(3:4)]);

% overlay the shoreline for land areas
zmax = max(max(s.CData));
plotBorders3d(trans.region, zmax);
end

function plotBorders3d(areaName, z)

% load list of filesnames
switch lower(areaName)
%     case 'conus'
%         basemapName = 'continental us';
%     case 'alaska'
%         basemapName = 'alaska';
%         %areaList = {'AK'}; % seemingly duplicated by 'borders'
    case 'pr'
        areaList = {'PR'};
    case 'vi'
        areaList = {'VQ'};
    case 'hawaii'
        areaList = {'HI'};
    case 'as'
        areaList = {'AS'};
    case 'guam'
        areaList = {'GQ'};
    case 'cnmi'
        areaList = {'CQ'};
    otherwise
        error('areaName not recognized\n does not work for conus or alaska')
end

% % draw state boundaries
% if ~isempty(basemapName)
%     borders(basemapName,'k');
% end

hold on; % hold figure

% supplement with GMT files, if needed for area
bb = strrep(userpath,';','\');
str = sprintf('%s//boundaries/',bb);
for k = 1:length(areaList)
    % get list of files for each area
    ff = dir(sprintf('%s%s*.gmt',str,areaList{k}));
    for i = 1:length(ff)
        % load and plot file
        x{i} = load([str ff(i).name]);
        plot3(x{i}(:,1), x{i}(:,2), ...
            ones(size(x{i},1),1)*z, 'k'); hold on
    end
end
end