Code for merging multibeam bathymetry data from the South Orkney Islands with background bathymetry grids

merge_grids_20190731.m

% MERGE_GRIDS_20190731
% by Povl Abrahamsen, British Antarctic Survey
% Updated 31 July 2019 (code cleaned up, comments added)
%
% This script merges gridded multibeam bathymetry from the South Orkney
% Islands with a background grid, in this case from GEBCO_2014.

%% define input and output file names/prefixes
grid_in='clean_swath_grid_20190731.grd';
grid_out='new_orkney_topo_20190731';

%% set up mbsystem prefix (system-dependent)
mbgrid_prefix = '/sw/bin/'; % Povl's configuration on MacOS, MB-system installed using Fink (http://www.finkproject.org/)
% mbgrid_prefix = 'setenv LD_LIBRARY_PATH /sw/lib;/sw/bin/'; % alternative: force using system NetCDF library, not Matlab's!
% mbgrid_prefix = ''; %default

%% note: input grid was generated using the following command:
% system([mbgrid_prefix,...
%     'mbgrid -Iupdated_swath_20190731.mb-1 -E100/100/meters -R-47/-37/-63/-59 ',...
%     '-Oclean_swath_grid_20190731 -G3 -A2 -W4 -N -M -C5/1 -T25 -P1 -F1 -V']);

%% load input swath grid
x=ncread(grid_in,'x');
y=ncread(grid_in,'y');
z=ncread(grid_in,'z')';

%% load background bathymetry (in this case, GEBCO_2014)
% grid downloaded from https://www.gebco.net - then subsetted using nco (http://nco.sourceforge.net/):
% ncks -d lon,15960,17161 -d lat,3240,3721 GEBCO_2014_2D.nc GEBCO_2014_2D_-47.0_-63.0_-37.0_-59.0.nc
grid_bg=fullfile('GEBCO_2014_2D_-47.0_-63.0_-37.0_-59.0.nc');
bgbathy.lon=ncread(grid_bg,'lon');
bgbathy.lat=ncread(grid_bg,'lat');
bgbathy.depth=ncread(grid_bg,'elevation')';

