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

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