if ~exist('raw')
    % load the ngdc file
    raw = load('cblastbathyraw');
    raw.h = -raw.h; % positive for water 
end

if ~exist('poly')
    % interactively create a polygon defining the area of bathymetry to patch
    poly = ginput;
end
w = poly([1:end 1],1)+sqrt(-1)*poly([1:end 1],2);

% find the points inside the polygon
ind = inside(raw.lon(:)+sqrt(-1)*raw.lat(:),w);

% don't do land points
trim = find(raw.h(ind)<0);
ind(trim) = [];

% convert pixel units to lon,lat
LON = c(1)*x-c(2);
LAT = d(1)*y-d(2);
% convert feet to meters
H = h*25.4*12/1000;

% interp new digitized bathymetry to original (raw) grid points that are
% inside the selected polygon
raw.hnew = raw.h;
new = griddata(LON,LAT,H,raw.lon(ind),raw.lat(ind));
raw.hnew(ind) = new;