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
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)