# VMT_unitQcont

## PURPOSE

VMT_UNITQCONT applies unit discharge continuity correction MCS velocities

## SYNOPSIS

function A = VMT_unitQcont(A,V,z)

## DESCRIPTION

     VMT_UNITQCONT applies unit discharge continuity correction MCS velocities
     This routine computes the correction, then applies it to the North and
     East velocity components of each ensemble and bin for the ADCP data. This
     routine runs BEFORE the data are gridded to the mean cross section, so
     if used, the correction applies to all subsequent computations (e.g., ZSD
     and Rozovskii decomposition).

     DISCUSSION
     In theory, the change in elevation in the streamwise direction
     should be zero (since we are creating a plane in the lateral direction).
     However, in practice, this may not be the case. Hopefully, the user has
     chosen an appropriate cross section location where the bed is not varying
     dramatically in the streamwise direction, but even still, some error may
     exist. One way to approach this potential error is to account for any
     changes in the mass flux between the sample location where data were
     collected (i.e., the actual ensemble) and the position on the mean cross
     section. This function corrects the mean cross section velocities,
     assuming that mass flux is constant in the streamwise direction.

     Definition Sketch
     ==============================================
                        ..     ^
           Shiptrack  .        |
                +     .        |
                |      .       |
                +-----> .      |
                        .      |
                        .      |
                       .       |
                      .        |
                     .         | u2h2
               u1h1 +----------+
                   .           |
                   .           |
                    .          |<-------+
                     .         |        |
                      .        |        |
                        .      |        +
                         .     |       MCS
                          .    |
                            .  v

     Under the assumption that mass flux (i.e., unit discharge) is constant in
     the streamwise direction:
           qs1 = u1h1 = qs2 = u2h2
     is independent of s, the streamwise direction

     Therefore, we can correct u2 by
                   h1
           u2 = u1----
                   h2

     For a further discussion, see:
           Hoitink, A. J. F., F. A. Buschmann, and B. Vermeulen. 2009.
               Continuous measurements of discharge from a horizontal acoustic
               Doppler current profiler in a tidal river. WRR 45, W11406,
               doi:10.1029/2009WR007791.

     There are several reasons to be cautious in choosing to apply this
     correction:
     1.It makes the assumption that mass flux is constant in the
       streamwise direction. This is probably true for well selected
       cross-sections. But it could break down in cases where flow is strongly
       convergent or divergent, or if the deviation of the shiptrack is too
       great.
     2.Streamwise flow in situations of changing depths will accel./decell.
       according to the Bernoulli principle- and this is the reason for the
       correction. However, because of the beam spread, it is likely that
       in these situations where streamwise depth is varying, the the depth
       will be incorrectly estimated, and will lead to high or low biased
       velocities (depending on whether depth is deepening/shoaling).
     3.Essentially, this amounts to correction of the data by estimating the
       depth at the mean cross section. If the coverage of the shiptracks are
       sparse for a certain region of the mean cross section, the quality of
       the depth estimation will be significantly affected, leading to a bias
       in the correction. A situation where shiptracks may not perfectly
       represent the bed is near the bank, with flow seperation. Here, ship
       tracks often drift off the best fit line computed by VMT, since the
       majority of the data are better clustered there. In this case, the
       correction will become biased by the estimation of the depth at that
       location.

     Frank L. Engel, USGS

     SEE ALSO: interp1, nanmean, smooth

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- VMT_GridData2MeanXS This routine generates a uniformly spaced grid for the mean cross section and

