VMT_GridData2MeanXS

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 This routine 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
 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] = VMT_GridData2MeanXS(z,A,V,unitQcorrection)
0002 % This routine 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 % Last modified: F.L. Engel, USGS 2/20/2013
0009 
0010 %% User Input
0011 
0012 xgdspc      = A(1).hgns; %Horizontal Grid node spacing in meters
0013 ygdspc      = A(1).vgns; %double(A(1).Sup.binSize_cm)/100; %Vertical Grid node spacing in meters
0014 
0015 % For now, take the minumum bin size as the vertical resolution.
0016 % Later, I will specify this in the GUI (FLE)
0017 %ygdspc      = min(min(double(A(1).Sup.binSize_cm)/100)); %Vertical Grid node spacing in meters
0018 if 0
0019     xgdspc  = V.meddens + V.stddens;  %Auto method should include 67% of the values
0020     %disp(['X Grid Node Auto Spacing = ' num2str(xgdspc) ' m'])
0021     log_text = ['X Grid Node Auto Spacing = ' num2str(xgdspc) ' m'];
0022 end
0023 
0024 
0025 %% Determine uniform mean c-s grid for vector interpolating
0026 
0027 % Determine mean cross-section velocity vector grid NEW: Allowed for
0028 % explicit specification of vertical grid node spacing Also, using linspace
0029 % doesn't necessarily give exactly hgns spaced grid nodes, so the method
0030 % has been adjusted (this is important for user that want to output to a
0031 % model grid). A fragment of length<xgdsoc may be truncated. The impact on
0032 % this for data analysis should be minor, but should be investigated [FLE
0033 % 3/25/2014]
0034 % V.mcsDist               = linspace(0,V.dl,floor(V.dl/xgdspc));
0035 V.mcsDist               = 0:xgdspc:V.dl;
0036 V.mcsDepth              = ...
0037     min(A(1).Wat.binDepth(:)):ygdspc:max(A(1).Wat.binDepth(:));
0038 % V.mcsDepth              = A(1).Wat.binDepth(:,1);
0039 [V.mcsDist, V.mcsDepth] = meshgrid(V.mcsDist,V.mcsDepth');
0040 
0041 
0042 % Define the MCS XY points. (REVISED PRJ, 10-18-12)
0043 % Coordinate assignments depend on the starting
0044 % point and the slope of the cross section. Theta is limited to 0 to 180
0045 % (geographic) and 90 to 270 (arithmetic).  For COS, arithmetic angles
0046 % between 90 and 270 are always negative so no need to add additional IF
0047 % statement based on the slope.  However, SIN theta (aritmetic) is positive
0048 % in MFD quadrants 2 and 4 and negative in 1 and 3. Therefore, we use the slope
0049 % (positive in MFD quadrants 1 and 3, negative in 2 and 4) to determine whether to add or
0050 % subtract the incremental distances from the start point.  (MFD = mean
0051 % flow direction, used to define quadrants above and below)
0052 
0053 if V.xLeftBank == V.xe % MFD Quadrants 2 and 3 (east start)
0054     V.mcsX = V.xLeftBank - V.mcsDist(1,:).*cosd(geo2arideg(V.theta));
0055 else % MFD Quadrants 1 and 4 (west start)
0056     V.mcsX = V.xLeftBank + V.mcsDist(1,:).*cosd(geo2arideg(V.theta));
0057 end
0058 
0059 if V.yLeftBank == V.yn % MFD Quadrants 1 and 2 (north start)
0060     if V.m >= 0 %MFD Quadrant 2
0061         V.mcsY = V.yLeftBank - V.mcsDist(1,:).*sind(geo2arideg(V.theta)); 
0062     else %MFD Quadrant 1
0063         V.mcsY = V.yLeftBank + V.mcsDist(1,:).*sind(geo2arideg(V.theta));
0064     end
0065 else % MFD Quadrants 3 and 4 (south start)
0066     if V.m >= 0 %MFD Quadrant 4
0067         V.mcsY = V.yLeftBank + V.mcsDist(1,:).*sind(geo2arideg(V.theta)); 
0068     else %MFD Quadrant 3
0069         V.mcsY = V.yLeftBank - V.mcsDist(1,:).*sind(geo2arideg(V.theta));
0070     end
0071 end
0072 
0073 V.mcsX = meshgrid(V.mcsX,V.mcsDepth(:,1));
0074 V.mcsY = meshgrid(V.mcsY,V.mcsDepth(:,1));
0075 
0076 
0077 % %Plot the MCS on figure 1
0078 % figure(1); hold on
0079 % plot(V.xLeftBank,V.yLeftBank,'gs','MarkerFaceColor','g'); hold on  %Green left bank start point
0080 % plot(V.xRightBank,V.yRightBank,'rs','MarkerFaceColor','r'); hold on %Red right bank end point
0081 % plot(V.mcsX(1,:),V.mcsY(1,:),'k+'); hold on
0082 % figure(1); set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1])
0083 % clear zi
0084 %
0085 % % Format the ticks for UTM and allow zooming and panning
0086 % figure(1);
0087 % ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM
0088 % hdlzm_fig1 = zoom;
0089 % set(hdlzm_fig1,'ActionPostCallback',@mypostcallback_zoom);
0090 % set(hdlzm_fig1,'Enable','on');
0091 % hdlpn_fig1 = pan;
0092 % set(hdlpn_fig1,'ActionPostCallback',@mypostcallback_pan);
0093 % set(hdlpn_fig1,'Enable','on');
0094 
0095 
0096 %% If specified, correct the streamwise velocity by enforcing mass
0097 % flux (capacitor) continuity
0098 if unitQcorrection
0099     A = VMT_unitQcont(A,V,z);
0100 end
0101 
0102 %% Interpolate individual transects onto uniform mean c-s grid
0103 % Fill in uniform grid based on individual transects mapped onto the mean
0104 % cross-section by interpolating between adjacent points
0105 %
0106 % Originally, the data matrices contained monotonic, finite values. With
0107 % the inclusion of RiverRay data, this assuption is no longer valid. To
0108 % interpolate the quasi-scattered data requires vectorizing and isolating
0109 % the valid data, and then using the TriScatteredInterp class. Ultimately,
0110 % this method should produce the same results with RioGrande data, but with
0111 % marked speed improvements (vecotrized operations are faster in Matlab).
0112 % NOTE: TriScattertedInterp has been replaced by scatteredInterpolant in
0113 % R2013+
0114 % [FLE 3/25/2014]
0115 
0116 XI = V.mcsDist(:);
0117 YI = V.mcsDepth(:);
0118 
0119 %ZI = interp2(X,Y,Z,XI,YI)
0120 for zi = 1 : z
0121     % Vectorize inputs to interp2, index valid data, and preallocate the
0122     % result vectors
0123     [nrows,ncols] = size(A(zi).Comp.itDist);
0124     X             = A(zi).Comp.itDist; 
0125     Y             = A(zi).Comp.itDepth;
0126     valid         = ~isnan(X) & ~isnan(Y);
0127            
0128     % Inputs
0129     bs  = A(zi).Clean.bs(:,A(zi).Comp.vecmap);
0130     vE  = A(zi).Clean.vEast(:,A(zi).Comp.vecmap);
0131     vN  = A(zi).Clean.vNorth(:,A(zi).Comp.vecmap);
0132     vV  = A(zi).Clean.vVert(:,A(zi).Comp.vecmap);
0133     vEr = A(zi).Clean.vError(:,A(zi).Comp.vecmap);
0134     
0135     % Create scatteredInterpolant class
0136     F = TriScatteredInterp(X(valid),Y(valid),bs(valid));
0137     
0138     % Interpolate to each output
0139     mcsBack  = F(XI,YI);
0140     F.V      = vE(valid);
0141     mcsEast  = F(XI,YI);
0142     F.V      = vN(valid);
0143     mcsNorth = F(XI,YI);
0144     F.V      = vV(valid);
0145     mcsVert  = F(XI,YI);
0146     F.V      = vEr(valid);
0147     mcsError = F(XI,YI);
0148     
0149     % Reshape and save to outputs
0150     A(zi).Comp.mcsBack  = reshape(mcsBack  ,size(V.mcsX));
0151     A(zi).Comp.mcsEast  = reshape(mcsEast  ,size(V.mcsX));
0152     A(zi).Comp.mcsNorth = reshape(mcsNorth ,size(V.mcsX));
0153     A(zi).Comp.mcsVert  = reshape(mcsVert  ,size(V.mcsX));
0154     A(zi).Comp.mcsError = reshape(mcsError ,size(V.mcsX));
0155     
0156     %A(zi).Comp.mcsBack = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
0157     %    A(zi).Clean.bs(:,A(zi).Comp.vecmap),V.mcsDist, V.mcsDepth);
0158     %A(zi).Comp.mcsBack(A(zi).Comp.mcsBack>=255) = NaN;
0159     %A(zi).Comp.mcsEast = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
0160     %    A(zi).Clean.vEast(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
0161     %A(zi).Comp.mcsNorth = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
0162     %    A(zi).Clean.vNorth(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
0163     %A(zi).Comp.mcsVert = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
0164     %    A(zi).Clean.vVert(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);
0165     
0166     %Compute magnitude
0167     A(zi).Comp.mcsMag = sqrt(A(zi).Comp.mcsEast.^2 + A(zi).Comp.mcsNorth.^2);
0168     
0169     
0170     %For direction, compute from the velocity components
0171     A(zi).Comp.mcsDir = ari2geodeg((atan2(A(zi).Comp.mcsNorth,A(zi).Comp.mcsEast))*180/pi); 
0172     
0173     A(zi).Comp.mcsBed  = interp1(A(zi).Comp.itDist(1,:),...
0174         nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2),V.mcsDist(1,:));
0175     
0176 end
0177 
0178 % clear zi
0179 
0180 %%%%%%%%%%%%%%%%
0181 % SUBFUNCTIONS %
0182 %%%%%%%%%%%%%%%%
0183 function mypostcallback_zoom(obj,evd)
0184 ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM (when zooming)
0185 
0186 function mypostcallback_pan(obj,evd)
0187 ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM (when panning)

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