# VMT_GridData2MeanXS_AVG

## PURPOSE

This routine generates a uniformly spaced grid for the mean cross section and

## SYNOPSIS

function [A,V] = VMT_GridData2MeanXS(z,A,V)

## DESCRIPTION

     This routine generates a uniformly spaced grid for the mean cross section and
     maps (interpolates) individual transects to this grid.

     AVG version averages the ensembles nearest each grid node (rather than
     interpolating).

     (adapted from code by J. Czuba)

     P.R. Jackson, USGS, 12-9-08

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [A,V] = VMT_GridData2MeanXS(z,A,V)
    0002 % This routine generates a uniformly spaced grid for the mean cross section and
    0003 % maps (interpolates) individual transects to this grid.
    0004 %
    0005 % AVG version averages the ensembles nearest each grid node (rather than
    0006 % interpolating).
    0007 %
    0008 % (adapted from code by J. Czuba)
    0009 %
    0010 % P.R. Jackson, USGS, 12-9-08
    0011
    0012 %% User Input
    0013
    0014 xgdspc = A(1).hgns; %Horizontal Grid node spacing in meters  (vertical grid spacing is set by bins)
    0015
    0016 %% Determine uniform mean c-s grid for vector interpolating
    0017 % Determine the end points of the mean cross-section line
    0018 % Initialize variable with mid range value
    0019 V.xe = mean(A(1).Comp.xm);
    0020 V.ys = mean(A(1).Comp.ym);
    0021 V.xw = mean(A(1).Comp.xm);
    0022 V.yn = mean(A(1).Comp.ym);
    0023
    0024 for zi = 1 : z
    0025
    0026     V.xe = max(max(A(zi).Comp.xm),V.xe);
    0027     V.ys = min(min(A(zi).Comp.ym),V.ys);
    0028     V.xw = min(min(A(zi).Comp.xm),V.xw);
    0029     V.yn = max(max(A(zi).Comp.ym),V.yn);
    0030
    0031 end
    0032
    0033 % Determine the distance between the mean cross-section endpoints
    0034 V.dx = V.xe-V.xw;
    0035 V.dy = V.yn-V.ys;
    0036
    0037 V.dl = sqrt(V.dx^2+V.dy^2);
    0038
    0039 % Determine mean cross-section velocity vector grid
    0040 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)
    0041 V.mcsDepth = A(1).Wat.binDepth(:,1);
    0042 [V.mcsDist V.mcsDepth] = meshgrid(V.mcsDist,V.mcsDepth);
    0043
    0044 % Determine the angle the mean cross-section makes with the
    0045 % x-axis (longitude)
    0046 % Plot mean cross-section line
    0047 if V.m >= 0
    0048     V.theta = atand(V.dy./V.dx);
    0049
    0050     figure(1); hold on
    0051     plot([V.xe V.xw],[V.yn V.ys],'ks'); hold on
    0052
    0053     if V.mfd >= 270 | V.mfd < 90 %Flow to the north
    0054         V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta);            %
    0055         V.mcsY = V.ys+V.mcsDist(1,:).*sind(V.theta);
    0056
    0057     elseif V.mfd >= 90 & V.mfd < 270 %Flow to the south
    0058         V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta);            %
    0059         V.mcsY = V.yn-V.mcsDist(1,:).*sind(V.theta);
    0060     end%
    0061
    0062     plot(V.mcsX,V.mcsY,'k+'); hold on
    0063     plot(V.mcsX(1),V.mcsY(1),'y*'); hold on
    0064
    0065 elseif V.m < 0
    0066     V.theta = -atand(V.dy./V.dx);
    0067
    0068     figure(1); hold on
    0069     plot([V.xe V.xw],[V.ys V.yn],'ks'); hold on
    0070
    0071     if V.mfd >= 270 | V.mfd < 90 %Flow to the north
    0072         V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta);            %
    0073         V.mcsY = V.yn+V.mcsDist(1,:).*sind(V.theta);
    0074
    0075     elseif V.mfd >= 90 & V.mfd < 270 %Flow to the south
    0076         V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta);
    0077         V.mcsY = V.ys-V.mcsDist(1,:).*sind(V.theta);
    0078     end%
    0079
    0080     plot(V.mcsX,V.mcsY,'k+'); hold on
    0081     plot(V.mcsX(1),V.mcsY(1),'y*'); hold on
    0082
    0083 end
    0084
    0085 V.mcsX = meshgrid(V.mcsX,V.mcsDepth(:,1));
    0086 V.mcsY = meshgrid(V.mcsY,V.mcsDepth(:,1));
    0087
    0088 clear zi
    0089
    0090 %% Determine location of mapped ensemble points for interpolating
    0091 % For all transects
    0092
    0093 %A = VMT_DxDyfromLB(z,A,V); %Computes dx and dy
    0094
    0095 for zi = 1 : z
    0096
    0097     % Determine if the mean cross-section line trends NW-SE or SW-NE
    0098     % Determine the distance in radians from the left bank mean
    0099     % cross-section point to the mapped ensemble point for an individual
    0100     % transect
    0101    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)
    0102
    0103     if V.m > 0
    0104         A(zi).Comp.dy = abs(V.yn-A(zi).Comp.ym);
    0105
    0106     elseif V.m < 0
    0107             A(zi).Comp.dy = abs(V.ys-A(zi).Comp.ym);
    0108
    0109     end
    0110
    0111     % Determine the distance in meters from the left bank mean
    0112     % cross-section point to the mapped ensemble point for an individual
    0113     % transect
    0114     A(zi).Comp.dl = sqrt(A(zi).Comp.dx.^2+A(zi).Comp.dy.^2);
    0115
    0116     % Sort vectors by dl
    0117     A(zi).Comp.dlsort = sort(A(zi).Comp.dl,'ascend');
    0118
    0119     % Map indices
    0120     for i = 1 : A(zi).Sup.noe
    0121         for k = 1 : A(zi).Sup.noe
    0122
    0123             if A(zi).Comp.dlsort(i,1) == A(zi).Comp.dl(k,1)
    0124                 A(zi).Comp.vecmap(i,1) = k;
    0125
    0126             end
    0127         end
    0128     end
    0129
    0130     % GPS position fix
    0131     % if distances from the left bank are the same for two ensembles the
    0132     % the position of the right most ensemble is interpolated from adjacent
    0133     % ensembles
    0134     % check for repeat values of distance
    0135     sbt(:,1)=diff(A(zi).Comp.dlsort);
    0136     chk(1,1)=1;
    0137     chk(2:A(zi).Sup.noe,1)=sbt(1:end,1);
    0138
    0139     % identify repeat values
    0140     A(zi).Comp.sd = (chk==0) > 0;
    0141
    0142     % if repeat values exist interpolate distances from adjacent ensembles
    0143     if sum(A(zi).Comp.sd) > 0
    0144
    0145         % bracket repeat sections
    0146         [I,J] = ind2sub(size(A(zi).Comp.sd),find(A(zi).Comp.sd==1));
    0147         df=diff(I);
    0148         nbrk=sum(df>1)+1;
    0149         [I2,J2] = ind2sub(size(df),find(df>1));
    0150
    0151         bg(1)=(I(1)-1);
    0152
    0153         for n = 2 : nbrk
    0154             bg(n)=(I(I2(n-1)+1)-1);
    0155         end
    0156
    0157         for n = 1 : nbrk -1
    0158             ed(n)=(I(I2(n))+1);
    0159         end
    0160
    0161         ed(nbrk)=I(end)+1;
    0162
    0163         % interpolate repeat values
    0164         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
    0165
    0166         for i = 1 : nbrk
    0167             for j = bg(i)+1 : ed(i)-1
    0168                 % interpolate
    0169                 if bg(i) > 0 && ed(i) < length(A(zi).Nav.lat_deg)
    0170
    0171                     den=(ed(i)-bg(i));
    0172                     num2=j-bg(i);
    0173                     num1=ed(i)-j;
    0174
    0175                     A(zi).Comp.dlsortgpsfix(j,1)=...
    0176                         (num1/den).*A(zi).Comp.dlsort(bg(i))+...
    0177                         (num2/den).*A(zi).Comp.dlsort(ed(i));
    0178
    0179                 end
    0180
    0181                 % extrapolate end
    0182                 if ed(i) > length(A(zi).Nav.lat_deg)
    0183
    0184                     numex=ed(i)-length(A(zi).Nav.lat_deg);
    0185
    0186                     A(zi).Comp.dlsortgpsfix(j,1)=numex.*...
    0187                         (A(zi).Comp.dlsort(bg(i))-...
    0188                         A(zi).Comp.dlsort(bg(i)-1))+...
    0189                         A(zi).Comp.dlsort(bg(i));
    0190
    0191                 end
    0192             end
    0193         end
    0194
    0195     else
    0196
    0197         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
    0198
    0199     end
    0200
    0201     % Determine velocity vector grid for individual transects
    0202     [A(zi).Comp.itDist A(zi).Comp.itDepth] = ...
    0203         meshgrid(A(zi).Comp.dlsortgpsfix,A(zi).Wat.binDepth(:,1));
    0204
    0205     clear I I2 J J2 bg chk df ed i j nbrk sbt xUTM yUTM n zi...
    0206         den num2 num1 numex
    0207
    0208 end
    0209
    0210 clear zi i k check
    0211 %% Average ensembles from individual transects using a spatial average to points on the uniform mean c-s grid
    0212 % Fill in uniform grid by averaging ensembles (of individual transects mapped onto the mean
    0213 % cross-section) within one-half a grid spacing on either side of a uniform mean
    0214 % c-s node. This uses all of the ensembles.
    0215
    0216 for zi = 1 : z
    0217 % reorder the ensembles so the distance increments across the c-s
    0218 A(zi).Clean.bsvecmap = A(zi).Clean.bs(:,A(zi).Comp.vecmap);
    0219 A(zi).Clean.vDirvecmap = A(zi).Clean.vDir(:,A(zi).Comp.vecmap);
    0220 A(zi).Clean.vMagvecmap = A(zi).Clean.vMag(:,A(zi).Comp.vecmap);
    0221 A(zi).Clean.vEastvecmap = A(zi).Clean.vEast(:,A(zi).Comp.vecmap);
    0222 A(zi).Clean.vNorthvecmap = A(zi).Clean.vNorth(:,A(zi).Comp.vecmap);
    0223 A(zi).Clean.vVertvecmap = A(zi).Clean.vVert(:,A(zi).Comp.vecmap);
    0224 A(zi).Clean.depthvecmap = nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2)';
    0225
    0226 % determine one-half grid spacing on either side of a node
    0227 plus=V.mcsDist(1,:)+(V.mcsDist(1,2)/2);
    0228 minus=V.mcsDist(1,:)-(V.mcsDist(1,2)/2);
    0229
    0230 for j=1:size(plus,2);%columns
    0231
    0232     indx = find(A(zi).Comp.itDist(1,:)<plus(1,j)&A(zi).Comp.itDist(1,:)>minus(1,j)); % Indices within bin
    0233
    0234     % spatial average ensembles on either side of the nodes
    0235     A(zi).Comp.mcsBack(:,j)=nanmean(A(zi).Clean.bsvecmap(:,indx),2);
    0236     % count the non-nan values that were averaged
    0237     A(zi).Comp.mcsBackContrib(:,j)=sum(~isnan(A(zi).Clean.bsvecmap(:,indx)),2);
    0238
    0239     A(zi).Comp.mcsDir(:,j)=nanmean(A(zi).Clean.vDirvecmap(:,indx),2);
    0240     A(zi).Comp.mcsDirContrib(:,j)=sum(~isnan(A(zi).Clean.vDirvecmap(:,indx)),2);
    0241
    0242     A(zi).Comp.mcsMag(:,j)=nanmean(A(zi).Clean.vMagvecmap(:,indx),2);
    0243     A(zi).Comp.mcsMagContrib(:,j)=sum(~isnan(A(zi).Clean.vMagvecmap(:,indx)),2);
    0244
    0245     A(zi).Comp.mcsEast(:,j)=nanmean(A(zi).Clean.vEastvecmap(:,indx),2);
    0246     A(zi).Comp.mcsEastContrib(:,j)=sum(~isnan(A(zi).Clean.vEastvecmap(:,indx)),2);
    0247
    0248     A(zi).Comp.mcsNorth(:,j)=nanmean(A(zi).Clean.vNorthvecmap(:,indx),2);
    0249     A(zi).Comp.mcsNorthContrib(:,j)=sum(~isnan(A(zi).Clean.vNorthvecmap(:,indx)),2);
    0250
    0251     A(zi).Comp.mcsVert(:,j)=nanmean(A(zi).Clean.vVertvecmap(:,indx),2);
    0252     A(zi).Comp.mcsVertContrib(:,j)=sum(~isnan(A(zi).Clean.vVertvecmap(:,indx)),2);
    0253
    0254     A(zi).Comp.mcsBed(:,j)=nanmean(A(zi).Clean.depthvecmap(:,indx),2);
    0255     A(zi).Comp.mcsBedContrib(:,j)=sum(~isnan(A(zi).Clean.depthvecmap(:,indx)),2);
    0256
    0257 end
    0258
    0259 A(zi).Comp.mcsBack(A(zi).Comp.mcsBack>=255) = NaN;
    0260
    0261 end
    0262
    0263 %% Interpolate individual transects onto uniform mean c-s grid
    0264 % Fill in uniform grid based on individual transects mapped onto the mean
    0265 % cross-section by interpolating between adjacent points
    0266
    0267 %ZI = interp2(X,Y,Z,XI,YI)
    0268 for zi = 1 : z
    0269
    0270     A(zi).Comp.mcsBackI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0271         A(zi).Clean.bs(:,A(zi).Comp.vecmap),V.mcsDist, V.mcsDepth);
    0272     A(zi).Comp.mcsBackI(A(zi).Comp.mcsBackI>=255) = NaN;
    0273
    0274     A(zi).Comp.mcsDirI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0275         A(zi).Clean.vDir(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
    0276     A(zi).Comp.mcsMagI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0277         A(zi).Clean.vMag(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
    0278     A(zi).Comp.mcsEastI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0279         A(zi).Clean.vEast(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
    0280     A(zi).Comp.mcsNorthI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0281         A(zi).Clean.vNorth(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
    0282     A(zi).Comp.mcsVertI = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
    0283         A(zi).Clean.vVert(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
    0284
    0285     A(zi).Comp.mcsBedI  = interp1(A(zi).Comp.itDist(1,:),...
    0286         nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2),V.mcsDist(1,:));
    0287
    0288 end
    0289
    0290 clear zi