VMT_MapEns2MeanXS

PURPOSE ^

Fits multiple transects at a single location with a single

SYNOPSIS ^

function [A,V,log_text] = VMT_MapEns2MeanXSV2(z,A,setends)

DESCRIPTION ^

 Fits multiple transects at a single location with a single
 line and maps individual ensembles to this line. Inputs are number of
 files (z) and data matrix (Z)(see ReadFiles.m). Output is the updated
 data matrix with new mapped variables.

 If supplied, the function has the capability to set the endpoints of the
 mean cross section

 (adapted from code by J. Czuba)

 P.R. Jackson, USGS, 12-9-08 
 Last modified: F.L. Engel, USGS, 2/20/2013

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [A,V,log_text] = VMT_MapEns2MeanXSV2(z,A,setends)
0002 % Fits multiple transects at a single location with a single
0003 % line and maps individual ensembles to this line. Inputs are number of
0004 % files (z) and data matrix (Z)(see ReadFiles.m). Output is the updated
0005 % data matrix with new mapped variables.
0006 %
0007 % If supplied, the function has the capability to set the endpoints of the
0008 % mean cross section
0009 %
0010 % (adapted from code by J. Czuba)
0011 %
0012 % P.R. Jackson, USGS, 12-9-08
0013 % Last modified: F.L. Engel, USGS, 2/20/2013
0014 
0015 
0016 
0017 %% Determine the best fit mean cross-section line from multiple transects
0018 % initialize vectors for concatenation
0019 
0020 x = [];
0021 y = [];
0022 % figure(1); clf
0023 % set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1])
0024 for zi = 1 : z
0025     
0026     % concatenate coords into a single column vector for regression
0027     x=cat(1,x,A(zi).Comp.xUTM);
0028     y=cat(1,y,A(zi).Comp.yUTM);
0029 
0030 %     figure(1); hold on
0031 %     plot(A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,'b'); hold on
0032     
0033     % Plot the various reject and/or adjusted GPS location data for
0034     % reference
0035     %plot(A(zi).Comp.xUTM,A(zi).Comp.yUTM,'r'); hold on
0036 %     plot(...
0037 %         ...A(zi).Comp.xUTMraw(A(zi).Comp.gps_reject_locations),...
0038 %         ...A(zi).Comp.yUTMraw(A(zi).Comp.gps_reject_locations),'g.',...
0039 %         ...A(zi).Comp.xUTMraw(A(zi).Comp.gps_repeat_locations),...
0040 %         ...A(zi).Comp.yUTMraw(A(zi).Comp.gps_repeat_locations),'y.',...
0041 %         A(zi).Comp.xUTMraw(A(zi).Comp.gps_fly_aways),...
0042 %         A(zi).Comp.yUTMraw(A(zi).Comp.gps_fly_aways),'r.')
0043     
0044     %Find the mean east and north velocity for eact transect (for mean flow
0045     %direction)
0046     mVe(zi) = nanmean(nanmean(A(zi).Clean.vEast,1));
0047     mVn(zi) = nanmean(nanmean(A(zi).Clean.vNorth,1));
0048               
0049 end
0050 
0051 % Gets a user text file with fixed cross section end points
0052 if setends
0053     [x,y] = loadUserSetEndpoints(); % subfunction
0054 %     figure(1); hold on
0055 %     plot(x,y,'go','MarkerSize',10); hold on
0056     
0057     % Save the shorepath
0058     % if exist('LastDir.mat') == 2
0059         % save('LastDir.mat','endspath','-append')
0060     % else
0061         % save('LastDir.mat','endspath')
0062     % end
0063 end
0064 
0065 % Compute the mean flow direction for all the transects in geographic angle
0066 mfdVe = mean(mVe);  %Mean east velocity for all the transects
0067 mfdVn = mean(mVn);  %Mean north velocity for all the transects
0068 V.mfd = ari2geodeg((atan2(mfdVn, mfdVe))*180/pi);  % Mean flow direction in geographic angle
0069 
0070 % find the equation of the best fit line
0071 xw = min(x); xe = max(x);
0072 ys = min(y); yn = max(y);
0073 xrng = xe - xw;
0074 yrng = yn - ys;
0075 
0076 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
0077     [P,S] = polyfit(x,y,1);
0078 %     figure(1); hold on;
0079 %     plot(x,polyval(P,x),'g-')
0080     V.m = P(1);
0081     V.b = P(2);
0082     dx = xe-xw;
0083     dy = polyval(P,xe)-polyval(P,xw);
0084 else
0085     [P,S] = polyfit(y,x,1);
0086 %     figure(1); hold on;
0087 %     plot(polyval(P,y),y,'g-')
0088     V.m = 1/P(1);           %Reformat slope and intercept in terms of y= fn(x) rather than x = fn(y)
0089     V.b = -P(2)/P(1);
0090     dx = polyval(P,yn)-polyval(P,ys);
0091     dy = yn-ys;
0092 %     if V.m >= 0
0093 %         dy = yn-ys;
0094 %     elseif V.m < 0
0095 %         dy = ys-yn;
0096 %     end
0097 end
0098 
0099 % Determine the distance between the mean cross-section
0100 % endpoints
0101 dl = sqrt(dx^2+dy^2);
0102 
0103 % Compute the angle of the MCS in geographic angle
0104 if V.m >= 0
0105     V.theta = ari2geodeg(atand(V.m)); 
0106 elseif V.m < 0
0107     V.theta = ari2geodeg(atand(V.m));
0108 end
0109 
0110 % Determine the direction of the streamwise coordinate, which
0111 % is taken as perpendicular to the mean cross section. Theta is
0112 % expressed in geographical (N = 0 deg, clockwise positive)
0113 % coordinates. This method uses a vector based approach which
0114 % is insensitive to orientation of the cross section.
0115 
0116 % First compute the normal unit vector to the mean
0117 % cross section
0118 N = [-dy/sqrt(dx^2+dy^2)...
0119     dx/sqrt(dx^2+dy^2)];
0120 
0121 % Compute the mean flow direction in the cross section. To do
0122 % this, we also have to convert from geographic angle to
0123 % arimetic angle
0124 arimfddeg = geo2arideg(V.mfd);
0125 [xmfd,ymfd] = pol2cart(arimfddeg*pi/180,1);
0126 M = [xmfd ymfd];
0127 
0128 % Now compute the angle between the normal and mean flow
0129 % direction unit vectors
0130 vdif = acos(dot(N,M)/(norm(N)*norm(M)))*180/pi;
0131 
0132 % If the angle is greater than 90 degs, the normal vector needs
0133 % to be reversed before resolving the u,v coordinates
0134 if vdif >= 90
0135     N = -N;
0136 end
0137 
0138 % Plot N and M to check (scale of the vectors is 10% of the
0139 % total length of the cross section)
0140 midy = ys+abs(yrng)/2;
0141 midx = xw+xrng/2;
0142 % figure(1); hold on;
0143 % quiver(...
0144 %     midx,midy,N(1)*dl*0.1,...
0145 %     N(2)*dl*0.1,1,'k')
0146 % quiver(...
0147 %     midx,midy,M(1)*dl*0.1,...
0148 %     M(2)*dl*0.1,1,'r')
0149 
0150 % Geographic angle of the normal vector to the cross section
0151 V.phi = ari2geodeg(cart2pol(N(1),N(2))*180/pi);
0152 
0153 clear x y stats whichstats zi
0154 
0155 %% Map ensembles to mean c-s line
0156 % Determine the point (mapped ensemble point) where the equation of the
0157 % mean cross-section line intercepts a line perpendicular to the mean
0158 % cross-section line passing through an ensemble from an individual
0159 % transect (see notes for equation derivation)
0160 
0161 for zi = 1 : z
0162     
0163     A(zi).Comp.xm = ((A(zi).Comp.xUTM-V.m.*V.b+V.m.*A(zi).Comp.yUTM)...
0164         ./(V.m.^2+1));
0165     A(zi).Comp.ym = ((V.b+V.m.*A(zi).Comp.xUTM+V.m.^2.*A(zi).Comp.yUTM)...
0166         ./(V.m.^2+1));
0167 
0168     %A(zi).Comp.h1 = nanmean(A(zi).Nav.depth,2)';
0169 end
0170 
0171 
0172 %Plot data to check
0173 xensall = [];
0174 yensall = [];
0175 for zi = 1 : z
0176 %   plot(A(zi).Comp.xm,A(zi).Comp.ym,'b.')
0177   xensall = [xensall; A(zi).Comp.xm];
0178   yensall = [yensall; A(zi).Comp.ym];
0179 end
0180 % plot(A(3).Comp.xm,A(3).Comp.ym,'xg')
0181 % plot(A(4).Comp.xm,A(4).Comp.ym,'oy')
0182 % xlabel('UTM Easting (m)')
0183 % ylabel('UTM Northing (m)')
0184 % box on
0185 % grid on
0186 %Plot a legend in Figure 1
0187 %figure(1); hold on
0188 %legend('Shoreline','GPS(corr)','GPS(raw)','Best Fit','Trans 1
0189 %(mapped)','Other Trans (mapped)')
0190 
0191 %Compute the median distance between mapped points
0192 Dmat = [xensall yensall];
0193 if xrng > yrng
0194     Dmat = sortrows(Dmat,1);
0195 else
0196     Dmat = sortrows(Dmat,2);
0197 end
0198 dxall = diff(Dmat(:,1));
0199 dyall = diff(Dmat(:,2));
0200 densall = sqrt(dxall.^2 + dyall.^2);
0201 V.meddens = median(densall);
0202 V.stddens = std(densall);
0203 % disp(['Median Spacing Between Mapped Ensembles = ' num2str(V.meddens) ' m'])
0204 % disp(['Standard Deviation of Spacing Between Mapped Ensembles = ' num2str(V.stddens) ' m'])
0205 % disp(['Recommended Grid Node Spacing > ' num2str(V.meddens + V.stddens) ' m'])
0206 
0207 %Display in message box for compiled version
0208 msg_string = {['Median Spacing Between Mapped Ensembles = ' num2str(V.meddens) ' m'],...
0209     ['Standard Deviation of Spacing Between Mapped Ensembles = ' num2str(V.stddens) ' m'],...
0210     ['Recommended Grid Node Spacing > ' num2str(V.meddens + V.stddens) ' m']};
0211 % msgbox(msg_string,'VMT Grid Node Spacing','help','replace');
0212 log_text = {...
0213     ['      Median Spacing Between Mapped Ensembles = ' num2str(V.meddens) ' m'];...
0214     ['      Standard Deviation of Spacing Between Mapped Ensembles = ' num2str(V.stddens) ' m'];...
0215     ['      Recommended Grid Node Spacing > ' num2str(V.meddens + V.stddens) ' m']};
0216     
0217 
0218 %% Determine location of mapped ensemble points for interpolating
0219 % For all transects
0220 
0221 % Determine the end points of the mean cross-section line
0222 % Initialize variable with mid range value
0223 V.xe = mean(A(1).Comp.xm);
0224 V.ys = mean(A(1).Comp.ym);
0225 V.xw = mean(A(1).Comp.xm);
0226 V.yn = mean(A(1).Comp.ym);
0227 
0228 for zi = 1 : z
0229     
0230     V.xe = max(max(A(zi).Comp.xm),V.xe);
0231     V.ys = min(min(A(zi).Comp.ym),V.ys);
0232     V.xw = min(min(A(zi).Comp.xm),V.xw);
0233     V.yn = max(max(A(zi).Comp.ym),V.yn);
0234     
0235 end
0236 
0237 % Determine the distance between the mean cross-section endpoints
0238 V.dx = V.xe-V.xw;
0239 V.dy = V.yn-V.ys;
0240 V.dl = sqrt(V.dx.^2+V.dy.^2);
0241 
0242 % Use the correctly oriented normal vector, or rather V.phi, to set the
0243 % start bank so we are always starting on the left bank looking
0244 % downstream (PRJ, 10-17-12)
0245 V = setStation(V); % Subfunction
0246 
0247 % Occasionally, the GPS location of an esemble will be the same between two
0248 % or more ensembles. Find those ensembles, and interpolate a new GPS
0249 % position based on adjacent good GPS positions.
0250 A = interpBadGPS(A,V,z); % Subfunction
0251 
0252 
0253 % clear zi i k check
0254 
0255 
0256 
0257 
0258 %%%%%%%%%%%%%%%%
0259 % SUBFUNCTIONS %
0260 %%%%%%%%%%%%%%%%
0261 function [x,y] = loadUserSetEndpoints()
0262 defaultpath = 'C:\';
0263 endspath = [];
0264 if 0 %exist('VMT\LastDir.mat') == 2
0265     load('VMT\LastDir.mat');
0266     if exist(endspath) == 7
0267         [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',endspath);
0268     else
0269         [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',defaultpath);
0270     end
0271 else
0272     [file,endspath] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File',defaultpath);
0273 end
0274 infile = [endspath file];
0275 %[file,path] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File');
0276 %infile = [path file];
0277 disp('Loading Endpoint File...' );
0278 disp(infile);
0279 data = dlmread(infile);
0280 x = data(:,1);
0281 y = data(:,2);
0282 
0283 function A = interpBadGPS(A,V,z)
0284 for zi = 1 : z
0285     
0286     % Compute the change in X and Y from the start point to each observation
0287     A(zi).Comp.dx = abs(V.xLeftBank - A(zi).Comp.xm);  %Revised (PRJ, 10-17-12)
0288     A(zi).Comp.dy = abs(V.yLeftBank - A(zi).Comp.ym);
0289     
0290     % Determine the distance in meters from the left bank mean
0291     % cross-section point to the mapped ensemble point for an individual
0292     % transect
0293     A(zi).Comp.dl = sqrt(A(zi).Comp.dx.^2+A(zi).Comp.dy.^2);
0294     
0295     % Sort vectors by dl
0296     [A(zi).Comp.dlsort,A(zi).Comp.vecmap] = sort(A(zi).Comp.dl,'ascend');
0297     
0298     % Map indices  %FIXME  This computation is VERY slow.  Suggest revising
0299     % for speed
0300     % FLE 12/10: This is fixed. The loop was to build a vecor map. This is
0301     % an included output of sort.
0302     %     for i = 1 : A(zi).Sup.noe
0303     %         for k = 1 : A(zi).Sup.noe
0304     %
0305     %             if A(zi).Comp.dlsort(i,1) == A(zi).Comp.dl(k,1)
0306     %                 A(zi).Comp.vecmap(i,1) = k;
0307     %
0308     %             end
0309     %         end
0310     %     end
0311     
0312     % GPS position fix
0313     % if distances from the left bank are the same for two ensembles the
0314     % the position of the right most ensemble is interpolated from adjacent
0315     % ensembles
0316     % check for repeat values of distance
0317     chk(:,1)=[1; diff(A(zi).Comp.dlsort)];
0318     
0319     % identify repeat values
0320     A(zi).Comp.sd = (chk==0) > 0;
0321     
0322     % if repeat values exist interpolate distances from adjacent ensembles
0323     if any(A(zi).Comp.sd)
0324         
0325         % bracket repeat sections
0326         [I,J] = ind2sub(size(A(zi).Comp.sd),find(A(zi).Comp.sd==1));
0327         df=diff(I);
0328         numberBreaks=sum(df>1)+1;
0329         [I2,J2] = ind2sub(size(df),find(df>1));
0330         
0331         idxBeginBracket(1)=(I(1)-1);
0332         
0333         for n = 2 : numberBreaks
0334             idxBeginBracket(n)=(I(I2(n-1)+1)-1);
0335         end
0336         
0337         for n = 1 : numberBreaks -1
0338             idxEndBracket(n)=(I(I2(n))+1);
0339         end
0340         
0341         idxEndBracket(numberBreaks)=I(end)+1;
0342         
0343         % interpolate repeat values
0344         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
0345         
0346         for i = 1 : numberBreaks
0347             for j = idxBeginBracket(i)+1 : idxEndBracket(i)-1
0348                 % interpolate
0349                 if idxBeginBracket(i) > 0 && idxEndBracket(i) < length(A(zi).Nav.lat_deg)
0350                     
0351                     den=(idxEndBracket(i)-idxBeginBracket(i));
0352                     num2=j-idxBeginBracket(i);
0353                     num1=idxEndBracket(i)-j;
0354                     
0355                     A(zi).Comp.dlsortgpsfix(j,1)=...
0356                         (num1/den).*A(zi).Comp.dlsort(idxBeginBracket(i))+...
0357                         (num2/den).*A(zi).Comp.dlsort(idxEndBracket(i));
0358                     
0359                 end
0360                 
0361                 % extrapolate end
0362                 if idxEndBracket(i) > length(A(zi).Nav.lat_deg)
0363                     
0364                     numex=idxEndBracket(i)-length(A(zi).Nav.lat_deg);
0365                     
0366                     A(zi).Comp.dlsortgpsfix(j,1)=numex.*...
0367                         (A(zi).Comp.dlsort(idxBeginBracket(i))-...
0368                         A(zi).Comp.dlsort(idxBeginBracket(i)-1))+...
0369                         A(zi).Comp.dlsort(idxBeginBracket(i));
0370                 end
0371             end
0372         end
0373         
0374     % No duplicate GPS points
0375     else 
0376         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
0377     end
0378     
0379     % Determine velocity vector grid for individual transects
0380     % Previous method used meshgrid. Now, tfile reads the binDepths
0381     % dynamically (for case of RiverRay), thus it's faster to use repmat
0382     % for itDist, and just assign itDepth. [FLE, 3/25/2014]
0383 %     [A(zi).Comp.itDist, ~] = ...
0384 %         meshgrid(A(zi).Comp.dlsortgpsfix,A(zi).Wat.binDepth(:,1));
0385     A(zi).Comp.itDist  = repmat(A(zi).Comp.dlsortgpsfix',size(A(zi).Wat.binDepth,1),1);
0386     A(zi).Comp.itDepth = A(zi).Wat.binDepth(:,A(zi).Comp.vecmap);
0387     %A(zi).Comp.itDepth = A(zi).Wat.binDepth;
0388     
0389     clear I I2 J J2 bg chk df ed i j nbrk xUTM yUTM n zi...
0390         den num2 num1 numex
0391     
0392 end
0393 
0394 function V = setStation(V)
0395 if V.phi > 0 && V.phi < 90 %PHI quadrant 1
0396     V.xLeftBank     = V.xw;
0397     V.yLeftBank     = V.yn;
0398     V.xRightBank    = V.xe;
0399     V.yRightBank    = V.ys;
0400 elseif V.phi > 90 && V.phi < 180 %PHI quadrant 2
0401     V.xLeftBank     = V.xe;
0402     V.yLeftBank     = V.yn;
0403     V.xRightBank    = V.xw;
0404     V.yRightBank    = V.ys;
0405 elseif V.phi > 180 && V.phi < 270 %PHI quadrant 3
0406     V.xLeftBank     = V.xe;
0407     V.yLeftBank     = V.ys;
0408     V.xRightBank    = V.xw;
0409     V.yRightBank    = V.yn;
0410 elseif V.phi > 270 && V.phi < 360 %PHI quadrant 4
0411     V.xLeftBank     = V.xw;
0412     V.yLeftBank     = V.ys;
0413     V.xRightBank    = V.xe;
0414     V.yRightBank    = V.yn;
0415 elseif V.phi == 0 %Set special cases
0416     V.xLeftBank     = V.xw;
0417     V.yLeftBank     = V.yn; %Does not matter if use N or S point (same)
0418     V.xRightBank    = V.xe;
0419     V.yRightBank    = V.ys;
0420 elseif V.phi == 90
0421     V.xLeftBank     = V.xe; %Does not matter if use E or W point (same)
0422     V.yLeftBank     = V.yn;
0423     V.xRightBank    = V.xw;
0424     V.yRightBank    = V.ys;
0425 elseif V.phi == 180
0426     V.xLeftBank     = V.xe;
0427     V.yLeftBank     = V.yn; %Does not matter if use N or S point (same)
0428     V.xRightBank    = V.xw;
0429     V.yRightBank    = V.ys;
0430 elseif V.phi == 270
0431     V.xLeftBank     = V.xe; %Does not matter if use E or W point (same)
0432     V.yLeftBank     = V.ys;
0433     V.xRightBank    = V.xw;
0434     V.yRightBank    = V.yn;
0435 end

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