% Example of usage of the Chesapeake Bay atlas (https://doi.org/10.17882/99441)
% MACAN webinar of 2024-09-05, pst-laurent@vims.edu
%
% A transect along a tributary of the Bay is visualized with surface/bottom
%   conditions and the bottom topography.
 
  clear;     % Empty the workspace.
  close all; % Close all existing figures (if there are any).

% Read from the atlas the different variables that we will need:

  atla = 'atlas_chesbay_v20240319.nc';
  long = ncread( atla, 'longitude'  ); % Size: 348
  lati = ncread( atla, 'latitude'   ); % Size:     567
  topo = ncread( atla, 'topography' ); % Size: 348x567 (units: meters).
  topo = double( topo );

% Select the variable and the MONTh that you want to illustrate. Here, TA is
%   selected but many other variables are available in the atlas.

  mont = 8;   % August (8th month of the year).
  bott = ncread( atla, 'TA_bottom',  [1, 1, mont], [inf, inf, 1] ); % Size: 348x567 (umol/kg).
  surf = ncread( atla, 'TA_surface', [1, 1, mont], [inf, inf, 1] );
  clear atla; % No longer necessary.

% Read the longitudes/latitudes corresponding to the COastLIne of the Bay:
 
  coli = ncread( 'coastlines_chesbay.nc', 'coli' ); % Size: 314078x2

% Read the COORdinates (longitudes/latitudes) corresponding to the Bay's MAin
%   STem as well as a few tributaries. Each of these files has 3 columns:
%   (1) longitude, (2) latitude, (3) position in kilometers relative to the
%   mouth of the Bay or tributary.

  coor = 'lon_lat_transects_chesbay.nc';
  mast = ncread( coor, 'main_stem_of_bay'   ); % 598x3
  york = ncread( coor, 'york_river'         ); % 112x3
  poto = ncread( coor, 'potomac_river'      ); % 340x3
  patu = ncread( coor, 'patuxent_river'     ); %  83x3
  rapp = ncread( coor, 'rappahannock_river' ); % 151x3
  clear coor; % No longer necessary.

% Illustrate on a map of the Bay the different transects for which coordinates
%   are available.
%   First draw the COastLInes, then set the geographical extent of the map,
%   and finally, overlay the lon/lat corresponding to the various transects.

  figure;
    hand = plot( coli(:, 2), coli(:, 1) ); axis ij;
    set( hand, 'color', 0.7 * [1 1 1] ); % RGB code for 'light gray'.
    set( gca,  'xlim', [36.8, 39.6], 'ylim', [- 77.3, - 75.5] );
  hold on;
    hand = plot( mast(:, 2), mast(:, 1), ...
                 york(:, 2), york(:, 1), ...
                 poto(:, 2), poto(:, 1), ...
                 patu(:, 2), patu(:, 1), ...
                 rapp(:, 2), rapp(:, 1) );
    set( hand, 'linewidth', 9 );         % Make the line thicker.
  hold off;
  clear hand; % No longer necessary.

% Select one specific transect that will be visualized:

  tran = mast; % Main stem of the Bay.
% tran = york; % York         River.
% tran = poto; % Potomac      River.
% tran = patu; % Patuxent     River.
% tran = rapp; % Rappahannock River.

% Interpolate the BATHymetry and the VaRiable (at the Surface/Bottom) along the
%   transect chosen above. The atlas information is defined on a regular lon/lat
%   grid and therefore a trivial bilinear interpolation is all that is needed.

  bath = interp2( long, lati, topo', tran(:, 1), tran(:, 2) ) * (- 1.);
  vr_b = interp2( long, lati, bott', tran(:, 1), tran(:, 2) ); % VaR. @ bott.
  vr_s = interp2( long, lati, surf', tran(:, 1), tran(:, 2) ); % VaR. @ surf.

% Create a closed POLygon defined by x,y points; the polygon will appear at the
%   bottom of the figure to illustrate the bathymetry. The horizontal (x) axis
%   of the figure will be the position in kilometers from the mouth:

  xpol = [ tran( :, 3 ); tran( end, 3 ); tran( 1, 3 ); tran( 1, 3 ) ];
  ypol = [ bath;         max( bath ) * 1.01 * [ 1; 1]; bath( 1    ) ]; 

% Note the MInimum/MAximum values of the variable along this transect (to set
%   the colorscale below):

  mima = [ min([ vr_b(:); vr_s(:) ]) - 1.e-3, ...
           max([ vr_b(:); vr_s(:) ]) + 1.e-3 ];
  
% Create the final figure. Select the 'jet' color palette with a size of 64 x 3
%   (the three columns corresponding to the Red/Green/Blue channels).

  figure; colormap( jet( 64 ) );
  hand = patch( xpol, ypol, 0.7 * [1 1 1] ); axis ij; % 0.7 means 'light gray'.
  set( gca, 'xlim', [min(xpol), max(xpol)], ...       % Set min/max x/y.
            'ylim', [       0., max(ypol)] );
  xlabel('Distance from mouth (km)');
  ylabel('Depth (m)'               );
  hold on; % Begin the overlay.

%   Overlay the selected variable as a series of colored circles matching a
%     color legend (colorbar).

    for i = 1 : length( bath ) % Loop through each point of the transect...

%     If a given point of the transect falls on land (and thus no values are
%       available at that point), then skip that point ('continue').

      if ( isnan( vr_b( i ) ) ); continue; end
  
%     Plot the bottom concentration near the sea floor (i.e., at 90% of the
%        bottom depth 'bath') as a circle ('o').

      hand = plot( tran( i, 3 ), bath( i ) * 0.9, 'o' );

%     Assign the color (Red/Green/Blue values) of that circle, knowing that the
%        color palette has 64 values and that its extremities correspond to the
%        MIn/MAx identified above:

      rgb  = interp1( (mima(1) : (mima(2) - mima(1)) / 63. : mima(2))', ...
                      jet( 64 ), vr_b( i ) );
      set( hand, 'markersize',      6,   ...
                 'markeredgecolor', rgb, ...
                 'markerfacecolor', rgb ); clear hand rgb;

%     Same, but for the surface (y = 0, vr_s instead of vr_b):
  
      hand = plot( tran( i, 3 ), 0.,              'o' );
      rgb  = interp1( (mima(1) : (mima(2) - mima(1)) / 63. : mima(2))', ...
                      jet( 64 ), vr_s( i ) );
      set( hand, 'markersize',      6,   ...
                 'markeredgecolor', rgb, ...
                 'markerfacecolor', rgb ); clear hand rgb;
    end; clear i;

%   Enforce a colorscale with the MIn/MAx defined previously:

    caxis( mima ); colorbar;

%   Add labels for longitude/latitude along the transect:

    if ( (max(tran( :, 1 )) -  min(tran( :, 1 ))) * cosd( 38. ) ...
       > (max(tran( :, 2 )) -  min(tran( :, 2 ))) ) % East -west  transect.
      post = tran( :, 1 ); unit = 'W';              % POSition along Transect.
    else                                            % North-south transect.
      post = tran( :, 2 ); unit = 'N';              % POSition along Transect.
    end
    labl = ( ceil( min( post ) * 10.) / 10. : 0.1 ...
           : floor(max( post ) * 10.) / 10. )';     % LABeLs to be added.

    for ilab = 1 : length( labl )        % Add labels one by one.
      clos = abs( post - labl( ilab ) ); % Find transect point CLOSest to label.
      clos = find( clos == min( clos(:) ) );
      clos = clos(1);
      hand = text( tran( clos, 3 ), max( bath(:) ) * 0.95, ...
                   [ num2str( labl( ilab ) ) '^{o}' unit ] );
      set( hand, 'rotation', 90 );       % Rotate label by 90deg.
      clear clos hand;
    end; clear ilab post unit labl;

  hold off; % Stop the overlay.

  clear bath vr_b vr_s xpol ypol mima; % No longer necessary.