Not used by current version This routine acts as a driver program to process multiple transects at a single cross-section for velocity mapping.
0001 function [A,V] = VMT_ProcessTransectsV3_new(z,A,setends) 0002 % Not used by current version 0003 %This routine acts as a driver program to process multiple transects at a 0004 %single cross-section for velocity mapping. 0005 0006 %V2 adds the cpability to set the endpoints of the mean cross section 0007 0008 %V3 adds the Rozovskii computations for secondary flow. 8/31/09 0009 0010 %Among other things, it: 0011 0012 % Determines the best fit mean cross-section line from multiple transects 0013 % Map ensembles to mean c-s line 0014 % Determine uniform mean c-s grid for vector interpolating 0015 % Determine location of mapped ensemble points for interpolating 0016 % Interpolate individual transects onto uniform mean c-s grid 0017 % Average mapped mean cross-sections from individual transects together 0018 % Rotate velocities into u, v, and w components 0019 0020 0021 %(adapted from code by J. Czuba) 0022 0023 %P.R. Jackson, USGS, 12-9-08 0024 0025 disp('Processing Data...') 0026 % warning off 0027 %% Map ensembles to mean cross-section 0028 0029 [A,V] = VMT_MapEns2MeanXSV2(z,A,setends); 0030 0031 %% Grid the measured data along the mean cross-section 0032 %[A,V] = VMT_GridData2MeanXS(z,A,V); 0033 [A,V] = VMT_GridData2MeanXS(z,A,V); 0034 0035 %% Computes the mean data for the mean cross-section 0036 %[A,V] = VMT_CompMeanXS(z,A,V); 0037 [A,V] = VMT_CompMeanXS_old(z,A,V); 0038 0039 %% Decompose the velocities into u, v, and w components 0040 [A,V] = VMT_CompMeanXS_UVW(z,A,V); 0041 0042 %% Decompose the velocities into primary and secondary components 0043 [A,V] = VMT_CompMeanXS_PriSec(z,A,V); 0044 0045 %% Perform the Rozovskii computations 0046 V = VMT_RozovskiiV2(V,A); 0047 0048 %% 0049 %figure(4); clf 0050 %plot3(V.mcsX,V.mcsY,V.mcsDepth(1)) 0051 disp('Processing Completed') 0052 0053 %% Notes: 0054 0055 %1. I removed scripts at the end of the original code that computes 0056 %standard deviations (12-9-08) 0057 0058 % else 0059 % 0060 % disp('Processing Data...') 0061 % %% Map ensembles to mean cross-section 0062 % hf = VMT_CreatePlotDisplay('shiptracks'); 0063 % 0064 % [A,V] = VMT_MapEns2MeanXSV2_PDF(hf,z,A,setends); 0065 % 0066 % %% Grid the measured data along the mean cross-section 0067 % %[A,V] = VMT_GridData2MeanXS(z,A,V); 0068 % [A,V] = VMT_GridData2MeanXS_PDF(hf,z,A,V); 0069 % 0070 % return 0071 % %% Computes the mean data for the mean cross-section 0072 % %[A,V] = VMT_CompMeanXS(z,A,V); 0073 % [A,V] = VMT_CompMeanXS_old(z,A,V); 0074 % 0075 % %% Decompose the velocities into u, v, and w components 0076 % [A,V] = VMT_CompMeanXS_UVW(z,A,V); 0077 % 0078 % %% Decompose the velocities into primary and secondary components 0079 % [A,V] = VMT_CompMeanXS_PriSec(z,A,V); 0080 % 0081 % %% Perform the Rozovskii computations 0082 % V = VMT_RozovskiiV2(V,A); 0083 % 0084 % %% 0085 % %figure(4); clf 0086 % %plot3(V.mcsX,V.mcsY,V.mcsDepth(1)) 0087 % disp('Processing Completed') 0088 % 0089 % %% Notes: 0090 % 0091 % %1. I removed scripts at the end of the original code that computes 0092 % %standard deviations (12-9-08) 0093 % 0094 % end 0095 0096 %========================================================================== 0097 function [A,V] = VMT_MapEns2MeanXSV2(z,A,setends) 0098 0099 %This routine fits multiple transects at a single location with a single 0100 %line and maps individual ensembles to this line. Inputs are number of files (z) and data matrix (Z)(see ReadFiles.m). 0101 %Output is the updated data matrix with new mapped variables. 0102 0103 %V2 adds the capability to set the endpoints of the mean cross section 0104 0105 %(adapted from code by J. Czuba) 0106 0107 %P.R. Jackson, USGS, 12-9-08 0108 0109 %% Determine the best fit mean cross-section line from multiple transects 0110 % initialize vectors for concatenation 0111 0112 x = []; 0113 y = []; 0114 % figure(1); clf 0115 % set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1]) 0116 mfd = zeros(1,z); 0117 for zi = 1 : z 0118 0119 % concatenate coords into a single column vector for regression 0120 x = cat(1,x,A(zi).Comp.xUTM); 0121 y = cat(1,y,A(zi).Comp.yUTM); 0122 0123 % figure(1); hold on 0124 % %plot(A(zi).Comp.xUTM,A(zi).Comp.yUTM,'r'); hold on 0125 % plot(A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,'b'); hold on 0126 0127 %Find the mean flow direction for each transect 0128 mfd(zi) = nanmean(nanmean(A(zi).Clean.vDir,1)); 0129 end 0130 V.mfd = nanmean(mfd); % Mean flow direction for all the transects 0131 0132 if setends %Gets a user text file with fixed cross section end points 0133 0134 defaultpath = 'C:\'; 0135 defaultpath = pwd; 0136 endspath = []; 0137 if exist('\VMT\LastDir.mat') == 2 0138 load('VMT\LastDir.mat'); 0139 if exist(endspath) == 7 0140 [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',endspath); 0141 else 0142 [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',defaultpath); 0143 end 0144 else 0145 [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',defaultpath); 0146 end 0147 0148 if ischar(file) 0149 infile = [endspath file]; 0150 %[file,path] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File'); 0151 %infile = [path file]; 0152 disp('Loading Endpoint File...' ); 0153 disp(infile); 0154 data = dlmread(infile); 0155 x = data(:,1); 0156 y = data(:,2); 0157 figure(1); hold on 0158 plot(x,y,'go','MarkerSize',10); hold on 0159 end 0160 end 0161 0162 % find the equation of the best fit line 0163 xrng = max(x) - min(x); 0164 yrng = max(y) - min(y); 0165 0166 if xrng >= yrng %Fit based on coordinate with larger range of values (original fitting has issues with N-S lines because of repeated X values), PRJ 12-12-08 0167 [P,~] = polyfit(x,y,1); 0168 % figure(1); hold on; 0169 % plot(x,polyval(P,x),'g-') 0170 V.m = P(1); 0171 V.b = P(2); 0172 else 0173 [P,~] = polyfit(y,x,1); 0174 % figure(1); hold on; 0175 % plot(polyval(P,y),y,'g-') 0176 V.m = 1/P(1); %Reformat slope and intercept in terms of y= fn(x) rather than x = fn(y) 0177 V.b = -P(2)/P(1); 0178 end 0179 0180 %Former method commented out 0181 % whichstats = {'tstat','yhat'}; 0182 % stats = regstats(y,x,'linear', whichstats); 0183 % 0184 % % mean cross-section line slope and intercept 0185 % V.m = stats.tstat.beta(2); 0186 % V.b = stats.tstat.beta(1); 0187 0188 clear x y stats whichstats zi 0189 0190 %% Map ensembles to mean c-s line 0191 % Determine the point (mapped ensemble point) where the equation of the 0192 % mean cross-section line intercepts a line perpendicular to the mean 0193 % cross-section line passing through an ensemble from an individual 0194 % transect (see notes for equation derivation) 0195 0196 for zi = 1:z 0197 A(zi).Comp.xm = ((A(zi).Comp.xUTM-V.m.*V.b+V.m.*A(zi).Comp.yUTM)... 0198 ./(V.m.^2+1)); 0199 A(zi).Comp.ym = ((V.b+V.m.*A(zi).Comp.xUTM+V.m.^2.*A(zi).Comp.yUTM)... 0200 ./(V.m.^2+1)); 0201 end 0202 0203 %Plot data to check 0204 xensall = []; 0205 yensall = []; 0206 for zi = 1:z 0207 % plot(A(zi).Comp.xm,A(zi).Comp.ym,'b.') 0208 xensall = [xensall; A(zi).Comp.xm]; 0209 yensall = [yensall; A(zi).Comp.ym]; 0210 end 0211 % % plot(A(3).Comp.xm,A(3).Comp.ym,'xg') 0212 % % plot(A(4).Comp.xm,A(4).Comp.ym,'oy') 0213 % xlabel('UTM Easting (m)') 0214 % ylabel('UTM Northing (m)') 0215 % box on 0216 % grid on 0217 % %Plot a legend in Figure 1 0218 % %figure(1); hold on 0219 % %legend('Shoreline','GPS(corr)','GPS(raw)','Best Fit','Trans 1 0220 % %(mapped)','Other Trans (mapped)') 0221 0222 %Compute the median distance between mapped points 0223 Dmat = [xensall yensall]; 0224 if xrng > yrng 0225 Dmat = sortrows(Dmat,1); 0226 else 0227 Dmat = sortrows(Dmat,2); 0228 end 0229 dxall = diff(Dmat(:,1)); 0230 dyall = diff(Dmat(:,2)); 0231 densall = sqrt(dxall.^2 + dyall.^2); 0232 V.meddens = median(densall); 0233 V.stddens = std(densall); 0234 disp(['Median Spacing Between Mapped Ensembles = ' num2str(V.meddens) ' m']) 0235 disp(['Standard Deviation of Spacing Between Mapped Ensembles = ' num2str(V.stddens) ' m']) 0236 disp(['Recommended Grid Node Spacing > ' num2str(V.meddens + V.stddens) ' m']) 0237 0238 %Display in message box for compiled version 0239 msg_string = {['Median Spacing Between Mapped Ensembles = ' num2str(V.meddens) ' m'],... 0240 ['Standard Deviation of Spacing Between Mapped Ensembles = ' num2str(V.stddens) ' m'],... 0241 ['Recommended Grid Node Spacing > ' num2str(V.meddens + V.stddens) ' m']}; 0242 msgbox(msg_string,'VMT Grid Node Spacing','help','replace'); 0243 0244 %%Save the shorepath 0245 if setends 0246 if exist('LastDir.mat','dir') 0247 save('LastDir.mat','endspath','-append') 0248 else 0249 save('LastDir.mat','endspath') 0250 end 0251 end 0252 0253 %========================================================================== 0254 function [A,V] = VMT_GridData2MeanXS(z,A,V) 0255 0256 %This routine generates a uniformly spaced grid for the mean cross section and 0257 %maps (interpolates) individual transects to this grid. 0258 0259 %(adapted from code by J. Czuba) 0260 0261 %P.R. Jackson, USGS, 12-9-08 0262 0263 %% User Input 0264 0265 xgdspc = A(1).hgns; %Horizontal Grid node spacing in meters (vertical grid spacing is set by bins) 0266 % if 0 0267 % xgdspc = V.meddens + V.stddens; %Auto method should include 67% of the values 0268 % disp(['X Grid Node Auto Spacing = ' num2str(xgdspc) ' m']) 0269 % end 0270 0271 %% Determine uniform mean c-s grid for vector interpolating 0272 % Determine the end points of the mean cross-section line 0273 % Initialize variable with mid range value 0274 V.xe = mean(A(1).Comp.xm); 0275 V.ys = mean(A(1).Comp.ym); 0276 V.xw = mean(A(1).Comp.xm); 0277 V.yn = mean(A(1).Comp.ym); 0278 0279 for zi = 1:z 0280 V.xe = max(max(A(zi).Comp.xm),V.xe); 0281 V.ys = min(min(A(zi).Comp.ym),V.ys); 0282 V.xw = min(min(A(zi).Comp.xm),V.xw); 0283 V.yn = max(max(A(zi).Comp.ym),V.yn); 0284 end 0285 0286 % Determine the distance between the mean cross-section endpoints 0287 V.dx = V.xe-V.xw; 0288 V.dy = V.yn-V.ys; 0289 0290 V.dl = sqrt(V.dx^2+V.dy^2); 0291 0292 % Determine mean cross-section velocity vector grid 0293 V.mcsDist = linspace(0,V.dl,floor(V.dl/xgdspc)); %%linspace(0,V.dl,V.dl); Changed to make it user selectable (PRJ, 12-12-08) 0294 V.mcsDepth = A(1).Wat.binDepth(:,1); 0295 [V.mcsDist V.mcsDepth] = meshgrid(V.mcsDist,V.mcsDepth); 0296 0297 % Determine the angle the mean cross-section makes with the 0298 % x-axis (longitude) 0299 % Plot mean cross-section line 0300 if V.m >= 0 0301 V.theta = atand(V.dy./V.dx); 0302 0303 % figure(1); hold on 0304 % plot([V.xe V.xw],[V.yn V.ys],'ks'); hold on 0305 0306 V.mcsX = V.xe - V.mcsDist(1,:).*cosd(V.theta); % 0307 V.mcsY = V.yn - V.mcsDist(1,:).*sind(V.theta); 0308 0309 % if V.mfd >= 270 | V.mfd < 90 %Flow to the north %This code was an attempt to auto detect left bank--did'nt work so removed. 0310 % V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta); % 0311 % V.mcsY = V.ys+V.mcsDist(1,:).*sind(V.theta); 0312 % 0313 % elseif V.mfd >= 90 & V.mfd < 270 %Flow to the south 0314 % V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta); % 0315 % V.mcsY = V.yn-V.mcsDist(1,:).*sind(V.theta); 0316 % end% 0317 0318 % plot(V.mcsX,V.mcsY,'k+'); hold on 0319 % plot(V.mcsX(1),V.mcsY(1),'y*'); hold on 0320 0321 elseif V.m < 0 0322 V.theta = -atand(V.dy./V.dx); 0323 %V.theta = atand(V.dy./V.dx); %Changed 9-28-10, PRJ (theta needs to be 0324 %negative--changed back to original) 0325 0326 % figure(1); hold on 0327 % plot([V.xe V.xw],[V.ys V.yn],'ks'); hold on 0328 0329 V.mcsX = V.xe - V.mcsDist(1,:).*cosd(V.theta); 0330 V.mcsY = V.ys - V.mcsDist(1,:).*sind(V.theta); 0331 %V.mcsY = V.ys+V.mcsDist(1,:).*sind(V.theta); %Changed 9-28-10, PRJ 0332 0333 % if V.mfd >= 270 | V.mfd < 90 %Flow to the north 0334 % V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta); % 0335 % V.mcsY = V.yn+V.mcsDist(1,:).*sind(V.theta); 0336 % 0337 % elseif V.mfd >= 90 & V.mfd < 270 %Flow to the south 0338 % V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta); 0339 % V.mcsY = V.ys-V.mcsDist(1,:).*sind(V.theta); 0340 % end% 0341 0342 % plot(V.mcsX,V.mcsY,'k+'); hold on 0343 % plot(V.mcsX(1),V.mcsY(1),'y*'); hold on 0344 end 0345 0346 V.mcsX = meshgrid(V.mcsX,V.mcsDepth(:,1)); 0347 V.mcsY = meshgrid(V.mcsY,V.mcsDepth(:,1)); 0348 % figure(1); set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1]) 0349 clear zi 0350 0351 % Format the ticks for UTM and allow zooming and panning 0352 % figure(1); 0353 % ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM 0354 % hdlzm_fig1 = zoom; 0355 % set(hdlzm_fig1,'ActionPostCallback',@mypostcallback_zoom); 0356 % set(hdlzm_fig1,'Enable','on'); 0357 % hdlpn_fig1 = pan; 0358 % set(hdlpn_fig1,'ActionPostCallback',@mypostcallback_pan); 0359 % set(hdlpn_fig1,'Enable','on'); 0360 0361 0362 %% Determine location of mapped ensemble points for interpolating 0363 % For all transects 0364 0365 %A = VMT_DxDyfromLB(z,A,V); %Computes dx and dy 0366 0367 for zi = 1:z 0368 % Determine if the mean cross-section line trends NW-SE or SW-NE 0369 % Determine the distance in radians from the left bank mean 0370 % cross-section point to the mapped ensemble point for an individual 0371 % transect 0372 A(zi).Comp.dx = abs(V.xe-A(zi).Comp.xm); %This assumes the easternmost bank is the left bank--changed PRJ 1-21-09 (now uses VMT_DxDyfromLB--not working 2/1/09) 0373 0374 if V.m > 0 0375 A(zi).Comp.dy = abs(V.yn-A(zi).Comp.ym); 0376 elseif V.m < 0 0377 A(zi).Comp.dy = abs(V.ys-A(zi).Comp.ym); 0378 end 0379 0380 % Determine the distance in meters from the left bank mean 0381 % cross-section point to the mapped ensemble point for an individual 0382 % transect 0383 A(zi).Comp.dl = sqrt(A(zi).Comp.dx.^2+A(zi).Comp.dy.^2); 0384 0385 % % Sort vectors by dl 0386 % A(zi).Comp.dlsort = sort(A(zi).Comp.dl,'ascend'); 0387 % 0388 % % Map indices 0389 % for i = 1:A(zi).Sup.noe 0390 % for k = 1 : A(zi).Sup.noe 0391 % if A(zi).Comp.dlsort(i,1) == A(zi).Comp.dl(k,1) 0392 % A(zi).Comp.vecmap(i,1) = k; 0393 % end 0394 % end 0395 % end 0396 0397 % Sort vectors by dl 0398 [A(zi).Comp.dlsort,A(zi).Comp.vecmap] = sort(A(zi).Comp.dl,'ascend'); 0399 0400 % GPS position fix 0401 % if distances from the left bank are the same for two ensembles the 0402 % the position of the right most ensemble is interpolated from adjacent 0403 % ensembles 0404 % check for repeat values of distance 0405 sbt(:,1) = diff(A(zi).Comp.dlsort); 0406 chk(1,1) = 1; 0407 chk(2:A(zi).Sup.noe,1) = sbt(1:end,1); 0408 0409 % identify repeat values 0410 A(zi).Comp.sd = (chk==0) > 0; 0411 0412 % if repeat values exist interpolate distances from adjacent ensembles 0413 if sum(A(zi).Comp.sd) > 0 0414 0415 % bracket repeat sections 0416 [I,~] = ind2sub(size(A(zi).Comp.sd),find(A(zi).Comp.sd==1)); 0417 df = diff(I); 0418 nbrk = sum(df>1)+1; 0419 [I2,~] = ind2sub(size(df),find(df>1)); 0420 0421 bg(1) = I(1)-1; 0422 for n = 2 : nbrk 0423 bg(n) = I(I2(n-1)+1)-1; 0424 end 0425 0426 for n = 1:nbrk -1 0427 ed(n) = I(I2(n))+1; 0428 end 0429 ed(nbrk) = I(end)+1; 0430 0431 % interpolate repeat values 0432 A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort; 0433 0434 for i = 1 : nbrk 0435 for j = bg(i)+1 : ed(i)-1 0436 % interpolate 0437 if bg(i) > 0 && ed(i) < length(A(zi).Nav.lat_deg) 0438 0439 den=(ed(i)-bg(i)); 0440 num2=j-bg(i); 0441 num1=ed(i)-j; 0442 0443 A(zi).Comp.dlsortgpsfix(j,1)=... 0444 (num1/den).*A(zi).Comp.dlsort(bg(i))+... 0445 (num2/den).*A(zi).Comp.dlsort(ed(i)); 0446 0447 end 0448 0449 % extrapolate end 0450 if ed(i) > length(A(zi).Nav.lat_deg) 0451 0452 numex = ed(i) - length(A(zi).Nav.lat_deg); 0453 0454 A(zi).Comp.dlsortgpsfix(j,1)=numex.*... 0455 (A(zi).Comp.dlsort(bg(i))-... 0456 A(zi).Comp.dlsort(bg(i)-1))+... 0457 A(zi).Comp.dlsort(bg(i)); 0458 0459 end 0460 end 0461 end 0462 else 0463 A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort; 0464 end % IF sum(A(zi).Comp.sd) > 0 0465 0466 % Determine velocity vector grid for individual transects 0467 [A(zi).Comp.itDist A(zi).Comp.itDepth] = ... 0468 meshgrid(A(zi).Comp.dlsortgpsfix,A(zi).Wat.binDepth(:,1)); 0469 0470 clear I I2 J J2 bg chk df ed i j nbrk sbt xUTM yUTM n zi... 0471 den num2 num1 numex 0472 0473 end 0474 0475 clear zi i k check 0476 0477 %% Interpolate individual transects onto uniform mean c-s grid 0478 % Fill in uniform grid based on individual transects mapped onto the mean 0479 % cross-section by interpolating between adjacent points 0480 0481 %ZI = interp2(X,Y,Z,XI,YI) 0482 for zi = 1:z 0483 A(zi).Comp.mcsBack = ... 0484 interp2(A(zi).Comp.itDist, ... 0485 A(zi).Comp.itDepth, ... 0486 A(zi).Clean.bs(:,A(zi).Comp.vecmap), ... 0487 V.mcsDist, ... 0488 V.mcsDepth); 0489 A(zi).Comp.mcsBack(A(zi).Comp.mcsBack>=255) = NaN; 0490 0491 %A(zi).Comp.mcsDir = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ... 0492 %A(zi).Clean.vDir(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth); %This interpolation scheme has issues when interpolating in a flow due north (0,360 interpolate to 180) 0493 0494 % For direction, must convert degrees to radians, take the sin of the 0495 % radians, and then interpolate. Following interpolation, convert 0496 % radians back to degrees. (PRJ, 9-28-10) ALSO BAD NEAR 180 0497 %A(zi).Comp.mcsDir = 180/pi*(interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ... 0498 %sin(pi/180*(A(zi).Clean.vDir(:,A(zi).Comp.vecmap))), V.mcsDist, V.mcsDepth)); 0499 0500 % A(zi).Comp.mcsMag = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ... 0501 % A(zi).Clean.vMag(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth); 0502 % (Recomputed from north and east components (PRJ, 3-21-11) 0503 0504 0505 A(zi).Comp.mcsEast = ... 0506 interp2(A(zi).Comp.itDist, ... 0507 A(zi).Comp.itDepth, ... 0508 A(zi).Clean.vEast(:,A(zi).Comp.vecmap), ... 0509 V.mcsDist, ... 0510 V.mcsDepth); 0511 A(zi).Comp.mcsNorth = ... 0512 interp2(A(zi).Comp.itDist, ... 0513 A(zi).Comp.itDepth, ... 0514 A(zi).Clean.vNorth(:,A(zi).Comp.vecmap), ... 0515 V.mcsDist, ... 0516 V.mcsDepth); 0517 A(zi).Comp.mcsVert = ... 0518 interp2(A(zi).Comp.itDist, ... 0519 A(zi).Comp.itDepth, ... 0520 A(zi).Clean.vVert(:,A(zi).Comp.vecmap), ... 0521 V.mcsDist, ... 0522 V.mcsDepth); 0523 0524 %Compute magnitude 0525 A(zi).Comp.mcsMag = sqrt(A(zi).Comp.mcsEast.^2 + A(zi).Comp.mcsNorth.^2); 0526 0527 %For direction, compute from the velocity components 0528 A(zi).Comp.mcsDir = 90 - (atan2(A(zi).Comp.mcsNorth, A(zi).Comp.mcsEast))*180/pi; %Compute the atan from the velocity componentes, convert to radians, and rotate to north axis 0529 qindx = find(A(zi).Comp.mcsDir < 0); 0530 if ~isempty(qindx) 0531 A(zi).Comp.mcsDir(qindx) = A(zi).Comp.mcsDir(qindx) + 360; %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis 0532 end 0533 0534 A(zi).Comp.mcsBed = ... 0535 interp1(A(zi).Comp.itDist(1,:), ... 0536 nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2), ... 0537 V.mcsDist(1,:)); 0538 end 0539 0540 clear zi 0541 0542 %% Embedded functions 0543 function mypostcallback_zoom(obj,evd) 0544 ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM (when zooming) 0545 0546 function mypostcallback_pan(obj,evd) 0547 ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM (when panning) 0548 0549 %========================================================================== 0550 function [A,V] = VMT_CompMeanXS_old(z,A,V) 0551 % Compute the mean cross section data from individual transects that have 0552 % been previously mapped to a common grid. 0553 0554 %(adapted from code by J. Czuba) 0555 0556 %P.R. Jackson, USGS, 12-9-08 0557 0558 %% Average mapped mean cross-sections from individual transects together 0559 % Assign mapped uniform grid vectors to the same array for averaging 0560 Back = zeros([size(A(1).Comp.mcsBack) z]); 0561 Dir = zeros([size(A(1).Comp.mcsDir) z]); 0562 Mag = zeros([size(A(1).Comp.mcsMag) z]); 0563 East = zeros([size(A(1).Comp.mcsEast) z]); 0564 North = zeros([size(A(1).Comp.mcsNorth) z]); 0565 Vert = zeros([size(A(1).Comp.mcsVert) z]); 0566 Bed = zeros([size(A(1).Comp.mcsBed) z]); 0567 for zi = 1 : z 0568 Back(:,:,zi) = A(zi).Comp.mcsBack; 0569 Dir(:,:,zi) = A(zi).Comp.mcsDir; 0570 Mag(:,:,zi) = A(zi).Comp.mcsMag; 0571 East(:,:,zi) = A(zi).Comp.mcsEast; 0572 North(:,:,zi) = A(zi).Comp.mcsNorth; 0573 Vert(:,:,zi) = A(zi).Comp.mcsVert; 0574 Bed(:,:,zi) = A(zi).Comp.mcsBed; 0575 end 0576 0577 % numavg = nansum(~isnan(Mag),3); 0578 % numavg(numavg==0) = NaN; 0579 % enscnt = nanmean(numavg,1); 0580 % [I,J] = ind2sub(size(enscnt),find(enscnt>=1)); %Changed to >= 1 PRJ 12-10-08 (uses data even if only one measurement) 0581 0582 Backone = Back; 0583 Magone = Mag; 0584 Vertone = Vert; 0585 Bedone = Bed; 0586 0587 Backone(~isnan(Back)) = 1; 0588 Magone(~isnan(Mag)) = 1; 0589 Vertone(~isnan(Vert)) = 1; 0590 Bedone(~isnan(Bed)) = 1; 0591 0592 V.countBack = nansum(Backone,3); 0593 V.countMag = nansum(Magone,3); 0594 V.countVert = nansum(Vertone,3); 0595 V.countBed = nansum(Bedone,3); 0596 0597 V.countBack(V.countBack==0) = NaN; 0598 V.countMag(V.countMag==0) = NaN; 0599 V.countVert(V.countVert==0) = NaN; 0600 V.countBed(V.countBed==0) = NaN; 0601 0602 % Average mapped mean cross-sections from individual transects together 0603 V.mcsBack = nanmean(Back,3); 0604 % V.mcsDir = nanmean(Dir,3); % Will not average correctly in all cases due to 0-360 0605 %wrapping (PRJ, 9-29-10) 0606 % V.mcsMag = nanmean(Mag,3); %Mag recomputed from north, east, up components(PRJ, 3-21-11) 0607 V.mcsEast = nanmean(East,3); 0608 V.mcsNorth = nanmean(North,3); 0609 V.mcsVert = nanmean(Vert,3); 0610 0611 % Average Magnitude 0612 V.mcsMag = sqrt(V.mcsEast.^2 + V.mcsNorth.^2 + V.mcsVert.^2); 0613 0614 % Average the flow direction 0615 V.mcsDir = 90 - (atan2(V.mcsNorth, V.mcsEast))*180/pi; %Compute the atan from the velocity componentes, convert to radians, and rotate to north axis 0616 qindx = find(V.mcsDir < 0); 0617 if ~isempty(qindx) 0618 V.mcsDir(qindx) = V.mcsDir(qindx) + 360; %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis 0619 end 0620 0621 V.mcsBed = nanmean(Bed,3); 0622 0623 % Compute the Bed Elevation in meters 0624 disp(['Assigned Water Surface Elevation (WSE; in meters) = ' num2str(A(1).wse)]) 0625 V.mcsBedElev = A(1).wse - V.mcsBed; 0626 0627 %========================================================================== 0628 function [A,V] = VMT_CompMeanXS_UVW(z,A,V) 0629 0630 %This routine computes the mean cross section velocity components (UVW) 0631 %from individual transects that have been previously mapped to a common grid and averaged. 0632 0633 %(adapted from code by J. Czuba) 0634 0635 %P.R. Jackson, USGS, 12-9-08 0636 0637 0638 %% Rotate velocities into u, v, and w components 0639 % Determine the direction of streamwise velocity (u) 0640 V.phi = 180 - V.theta; %Taken as perpendicular to the mean XS 0641 0642 % Determine the deviation of a vector from streamwise velocity 0643 V.psi = V.phi - V.mcsDir; 0644 0645 % Determine streamwise (u), transverse (v), and vertical (w) velocities 0646 V.u = cosd(V.psi).*V.mcsMag; 0647 V.v = sind(V.psi).*V.mcsMag; 0648 V.w = V.mcsVert; 0649 0650 for zi = 1:z 0651 A(zi).Comp.u = cosd(V.psi).*A(zi).Comp.mcsMag; 0652 A(zi).Comp.v = sind(V.psi).*A(zi).Comp.mcsMag; 0653 A(zi).Comp.w = A(zi).Comp.mcsVert; 0654 end 0655 0656 %========================================================================== 0657 function [A,V] = VMT_CompMeanXS_PriSec(z,A,V) 0658 0659 %This routine computes the mean cross section velocity components (Primary and secondary) 0660 %from individual transects that have been previously mapped to a common grid and averaged. 0661 %The Primary velocity is defined as the component of the flow in the direction of the discharge 0662 %(i.e. rotated from the streamwise direction so the secrondary discharge is 0663 %zero). 0664 0665 % This is referred to as the "zero net cross-stream discharge definition" 0666 % (see Lane et al. 2000, Hydrological Processes 14, 2047-2071) 0667 0668 %(adapted from code by J. Czuba) 0669 0670 %P.R. Jackson, USGS, 12-9-08 0671 0672 %% Rotate velocities into p and s components for the mean transect 0673 % calculate dy and dz for each meaurement point 0674 dy = mean(diff(V.mcsDist(1,:))); % m 0675 dz = mean(diff(V.mcsDepth(:,1))); % m 0676 dydz = dy.*dz; 0677 0678 % calculate the bit of discharge for each imaginary cell around the 0679 % velocity point 0680 qyi = V.v.*dydz; % cm*m^2/s 0681 qxi = V.u.*dydz; % cm*m^2/s 0682 % qyi = V.v.*dy.*dz;%cm*m^2/s 0683 % qxi = V.u.*dy.*dz;%cm*m^2/s 0684 0685 % sum the streamwise and transverse Q and calculate the angle of the 0686 % cross section 0687 V.Qy = nansum(nansum(qyi));%cm*m^2/s 0688 V.Qx = nansum(nansum(qxi));%cm*m^2/s 0689 % V.Qy = nansum(qyi(:)); % cm*m^2/s % PDF? 0690 % V.Qx = nansum(qxi(:)); % cm*m^2/s 0691 0692 V.alphasp = atand(V.Qy./V.Qx); 0693 V.phisp = V.phi - V.alphasp; 0694 0695 % Rotate the velocities so that Qy is effectively zero 0696 qpi = qxi.*cosd(V.alphasp) + qyi.*sind(V.alphasp); 0697 qsi = -qxi.*sind(V.alphasp) + qyi.*cosd(V.alphasp); 0698 % R = [ cosd(V.alphasp) sind(V.alphasp) % PDF? Depends on dimensions 0699 % -sind(V.alphasp) cosd(V.alphasp)]; 0700 0701 V.Qp = nansum(nansum(qpi));%cm*m^2/s 0702 V.Qs = nansum(nansum(qsi));%cm*m^2/s 0703 % V.Qp = nansum(qpi(:)));%cm*m^2/s % PDF? 0704 % V.Qs = nansum(qsi(:)));%cm*m^2/s 0705 disp(['Secondary Discharge after Rotation (ZSD definition; m^3/s) = ' num2str(V.Qs/100)]) 0706 0707 V.vp = qpi./dydz; % cm/s 0708 V.vs = qsi./dydz; % cm/s 0709 % V.vp = qpi./(dy.*dz);%cm/s 0710 % V.vs = qsi./(dy.*dz);%cm/s 0711 0712 %% Transform each individual transect 0713 0714 for zi = 1 : z 0715 % calculate the bit of discharge for each imaginary cell around the 0716 % velocity point 0717 A(zi).Comp.qyi = A(zi).Comp.v.*dydz; % cm*m^2/s 0718 A(zi).Comp.qxi = A(zi).Comp.u.*dydz; % cm*m^2/s 0719 % A(zi).Comp.qyi=A(zi).Comp.v.*dy.*dz;%cm*m^2/s 0720 % A(zi).Comp.qxi=A(zi).Comp.u.*dy.*dz;%cm*m^2/s 0721 0722 % rotate the velocities so that Qy is effcetively zero 0723 A(zi).Comp.qpi = A(zi).Comp.qxi.*cosd(V.alphasp) + A(zi).Comp.qyi.*sind(V.alphasp); 0724 A(zi).Comp.qsi = -A(zi).Comp.qxi.*sind(V.alphasp) + A(zi).Comp.qyi.*cosd(V.alphasp); 0725 0726 A(zi).Comp.Qp=nansum(nansum(A(zi).Comp.qpi));%cm*m^2/s 0727 A(zi).Comp.Qs=nansum(nansum(A(zi).Comp.qsi));%cm*m^2/s 0728 % A(zi).Comp.Qp = nansum(A(zi).Comp.qpi(:)); % cm*m^2/s % PDF? 0729 % A(zi).Comp.Qs = nansum(A(zi).Comp.qsi(:)); % cm*m^2/s 0730 0731 A(zi).Comp.vp = A(zi).Comp.qpi./dydz; % cm/s 0732 A(zi).Comp.vs =A(zi).Comp.qsi./dydz; % cm/s 0733 % A(zi).Comp.vp=A(zi).Comp.qpi./(dy.*dz);%cm/s 0734 % A(zi).Comp.vs=A(zi).Comp.qsi./(dy.*dz);%cm/s 0735 end 0736 0737 0738 %% Determine velocity deviations from the p direction 0739 0740 V.mcsDirDevp = V.phisp - V.mcsDir; 0741 0742 for zi = 1:z 0743 A(zi).Comp.mcsDirDevp = V.phisp - A(zi).Comp.mcsDir; 0744 end 0745 0746 %========================================================================== 0747 function [V] = VMT_RozovskiiV2(V,A) 0748 0749 % Computes the frame of reference transpositon as described in Kenworthy 0750 % and Rhoads (1998) ESPL using a Rozovskii type analysis. 0751 0752 % Notes: 0753 % -The depth averaging currently linearly interpolates to the bed, 0754 % however we may want some other approach such as log law, etc. 0755 % -I extrapolate the velocity at the near surface bin to the water 0756 % surface for the depth averaging (ie, BC at u(z=0) = u(z=bin1)) 0757 % -There are cases where the bin corresponding with the bed actually 0758 % contains flow data (i.e., not NaN or zero). For cases where the 0759 % blanking distance DOES exists, I have imposed a BC of U=0 at the bed, 0760 % -In cases where data goes all of the way to the bed, I have divided 0761 % the last bin's velocity by 2, essentially imposing a U=0 at the 0762 % boundary by extrapolating to the bottom of the bin. 0763 % -This function still needs to be incorporated into the GUI. 0764 0765 % V2 modifies the code for integration into the VMT GUI. Adds the Rozovskii output to 0766 % the V structure and computes the X components of the primary and secondary 0767 % velocities (in addition to cross stream Y components). 0768 % P.R. Jackson, USGS 0769 0770 0771 % Written by: 0772 % Frank L. Engel (fengel2@illinois.edu) 0773 0774 % Last edited: 6/10/2010 FE 0775 % FE- I fixed an error in computing theta for data in quadrants 3 & 4. The 0776 % linear interpolation of velocity to the bed BC was creating errors 0777 % in the computation of us, so I removed it. Also, I made a new 0778 % variable which is the sum of all vs as an error check (it should 0779 % always sum to zero) 0780 0781 disp('Performing Rozovskii analysis...') 0782 bin_size = A(1,1).Sup.binSize_cm/100; % in meters 0783 0784 for i = 1:size(V.mcsMag,2) 0785 % Finds closest bin to beam avg. depth (ie from V.mcsBed) 0786 [min_difference(i), array_position(i)]... 0787 = min(abs(V.mcsDepth(:,i) - V.mcsBed(i))); 0788 %disp(['ap = ' num2str(array_position(i))]) 0789 % Create a seperate version of the velocity data which can be modified, 0790 % preserving the VMT processing. Replaces all of the NaNs with u=0. 0791 temp_u = V.u; 0792 temp_v = V.v; 0793 n = find(isnan(V.u)); 0794 temp_u(n) = 0; 0795 temp_v(n) = 0; 0796 0797 % Compute Depth-averaged velocities and angles (using a difference 0798 % scheme) 0799 for j = 1:array_position(i) 0800 if j == 1 % Near water surface 0801 % Compute first bin by exprapolating velocity to the water 0802 % surface. WSE = 0. Imposes BC u(z=0) = u(z=bin1) 0803 du_i(j,i) = temp_u(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))... 0804 + temp_u(j,i)*(V.mcsDepth(j,i)-bin_size/2-0); 0805 dv_i(j,i) = temp_v(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))... 0806 + temp_v(j,i)*(V.mcsDepth(j,i)-bin_size/2-0); 0807 elseif j < array_position(i) % Inbetween 0808 du_i(j,i) = temp_u(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2; 0809 dv_i(j,i) = temp_v(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2; 0810 0811 % Got rid of the linear interpolation to the bed- it was 0812 % messing up the Rozovskii secondary velocities (FE 6/10/2010) 0813 elseif j == array_position(i) % Near bed 0814 % k=0; 0815 % % Search bins above the bed for the first bin containing flow 0816 % % data. 0817 % 0818 % % while temp_u(j-k,i) == 0 0819 % % if temp_u(j-k,i) == 0; k = k + 1; else end % find next good bin 0820 % % end 0821 % 0822 indx = find(temp_u(:,i) ~= 0); %Revision PRJ 9-1-09 0823 if isempty(indx) 0824 du_i(:,i) = NaN; 0825 dv_i(:,i) = NaN; 0826 else 0827 l = indx(end); 0828 k = j - l; 0829 % Computes du from last good bin to the bed by linear 0830 % interpolation. IMPOSES BC: u=0 at the bed 0831 % du_i(j-k+1,i) = (temp_u(j-k,i)-temp_u(j,i))/k... 0832 % *(V.mcsDepth(j,i)-V.mcsDepth(j-k,i))/k; 0833 % dv_i(j-k+1,i) = (temp_v(j-k,i)-temp_v(j,i))/k... 0834 % *(V.mcsDepth(j,i)-V.mcsDepth(j-k,i))/k; 0835 0836 % Paints everything below last bin as NaN 0837 du_i(j-k+2:size(V.u,2),i) = NaN; 0838 dv_i(j-k+2:size(V.u,2),i) = NaN; 0839 end 0840 end 0841 end 0842 0843 % Depth averaged quantities 0844 U(i) = nansum(du_i(:,i))/V.mcsDepth(array_position(i),i); 0845 V1(i) = nansum(dv_i(:,i))/V.mcsDepth(array_position(i),i); 0846 U_mag(i) = sqrt(U(i)^2+V1(i)^2); % resultant vector 0847 0848 % Angle of resultant vector from a perpendicular line along the 0849 % transect 0850 phi(i) = atan(V1(i)/U(i)); 0851 phi_deg(i) = rad2deg(phi(i)); 0852 0853 % Compute Rozovskii variables at each bin 0854 for j = 1:array_position(i) 0855 u(j,i) = sqrt(V.u(j,i)^2+V.v(j,i)^2); 0856 if (V.u(j,i) < 0) && (V.v(j,i) < 0) 0857 theta(j,i) = atan(V.v(j,i)/V.u(j,i)) - pi(); 0858 elseif (V.u(j,i) < 0) && (V.v(j,i) > 0) 0859 theta(j,i) = atan(V.v(j,i)/V.u(j,i)) + pi(); 0860 else 0861 theta(j,i) = atan(V.v(j,i)/V.u(j,i)); 0862 end 0863 theta_deg(j,i) = rad2deg(theta(j,i)); 0864 up(j,i) = u(j,i)*cos(theta(j,i)-phi(i)); 0865 us(j,i) = u(j,i)*sin(theta(j,i)-phi(i)); 0866 upy(j,i) = up(j,i)*sin(phi(i)); 0867 upx(j,i) = up(j,i)*cos(phi(i)); 0868 usy(j,i) = us(j,i)*cos(phi(i)); 0869 usx(j,i) = us(j,i)*sin(phi(i)); 0870 depths(j,i) = V.mcsDepth(j,i); 0871 0872 % Compute d_us to check for zero secondary discharge constraint 0873 if j == 1 % Near water surface 0874 dus_i(j,i) = us(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))... 0875 + us(j,i)*(V.mcsDepth(j,i)-bin_size/2-0); 0876 elseif j < array_position(i) % Inbetween 0877 dus_i(j,i) = us(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2; 0878 end 0879 % Sum dus_i - this should be near zero for each ensemble 0880 q_us(i) = nansum(dus_i(:,i)); 0881 end 0882 0883 % Resize variables to be the same as V structure array 0884 indices = j+1:size(V.mcsMag,1); 0885 u(indices,i) = NaN; 0886 theta(indices,i) = NaN; 0887 theta_deg(indices,i) = NaN; 0888 up(indices,i) = NaN; 0889 us(indices,i) = NaN; 0890 upy(indices,i) = NaN; 0891 usy(indices,i) = NaN; 0892 upx(indices,i) = NaN; 0893 usx(indices,i) = NaN; 0894 depths(indices,i) = NaN; 0895 dus_i(indices,i) = NaN; 0896 % u(j+1:size(V.mcsMag,1),i) = NaN; 0897 % theta(j+1:size(V.mcsMag,1),i) = NaN; 0898 % theta_deg(j+1:size(V.mcsMag,1),i) = NaN; 0899 % up(j+1:size(V.mcsMag,1),i) = NaN; 0900 % us(j+1:size(V.mcsMag,1),i) = NaN; 0901 % upy(j+1:size(V.mcsMag,1),i) = NaN; 0902 % usy(j+1:size(V.mcsMag,1),i) = NaN; 0903 % upx(j+1:size(V.mcsMag,1),i) = NaN; 0904 % usx(j+1:size(V.mcsMag,1),i) = NaN; 0905 % depths(j+1:size(V.mcsMag,1),i) = NaN; 0906 % dus_i(j+1:size(V.mcsMag,1),i) = NaN; 0907 end 0908 0909 % Display error message if rozovskii computation of q_us doesn't sum to 0910 % zero 0911 if q_us > 1e-4 0912 disp('Warning: Rozovskii secondary velocities not satisfying continuity!') 0913 else 0914 disp('Computation successfull: Rozovskii secondary velocities satisfy continuity') 0915 end 0916 0917 % Rotate local velocity vectors into global coordinate system by 0918 % determining the angle of the transect using endpoint locations. The 0919 % function "vrotation" is a standard rotation matrix 0920 XStheta = atan((V.mcsY(1,end)-V.mcsY(1,1))/(V.mcsX(1,end)-V.mcsX(1,1))); 0921 XSalpha = XStheta - pi/2; 0922 [ux, uy, uz] = vrotation(V.u,V.v,V.w,XSalpha); 0923 0924 % Put results into the V structure 0925 0926 V.Roz.U = U; 0927 V.Roz.V = V1; 0928 V.Roz.U_mag = U_mag; 0929 V.Roz.phi = phi; 0930 V.Roz.phi_deg = phi_deg; 0931 V.Roz.u = V.u; 0932 V.Roz.v = V.v; 0933 V.Roz.u_mag = u; 0934 V.Roz.depth = depths; 0935 V.Roz.theta = theta; 0936 V.Roz.theta_deg = theta_deg; 0937 V.Roz.up = up; 0938 V.Roz.us = us; 0939 V.Roz.upy = upy; 0940 V.Roz.usy = usy; 0941 V.Roz.upx = upx; 0942 V.Roz.usx = usx; 0943 V.Roz.ux = ux; 0944 V.Roz.uy = uy; 0945 V.Roz.uz = uz; 0946 V.Roz.alpha = XSalpha; 0947 0948 % Rozovskii = struct('Ux', {Ux}, 'Uy', {Uy}, 'U', {U}, 'phi', {phi},... 0949 % 'phi_deg', {phi_deg},'ux', {V.u}, 'uy', {V.v}, 'u', {u}, 'depth', {depths},... 0950 % 'theta', {theta}, 'theta_deg', {theta_deg}, 'up', {up}, 'us', {us},... 0951 % 'upy', {upy}, 'usy', {usy}); 0952 0953 % Name of output file needs to be modified to take handle args from GUI 0954 % Save variable into the VMTProcFile folder 0955 % outfolder = ['VMTProcFiles\']; 0956 % outfile=['Rozovskii']; 0957 % filename = [outfolder outfile]; 0958 % save(filename, 'Rozovskii'); 0959 0960 disp('Rozovskii analysis complete. Added .Roz structure to V data structure.') 0961 % directory = pwd; 0962 % fileloc = [directory '\' filename '.mat']; 0963 % disp(fileloc) 0964