VMT_GridData2MeanXS

PURPOSE ^

Generates a uniformly spaced grid for the mean cross section and

SYNOPSIS ^

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

DESCRIPTION ^

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

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

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