%% interpolate background bathymetry onto swath grid
bg_depth=interp2(bgbathy.lon*.5,bgbathy.lat',bgbathy.depth,...
    x*.5,y','cubic');

%% extend coordinates by one grid cell in each direction
x0=x(1)-median(diff(x));
xend=x(end)+median(diff(x));
xtemp=[x0;x;xend];

y0=y(1)-median(diff(y));
yend=y(end)+median(diff(y));
ytemp=[y0;y;yend];

%% calculate difference between swath and background; pad with zeros around edges
ztemp=zeros(size(z)+2);
ztemp(2:end-1,2:end-1)=z-bg_depth;
[xtemp,ytemp]=meshgrid(xtemp,ytemp);
goodind=find(~isnan(ztemp));
ztemp=ztemp(goodind);
xtemp=xtemp(goodind);
ytemp=ytemp(goodind);

%% write differences to temporary ASCII file - and temporary file list
fid=fopen('temp_file.mb162','wt'); % XYZ (lon lat depth) ASCII soundings, generic
fprintf(fid,'%.6f %.6f %d\n',[xtemp,ytemp,round(ztemp)]');
fclose(fid);

fid=fopen('temp_filelist.mb-1','wt');
fprintf(fid,'temp_file.mb162 162\n');
fclose(fid);

%% grid differences using mbgrid

%warning: the next line might take a while to execute!
system(sprintf([mbgrid_prefix,...
    'mbgrid -Itemp_filelist.mb-1 -E%.12f/%.12f/degrees ',...
    '-R%.8f/%.8f/%.8f/%.8f -Odiff_gridded -G3 -A2 -C9999/3 -T50 -W1 -V'],...
    median(diff(x))*10,median(diff(y))*10,x0,xend,y0,yend));

%% load difference grid
diff_grid.x=ncread('diff_gridded.grd','x');
diff_grid.y=ncread('diff_gridded.grd','y');
diff_grid.z=ncread('diff_gridded.grd','z')';

%% interpolate difference grid back onto swath grid
diff_regridded=interp2(diff_grid.x,diff_grid.y,diff_grid.z,x,y');

%% subtract difference grid from background grid and insert into swath grid gaps!
z(isnan(z))=bg_depth(isnan(z))-diff_regridded(isnan(z));

%% write grid to NetCDF file
create_gmt_grid([grid_out,'.grd'],x,y,z);

%% alternative: write grid to ArcGIS-compatible ASCII file
% create_arc_grid([grid_out,'.asc'],x,y,z);


function create_gmt_grid(filename,x,y,z)
% write a GMT-compatible NetCDF grid

if nargin<1 || ~ischar(filename)
    error('First argument must be file name');
end

nc_create_empty(filename,nc_clobber_mode); 
% or bitor(nc_clobber_mode,nc_64bit_offset_mode) for really big files

nc_attput(filename,nc_global,'title','Bathymetry');
nc_attput(filename,nc_global,'history',...
    sprintf('Written by Matlab on %s',datestr(now)));

nc_add_dimension(filename,'x',length(x));
nc_add_dimension(filename,'y',length(y));

nc_addvar(filename,struct('Name','x','Nctype',nc_double,'Dimension',{{'x'}}));
nc_attput(filename,'x','long_name','Longitude');

nc_addvar(filename,struct('Name','y','Nctype',nc_double,'Dimension',{{'y'}}));
nc_attput(filename,'y','long_name','Latitude');

nc_addvar(filename,struct('Name','z','Nctype',nc_float,...
    'Dimension',{{'y','x'}}));
nc_attput(filename,'z','long_name','Topography (m)');
nc_attput(filename,'z','_FillValue',nan);

nc_varput(filename,'x',x);
nc_varput(filename,'y',y);
nc_varput(filename,'z',z);

end

function create_arc_grid(filename,x,y,z)
% write an ArcGIS-compatible ASCII grid

%% interpolate longitudes onto same resolution as latitude
xnew=[x(1):median(diff(y)):x(end)]';
prevwarn=warning('off','MATLAB:interp1:NaNinY');
znew=interp1(x,z',xnew,'linear')';
warning(prevwarn);
clear prevwarn

%% write ascii grid file
fid=fopen(filename,'wt');
fprintf(fid,'ncols %d\n',length(xnew));
fprintf(fid,'nrows %d\n',length(y));
fprintf(fid,'xllcenter %.4f\n',xnew(1));
fprintf(fid,'yllcenter %.4f\n',y(1));
fprintf(fid,'cellsize %.4f\n',median(diff(xnew)));
fprintf(fid,'nodata_value %d',-32768);
znew(isnan(znew))=-32768;
for n=length(y):-1:1
    fprintf(fid,'\n');
    fprintf(fid,'%d ',round(znew(n,:)));
end
fclose(fid);

end

merge_grids_20190731_gebco2019.m

% MERGE_GRIDS_20190731_GEBCO2019
% by Povl Abrahamsen, British Antarctic Survey
% Updated 31 July 2019 (code cleaned up, comments added)
%
% This script merges gridded multibeam bathymetry from the South Orkney
% Islands with a background grid, in this case from GEBCO_2019.

%% define input and output file names/prefixes
grid_in='clean_swath_grid_20190731.grd';
grid_out='new_orkney_topo_20190731_gebco2019';

%% set up mbsystem prefix (system-dependent)
mbgrid_prefix = '/sw/bin/'; % Povl's configuration on MacOS, MB-system installed using Fink (http://www.finkproject.org/)
% mbgrid_prefix = 'setenv LD_LIBRARY_PATH /sw/lib;/sw/bin/'; % alternative: force using system NetCDF library, not Matlab's!
% mbgrid_prefix = ''; %default

%% note: input grid was generated using the following command:
% system([mbgrid_prefix,...
%     'mbgrid -Iupdated_swath_20190731.mb-1 -E100/100/meters -R-47/-37/-63/-59 ',...
%     '-Oclean_swath_grid_20190731 -G3 -A2 -W4 -N -M -C5/1 -T25 -P1 -F1 -V']);

%% load input swath grid
x=ncread(grid_in,'x');
y=ncread(grid_in,'y');
z=ncread(grid_in,'z')';

%% load background bathymetry (in this case, GEBCO_2019)
% file downloaded from https://www.gebco.net
grid_bg=fullfile('GEBCO_2019_-47.01_-59.0_-37.0_-63.01.nc');
bgbathy.lon=ncread(grid_bg,'lon');
bgbathy.lat=ncread(grid_bg,'lat');
bgbathy.depth=double(ncread(grid_bg,'elevation'))';

%% interpolate background bathymetry onto swath grid
bg_depth=interp2(bgbathy.lon*.5,bgbathy.lat',bgbathy.depth,...
    x*.5,y','cubic');

%% extend coordinates by one grid cell in each direction
x0=x(1)-median(diff(x));
xend=x(end)+median(diff(x));
xtemp=[x0;x;xend];

y0=y(1)-median(diff(y));
yend=y(end)+median(diff(y));
ytemp=[y0;y;yend];

%% calculate difference between swath and background; pad with zeros around edges
ztemp=zeros(size(z)+2);
ztemp(2:end-1,2:end-1)=z-bg_depth;
[xtemp,ytemp]=meshgrid(xtemp,ytemp);
goodind=find(~isnan(ztemp));
ztemp=ztemp(goodind);
xtemp=xtemp(goodind);
ytemp=ytemp(goodind);

%% write differences to temporary ASCII file - and temporary file list
fid=fopen('temp_file.mb162','wt'); % XYZ (lon lat depth) ASCII soundings, generic
fprintf(fid,'%.6f %.6f %d\n',[xtemp,ytemp,round(ztemp)]');
fclose(fid);

fid=fopen('temp_filelist.mb-1','wt');
fprintf(fid,'temp_file.mb162 162\n');
fclose(fid);

%% grid differences using mbgrid

%warning: the next line might take a while to execute!
system(sprintf([mbgrid_prefix,...
    'mbgrid -Itemp_filelist.mb-1 -E%.12f/%.12f/degrees ',...
    '-R%.8f/%.8f/%.8f/%.8f -Odiff_gridded -G3 -A2 -C9999/3 -T50 -W1 -V'],...
    median(diff(x))*10,median(diff(y))*10,x0,xend,y0,yend));

%% load difference grid
diff_grid.x=ncread('diff_gridded.grd','x');
diff_grid.y=ncread('diff_gridded.grd','y');
diff_grid.z=ncread('diff_gridded.grd','z')';

%% interpolate difference grid back onto swath grid
diff_regridded=interp2(diff_grid.x,diff_grid.y,diff_grid.z,x,y');

%% subtract difference grid from background grid and insert into swath grid gaps!
z(isnan(z))=bg_depth(isnan(z))-diff_regridded(isnan(z));

%% write grid to NetCDF file
create_gmt_grid([grid_out,'.grd'],x,y,z);

%% alternative: write grid to ArcGIS-compatible ASCII file
% create_arc_grid([grid_out,'.asc'],x,y,z);


function create_gmt_grid(filename,x,y,z)
% write a GMT-compatible NetCDF grid

if nargin<1 || ~ischar(filename)
    error('First argument must be file name');
end

nc_create_empty(filename,nc_clobber_mode); 
% or bitor(nc_clobber_mode,nc_64bit_offset_mode) for really big files

nc_attput(filename,nc_global,'title','Bathymetry');
nc_attput(filename,nc_global,'history',...
    sprintf('Written by Matlab on %s',datestr(now)));

nc_add_dimension(filename,'x',length(x));
nc_add_dimension(filename,'y',length(y));

nc_addvar(filename,struct('Name','x','Nctype',nc_double,'Dimension',{{'x'}}));
nc_attput(filename,'x','long_name','Longitude');

nc_addvar(filename,struct('Name','y','Nctype',nc_double,'Dimension',{{'y'}}));
nc_attput(filename,'y','long_name','Latitude');

nc_addvar(filename,struct('Name','z','Nctype',nc_float,...
    'Dimension',{{'y','x'}}));
nc_attput(filename,'z','long_name','Topography (m)');
nc_attput(filename,'z','_FillValue',nan);

nc_varput(filename,'x',x);
nc_varput(filename,'y',y);
nc_varput(filename,'z',z);

end

function create_arc_grid(filename,x,y,z)
% write an ArcGIS-compatible ASCII grid

%% interpolate longitudes onto same resolution as latitude
xnew=[x(1):median(diff(y)):x(end)]';
prevwarn=warning('off','MATLAB:interp1:NaNinY');
znew=interp1(x,z',xnew,'linear')';
warning(prevwarn);
clear prevwarn

%% write ascii grid file
fid=fopen(filename,'wt');
fprintf(fid,'ncols %d\n',length(xnew));
fprintf(fid,'nrows %d\n',length(y));
fprintf(fid,'xllcenter %.4f\n',xnew(1));
fprintf(fid,'yllcenter %.4f\n',y(1));
fprintf(fid,'cellsize %.4f\n',median(diff(xnew)));
fprintf(fid,'nodata_value %d',-32768);
znew(isnan(znew))=-32768;
for n=length(y):-1:1
    fprintf(fid,'\n');
    fprintf(fid,'%d ',round(znew(n,:)));
end
fclose(fid);

end