Bathymetry plotting along a transect

[kakearney]

Both bathymetry_section and bathymetry_chord suffer from issues regarding how they order the output points when plotting against a transect.

For bathymetry_section, there are two issues. The first is that it uses the coordinates of the closest-matching-grid cell, rather than the input coordinates, along the x axis. Unless the input line is parallel to the specified axis (i.e. an east/west line paired with lon, or a north/south line paired with lat), this may cause strong reduction in the plotted x-coordinate resolution relative to the input coordinates. The points are also plotted in the order they appear in the input, rather than in increasing order with longitude and/or latitude, which can be an issue for any set of points whose x-coordinate dataset (x for xref='lon' or y for xref='lat') is not monotonically increasing.

For bathymetry_chord, the output can get even messier. The general idea behind this function appears to be to plot the bathymetry along a transect line. However, it only does so in a very limited context: when the line is parallel to the specified x axis, and the specified width is approximately equal to or less than to the resolution of the bathymetry dataset. Anything else results in a mess of lines doubling back on themselves, due to the way grid cell within the polygon are collected.

A demonstration:

% A few lines near the Aleutian Islands

latlon = [...
       53.51       176.23
       53.74         -177
       55.055       179.68
       50.623       -179.8
       53.165       171.64
       51.687       178.25];
latlon(:,2) = wrapTo360(latlon(:,2));

latlim = minmax(latlon(:,1), 'expand');
lonlim = minmax(latlon(:,2), 'expand');

% Read bathymetry

bdir = '/Volumes/LacieKK2019/SmithSandwell/topo_20.1.nc'; % user-specific
[z,lt,ln] = bathymetry_extract(bdir, [latlim lonlim]);

cmap = get(0, 'defaultAxesColorOrder');

% Plot map of bathymetry

figure; axes;
pcolor(ln,lt,z');
shading flat;
colormap('gray');
hold on;
plot(latlon(1:2,2), latlon(1:2,1), 'color', cmap(1,:));
plot(latlon(3:4,2), latlon(3:4,1), 'color', cmap(2,:));
plot(latlon(5:6,2), latlon(5:6,1), 'color', cmap(3,:));
set(gcf, 'color', 'w');

d = {'lat', 'lon'};

% Plot transects

for ii = 1:3

    idx = (ii-1)*2+[1 2];
    [y, x] = interpm(latlon(idx,1), latlon(idx,2), 1/60);

    zline = interp2(ln,lt,z',x,y, 'nearest');

    h = plotgrid('size', [4 2], 'sp', 0.05);

    for id = 1:2
        % Wider chord
        axes(h.ax(1,id));
        bathymetry_chord(bdir, x(1), y(1), x(end), y(end), d{id},  0, 1/20);
        % Narrower chord
        axes(h.ax(2,id));
        bathymetry_chord(bdir, x(1), y(1), x(end), y(end), d{id},  0, 1/60);
        % Section
        axes(h.ax(3,id));
        bathymetry_section(bdir, x, y, d{id}, 0);
        % Section (my way)
        switch id
            case 1
                plot(h.ax(4,id), y, zline, 'color', 'k');
            case 2
                plot(h.ax(4,id), x, zline, 'color', 'k');
        end
    end

    set(findall(h.ax, 'type', 'line'), 'linewidth', 1, 'color', cmap(ii,:));
    set(h.ax(:,1), 'xlim', minmax(y));
    set(h.ax(:,2), 'xlim', minmax(x), 'yaxisloc', 'right');
    set(h.ax, 'ylim', minmax(zline));
    xlabel(h.ax(end,1), 'Latitude');
    xlabel(h.ax(end,2), 'Longitude');
    set(h.fig, 'color', 'w');
    set(h.ax(1:end-1,:), 'xticklabel', '');

    h.lbl = multitextloc(h.ax(:,1), {...
        'bathymetry_chord (width=1/20)'
        'bathymetry_chord (width=1/60)'
        'bathymetry_section'
        'the transect'}, 'northwestoutsideabove');
    set(h.lbl, 'interpreter', 'none');

end

(Note: this code uses some of my personal toolbox: plotgrid, minmax, multitextloc)

[lnferris]

I have removed bathymetry_chord(). You are correct that it serves narrow purpose and I think it will get users into more trouble than it is worth. For bathymetry_search I have replaced the search for nearest grid cells with grid densification and interpolation as you have suggested. I think the points are ordered by the chosen reference variable rather than the order of the input (see below) but please let me know if I have misunderstood this comment.

% Order bathymetry_section by xref.
[xvar,xvar_inds] = sort(xvar);
bath_section = bath_section(xvar_inds);
lon_section = lon_section(xvar_inds);
lat_section = lat_section(xvar_inds);
time_section = time_section(xvar_inds);

I am closing this issue but feel free to reopen it if these fixes are not sufficient. Thank you for the helpful advice on this.