## SOURCE CODE

    0001 function A = VMT_unitQcont(A,V,z)
    0002 % VMT_UNITQCONT applies unit discharge continuity correction MCS velocities
    0003 % This routine computes the correction, then applies it to the North and
    0004 % East velocity components of each ensemble and bin for the ADCP data. This
    0005 % routine runs BEFORE the data are gridded to the mean cross section, so
    0006 % if used, the correction applies to all subsequent computations (e.g., ZSD
    0007 % and Rozovskii decomposition).
    0008 %
    0009 % DISCUSSION
    0010 % In theory, the change in elevation in the streamwise direction
    0011 % should be zero (since we are creating a plane in the lateral direction).
    0012 % However, in practice, this may not be the case. Hopefully, the user has
    0013 % chosen an appropriate cross section location where the bed is not varying
    0014 % dramatically in the streamwise direction, but even still, some error may
    0015 % exist. One way to approach this potential error is to account for any
    0016 % changes in the mass flux between the sample location where data were
    0017 % collected (i.e., the actual ensemble) and the position on the mean cross
    0018 % section. This function corrects the mean cross section velocities,
    0019 % assuming that mass flux is constant in the streamwise direction.
    0020 %
    0021 % Definition Sketch
    0022 % ==============================================
    0023 %                    ..     ^
    0024 %       Shiptrack  .        |
    0025 %            +     .        |
    0026 %            |      .       |
    0027 %            +-----> .      |
    0028 %                    .      |
    0029 %                    .      |
    0030 %                   .       |
    0031 %                  .        |
    0032 %                 .         | u2h2
    0033 %           u1h1 +----------+
    0034 %               .           |
    0035 %               .           |
    0036 %                .          |<-------+
    0037 %                 .         |        |
    0038 %                  .        |        |
    0039 %                    .      |        +
    0040 %                     .     |       MCS
    0041 %                      .    |
    0042 %                        .  v
    0043 %
    0044 % Under the assumption that mass flux (i.e., unit discharge) is constant in
    0045 % the streamwise direction:
    0046 %       qs1 = u1h1 = qs2 = u2h2
    0047 % is independent of s, the streamwise direction
    0048 %
    0049 % Therefore, we can correct u2 by
    0050 %               h1
    0051 %       u2 = u1----
    0052 %               h2
    0053 %
    0054 % For a further discussion, see:
    0055 %       Hoitink, A. J. F., F. A. Buschmann, and B. Vermeulen. 2009.
    0056 %           Continuous measurements of discharge from a horizontal acoustic
    0057 %           Doppler current profiler in a tidal river. WRR 45, W11406,
    0058 %           doi:10.1029/2009WR007791.
    0059 %
    0060 % There are several reasons to be cautious in choosing to apply this
    0061 % correction:
    0062 % 1.It makes the assumption that mass flux is constant in the
    0063 %   streamwise direction. This is probably true for well selected
    0064 %   cross-sections. But it could break down in cases where flow is strongly
    0065 %   convergent or divergent, or if the deviation of the shiptrack is too
    0066 %   great.
    0067 % 2.Streamwise flow in situations of changing depths will accel./decell.
    0068 %   according to the Bernoulli principle- and this is the reason for the
    0069 %   correction. However, because of the beam spread, it is likely that
    0070 %   in these situations where streamwise depth is varying, the the depth
    0071 %   will be incorrectly estimated, and will lead to high or low biased
    0072 %   velocities (depending on whether depth is deepening/shoaling).
    0073 % 3.Essentially, this amounts to correction of the data by estimating the
    0074 %   depth at the mean cross section. If the coverage of the shiptracks are
    0075 %   sparse for a certain region of the mean cross section, the quality of
    0076 %   the depth estimation will be significantly affected, leading to a bias
    0077 %   in the correction. A situation where shiptracks may not perfectly
    0078 %   represent the bed is near the bank, with flow seperation. Here, ship
    0079 %   tracks often drift off the best fit line computed by VMT, since the
    0080 %   majority of the data are better clustered there. In this case, the
    0081 %   correction will become biased by the estimation of the depth at that
    0082 %   location.
    0083 %
    0084 % Frank L. Engel, USGS
    0085 %
    0086 % SEE ALSO: interp1, nanmean, smooth
    0087
    0088 % disp('Applying streamwise unit discharge correction...')
    0089 % Create matrix of all depths, with their stations
    0090 h2_interp = [];
    0091 for zi = 1:z
    0092     A(zi).Comp.h1 = nanmean(A(zi).Nav.depth,2)';
    0093     h2_interp = [h2_interp; A(zi).Comp.dl A(zi).Comp.h1'];
    0094 end
    0095
    0096 % Sort it by station
    0097 [~,idx] = sort(h2_interp(:,1));
    0098 h2_interp = h2_interp(idx,:);
    0099
    0100 % Remove duplicate stations (to avoid issues during interp1)
    0101 uniqs = find(diff(h2_interp(:,1)) ~= 0);
    0102 h2_interp = h2_interp([1 (uniqs+1)'],:);
    0103
    0104 % Remove nans in depths to aviod issues during smoothing
    0105 h2_interp = h2_interp(~isnan(h2_interp(:,2)),:);
    0106
    0107 % Smooth the bed, to get a representation of the MCS bed Estimate a
    0108 % smoothing span which will approximate that created by averaging the bed
    0109 % at each grid node.
    0110 span = 5*ceil(A(1).hgns/(V.meddens+V.stddens));
    0111 %h2_interp(:,3) = smooth(h2_interp(:,2),span);
    0112 %Replaced smooth as it did not handle edges well (PRJ, 10-19-12)
    0113 [h2_interp(:,3),~] = nanmoving_average(h2_interp(:,2),span/2);
    0114
    0115 if 0 %Plot to check
    0116     figure(8); clf
    0117     plot(h2_interp(:,1),h2_interp(:,2),'r.'); hold on
    0118     plot(h2_interp(:,1),h2_interp(:,3),'k-');
    0119     set(gca,'Ydir','reverse')
    0120 end
    0121
    0122 % Do the computations
    0123 PD = [];  %Percent diff
    0124 for zi = 1:z
    0125     % Ensure that the stationing is in the right direction, and then write
    0126     % h1 to the A.Comp struct
    0127     h1_sort = [A(zi).Comp.dl nanmean(A(zi).Nav.depth,2)];
    0128     [~,h1_idx] = sort(h1_sort(:,1));
    0129     h1_sort = h1_sort(h1_idx,:);
    0130     A(zi).Comp.h1 = h1_sort(:,2)';
    0131
    0132     % Interpolate the bed at the MCS for each ensemble in the current
    0133     % transect
    0134     A(zi).Comp.h2 = interp1(...
    0135         h2_interp(:,1),...
    0136         h2_interp(:,3),...
    0137         A(zi).Comp.itDist(1,:));
    0138
    0139     % To compute the correction, we have to get the intrinsic coordinates
    0140     % of the velocity vectors (neglecting vertical component). Technically,
    0141     % we are not correcting for streamwise changes, assuming that
    0142     % streamwise is taken exactly perpendicular to the MCS
    0143     A(zi).Comp.psi = V.phi-A(zi).Clean.vDir;
    0144     A(zi).Comp.u1 = cosd(A(zi).Comp.psi).*A(zi).Clean.vMag;
    0145     A(zi).Comp.v1 = sind(A(zi).Comp.psi).*A(zi).Clean.vMag;
    0146     A(zi).Comp.hratio = repmat(A(zi).Comp.h1./A(zi).Comp.h2,size(A(zi).Comp.itDist,1),1);
    0147
    0148     % Now we can compute the corrected velocities
    0149     A(zi).Comp.u2 = A(zi).Comp.u1.*A(zi).Comp.hratio;
    0150
    0151     % Save the percent differnce for statisical analysis
    0152     PD = [PD; (A(zi).Comp.hratio(1,:)' - 1)*100];
    0153     %A(zi).Comp.w1 = A(zi).Wat.vVert;
    0154
    0155     if 0 % Plot to check
    0156         figure(9); clf
    0157         plot(A(zi).Comp.itDist(1,:),A(zi).Comp.h1,'r.-'); hold on
    0158         plot(A(zi).Comp.itDist(1,:),A(zi).Comp.h2,'k.-');
    0159         set(gca,'Ydir','reverse')
    0160         pause
    0161     end
    0162
    0163 end
    0164
    0165 % Report the percent difference stats & plot profiles (for testing)
    0166 if 0
    0167     PDmed = nanmedian(PD);
    0168     PDmax = nanmax(PD);
    0169     PDmin = nanmin(PD);
    0170     scrnwrt = {['Median Percent Difference (%) = ' num2str(PDmed)];...
    0171                ['Max Percent Difference (%) = ' num2str(PDmax)];...
    0172                ['Min Percent Difference (%) = ' num2str(PDmin)]};
    0173     figure(10); clf
    0174     for zi = 1:z
    0175         plot(A(zi).Comp.itDist(1,:),(A(zi).Comp.hratio(1,:) - 1)*100,'k-')
    0176         hold on;
    0177         xlabel('Distance from Left Bank, in meters')
    0178         ylabel('Percent Difference in Velocity')
    0179     end
    0180     grid on
    0181     ylim([min([min(PD) -50]) max([max(PD) 50])])
    0182     text(0.1*mean(A(1).Comp.itDist(1,:)),25,scrnwrt)
    0183     set(gca,'YMinorGrid','on')
    0184 end
    0185
    0186 % Rotate the new corrected u,v back into the East,North frame of reference.
    0187 % Note that to recover the components, we need to rotate to -V.theta
    0188 % degrees.
    0189
    0190 % Rotation matrix
    0191 R = [cos(-geo2arideg(V.theta)*pi/180) sin(-geo2arideg(V.theta)*pi/180);...
    0192     -sin(-geo2arideg(V.theta)*pi/180) cos(-geo2arideg(V.theta)*pi/180)]';
    0193
    0194 % Overwrite the original East,North components from the original ASC inputs
    0195 for zi = 1:z
    0196     % The vectors to be rotated
    0197     N = [A(zi).Comp.u2(:) A(zi).Comp.v1(:)];
    0198     % Do the rotation
    0199     T = N*R;
    0200
    0201 %     if V.phi >= 0 && V.phi <= 180
    0202 %         East = -reshape(T(:,2),size(A(zi).Clean.vEast));
    0203 %     elseif V.phi > 180 && V.phi < 360
    0204 %         East = reshape(T(:,2),size(A(zi).Clean.vEast));
    0205 %     end
    0206
    0207     if V.phi >90 && V.phi < 270
    0208         North = -reshape(T(:,1),size(A(zi).Clean.vNorth));
    0209         East  = reshape(T(:,2),size(A(zi).Clean.vEast));
    0210     elseif V.phi >= 270 || V.phi <= 90
    0211         North = reshape(T(:,1),size(A(zi).Clean.vNorth));
    0212         East = -reshape(T(:,2),size(A(zi).Clean.vEast));
    0213     end
    0214     % Write the result back to the A struct
    0215     A(zi).Clean.vEast = East;
    0216     A(zi).Clean.vNorth = North;
    0217     clear N T East North
    0218 end