function bathymetry = extractBathymetry(lat, lon, bathy_db_path)
% Extract bathymetry data for specific lat/lon grid
% Author: Guillaume Bourdin
% Date: 2023-01-19
% Input:
%   - lat nxm <double> latitude for which to extract bathymetry
%   - lon nxm <double> longitude for which to extract bathymetry
%   - bathy_db_path <char> path to GEBCO bathymetry data (netCDF format)
%%
if nargin < 3
  error('Not enough input argument')
end
fprintf('Interpolating bathymetry ... ')

% load bathymetric database lat/lon
lat_bath_v = ncread(bathy_db_path, '/lat');
lon_bath_v = ncread(bathy_db_path, '/lon');

% wrap longitudes around 360°
foo_lon = lon;
foo_lonbath = lon_bath_v;
if (min(lon(:)) - 0.5) * (max(lon(:)) + 0.5) < 0
  foo_lon(foo_lon < 0) = foo_lon(foo_lon < 0) + 360;
  foo_lonbath(foo_lonbath < 0) = foo_lonbath(foo_lonbath < 0) + 360;
end

% reduce bathymetry lat/lon around output lat/lon to increase speed
idx_lat = lat_bath_v > min(lat(:)) - 0.5 & lat_bath_v < max(lat(:)) + 0.5;
idx_lon = foo_lonbath > min(foo_lon(:)) - 0.5 & foo_lonbath < max(foo_lon(:)) + 0.5;

% get elevation variable id in nectcdf
ncid = netcdf.open(bathy_db_path,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid, 'elevation');

% get lat columns to extract
start_col = find(idx_lat, 1, 'first') - 1;
ncol = find(idx_lat, 1,'last') - start_col;

% get lon rows to extract
foo = find(abs(diff([0; idx_lon; 0]))) - 1;
if size(foo, 1) == 4 % if image is over -180/180° meridian
  start_row1 = foo(1);
  nrow1 = foo(2) - start_row1;
  bath_subset1 = netcdf.getVar(ncid, varid,  [start_row1 start_col], [nrow1 ncol], 'double');
  start_row2 = foo(3);
  nrow2 = foo(4) - start_row2;
  bath_subset2 = netcdf.getVar(ncid, varid, [start_row2 start_col], [nrow2 ncol], 'double');
  bath_subset = [bath_subset1; bath_subset2]; 
else
  start_row = find(idx_lon, 1, 'first') - 1;
  nrow = find(idx_lon, 1,'last') - start_row;
  bath_subset = netcdf.getVar(ncid, varid, [start_row start_col], [nrow ncol], 'double');
end

% grid bathymetry lat/lon onto matrices
lat_bath_subset = repmat(lat_bath_v(idx_lat)', size(foo_lonbath(idx_lon),1), 1);
lon_bath_subset = repmat(foo_lonbath(idx_lon), 1, size(lat_bath_v(idx_lat), 1));

% interpolate bathymetry onto output lat/lon [linear, cubic, natural, nearest]
bathymetry = griddata(lon_bath_subset, lat_bath_subset, bath_subset, double(foo_lon), double(lat), 'cubic');

fprintf('Done\n')
