VMT_ProcessTransectsV3_new

PURPOSE ^

Not used by current version

SYNOPSIS ^

function [A,V] = VMT_ProcessTransectsV3_new(z,A,setends)

DESCRIPTION ^

 Not used by current version
This routine acts as a driver program to process multiple transects at a
single cross-section for velocity mapping.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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