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

- VMT_ProcessTransects Driver program to process multiple transects at a single cross-section

## SUBFUNCTIONS

- function [x,y] = loadUserSetEndpoints()
- function A = interpBadGPS(A,V,z)
- function V = setStation(V)

## 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     [A(zi).Comp.itDist, A(zi).Comp.itDepth] = ...
    0381         meshgrid(A(zi).Comp.dlsortgpsfix,A(zi).Wat.binDepth(:,1));
    0382     %A(zi).Comp.itDepth = A(zi).Wat.binDepth;
    0383
    0384     clear I I2 J J2 bg chk df ed i j nbrk xUTM yUTM n zi...
    0385         den num2 num1 numex
    0386
    0387 end
    0388
    0389 function V = setStation(V)
    0390 if V.phi > 0 && V.phi < 90 %PHI quadrant 1
    0391     V.xLeftBank     = V.xw;
    0392     V.yLeftBank     = V.yn;
    0393     V.xRightBank    = V.xe;
    0394     V.yRightBank    = V.ys;
    0395 elseif V.phi > 90 && V.phi < 180 %PHI quadrant 2
    0396     V.xLeftBank     = V.xe;
    0397     V.yLeftBank     = V.yn;
    0398     V.xRightBank    = V.xw;
    0399     V.yRightBank    = V.ys;
    0400 elseif V.phi > 180 && V.phi < 270 %PHI quadrant 3
    0401     V.xLeftBank     = V.xe;
    0402     V.yLeftBank     = V.ys;
    0403     V.xRightBank    = V.xw;
    0404     V.yRightBank    = V.yn;
    0405 elseif V.phi > 270 && V.phi < 360 %PHI quadrant 4
    0406     V.xLeftBank     = V.xw;
    0407     V.yLeftBank     = V.ys;
    0408     V.xRightBank    = V.xe;
    0409     V.yRightBank    = V.yn;
    0410 elseif V.phi == 0 %Set special cases
    0411     V.xLeftBank     = V.xw;
    0412     V.yLeftBank     = V.yn; %Does not matter if use N or S point (same)
    0413     V.xRightBank    = V.xe;
    0414     V.yRightBank    = V.ys;
    0415 elseif V.phi == 90
    0416     V.xLeftBank     = V.xe; %Does not matter if use E or W point (same)
    0417     V.yLeftBank     = V.yn;
    0418     V.xRightBank    = V.xw;
    0419     V.yRightBank    = V.ys;
    0420 elseif V.phi == 180
    0421     V.xLeftBank     = V.xe;
    0422     V.yLeftBank     = V.yn; %Does not matter if use N or S point (same)
    0423     V.xRightBank    = V.xw;
    0424     V.yRightBank    = V.ys;
    0425 elseif V.phi == 270
    0426     V.xLeftBank     = V.xe; %Does not matter if use E or W point (same)
    0427     V.yLeftBank     = V.ys;
    0428     V.xRightBank    = V.xw;
    0429     V.yRightBank    = V.yn;
    0430 end