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:

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

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005