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
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