VMT_MBBathy

PURPOSE ^

Computes the multibeam bathymetry from the four beams of the ADCP

SYNOPSIS ^

function [A] = VMT_MBBathy(z,A,savefile,beamAng,magVar,wsedata,saveaux)

DESCRIPTION ^

 Computes the multibeam bathymetry from the four beams of the ADCP
 using a script by D.Mueller (USGS). Beam depths are computed for each
 transect prior to any averaging or mapping.

 Added the ability to export timestamps, pitch, roll, and heading
 (2/1/10)

 Added an identifier for each individual transect to the csv output
(FEL, 6/14/12)

 V3 Adds capability to handle timeseries of WSE, PRJ 8-7-12

 P.R. Jackson, USGS, 8-7-12

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [A] = VMT_MBBathy(z,A,savefile,beamAng,magVar,wsedata,saveaux)
0002 % Computes the multibeam bathymetry from the four beams of the ADCP
0003 % using a script by D.Mueller (USGS). Beam depths are computed for each
0004 % transect prior to any averaging or mapping.
0005 %
0006 % Added the ability to export timestamps, pitch, roll, and heading
0007 % (2/1/10)
0008 %
0009 % Added an identifier for each individual transect to the csv output
0010 %(FEL, 6/14/12)
0011 %
0012 % V3 Adds capability to handle timeseries of WSE, PRJ 8-7-12
0013 %
0014 % P.R. Jackson, USGS, 8-7-12
0015 
0016 %% Start
0017 try
0018     %disp('Computing corrected beam depths')
0019     if isstruct(wsedata)
0020         if length(wsedata.elev) == 1
0021             %disp('WSE is a constant value')
0022             wsefiletype = 0;
0023         else
0024             %disp('WSE is a timeseries')
0025             wsefiletype = 1;
0026         end
0027     else
0028         %disp('WSE is a constant value')
0029         warning off
0030         wsedata.elev = wsedata;
0031         wsefiletype = 0;
0032         warning on
0033     end
0034     
0035     %% Step through each transect in the given set
0036     %figure(1); clf
0037     for zi = 1 : z
0038         
0039         % Work on the WSE if a vector
0040         %WSE vector must have a value for each ensemble, so interpolate given
0041         %values to ensemble times
0042         
0043         if wsefiletype  %only process as vector if loaded file rather than single value
0044             %Build an ensemble time vector
0045             enstime = datenum([A(zi).Sup.year+2000 A(zi).Sup.month A(zi).Sup.day...
0046                 A(zi).Sup.hour A(zi).Sup.minute (A(zi).Sup.second+A(zi).Sup.sec100./100)]);
0047             %Had to add 2000 to year--will not work for years < 2000
0048             %Check the times (for debugging)
0049             if 0
0050                 obs_start = datestr(wsedata.obstime(1))
0051                 obs_end = datestr(wsedata.obstime(end))
0052                 
0053                 ens_start = datestr(enstime(1))
0054                 ens_end = datestr(enstime(end))
0055             end
0056             
0057             % Interpolate the WSE values to the ENS times
0058             wse = interp1(wsedata.obstime,wsedata.elev,enstime);
0059             % Plot for debugging
0060             if 0
0061                 figure(1); hold on
0062                 plot(enstime,wse,'k-')
0063                 datetick('x',13)
0064                 ylabel('WSE, in meters')
0065             end
0066         else
0067             wse = wsedata.elev; %Single value
0068         end
0069         
0070         % Compute position and elevation of each beam depth
0071         [exyz] = depthxyz(A(zi).Nav.depth,A(zi).Sup.draft_cm,...
0072             A(zi).Sensor.pitch_deg,A(zi).Sensor.roll_deg,....
0073             A(zi).Sensor.heading_deg,beamAng,...
0074             'm',A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,wse,A(zi).Sup.ensNo);  %magVar,removed 4-8-10
0075         
0076         %Build the auxillary data matrix
0077         if saveaux
0078             auxmat = [A(zi).Sup.year+2000 A(zi).Sup.month A(zi).Sup.day...
0079                 A(zi).Sup.hour A(zi).Sup.minute (A(zi).Sup.second+A(zi).Sup.sec100./100) ...
0080                 A(zi).Sensor.heading_deg A(zi).Sensor.pitch_deg A(zi).Sensor.roll_deg ...
0081                 repmat(zi,size(A(zi).Sup.year))]; % Added transect #(zi) FLE 6/14/12    %Had to add 2000 to year--will not work for years < 2000
0082             auxmat2 = [];
0083             for i = 1:length(A(zi).Sup.ensNo);
0084                 dum = repmat(auxmat(i,:),4,1);
0085                 auxmat2 = cat(1,auxmat2,dum);
0086             end
0087             clear auxmat dum enstime wse
0088         end
0089         
0090         % Store results
0091         idxmbb = find(~isnan(exyz(:,4))& ~isnan(exyz(:,2)));
0092         if zi==1
0093             zmbb=[exyz(idxmbb,1) exyz(idxmbb,2)...
0094                 exyz(idxmbb,3) exyz(idxmbb,4)];
0095             if saveaux
0096                 auxmbb = auxmat2(idxmbb,:);
0097             end
0098         else
0099             zmbb=cat(1,zmbb,[exyz(idxmbb,1)...
0100                 exyz(idxmbb,2) exyz(idxmbb,3) exyz(idxmbb,4)]);
0101             if saveaux
0102                 auxmbb = cat(1,auxmbb,auxmat2(idxmbb,:));
0103             end
0104         end
0105         
0106         A(zi).Comp.exyz = exyz(idxmbb,:);
0107         
0108         
0109         clear idxmbb exyz;
0110         %pack;
0111     end
0112     
0113     
0114     %% Save the data
0115     
0116     if 1
0117         %disp('Exporting multibeam bathymetry')
0118         %disp([savefile(1:end-4) '_mbxyz.csv'])
0119         outfile = [savefile(1:end-4) '.csv'];
0120         if saveaux
0121             outmat = [zmbb auxmbb];
0122             ofid   = fopen(outfile, 'wt');
0123             outcount = fprintf(ofid,'EnsNo,     Easting_WGS84_m,    Northing_WGS84_m,  Elev_m,  Year,  Month,  Day,  Hour,  Minute,  Second,  Heading_deg,  Pitch_deg,  Roll_deg,  Transect\n'); % Modified to output transect # FLE 6/14/12
0124             outcount = fprintf(ofid,'%6.0f, %14.2f, %14.2f, %8.2f, %4.0f, %2.0f, %2.0f, %2.0f, %2.0f, %2.2f, %3.3f, %3.3f, %3.3f, %3.0f\n',outmat');
0125             fclose(ofid);
0126         else
0127             outmat = zmbb;
0128             ofid   = fopen(outfile, 'wt');
0129             outcount = fprintf(ofid,'EnsNo,     Easting_WGS84_m,    Northing_WGS84_m,  Elev_m\n');
0130             outcount = fprintf(ofid,'%6.0f, %14.2f, %14.2f, %8.2f\n',outmat');
0131             fclose(ofid);
0132         end
0133         %dlmwrite([savefile(1:end-4) '_mbxyz.csv'],outmat,'precision',15);
0134     end
0135 catch err
0136      if isdeployed
0137         errLogFileName = fullfile(pwd,...
0138             ['errorLog' datestr(now,'yyyymmddHHMMSS') '.txt']);
0139         msgbox({['An unexpected error occurred. Error code: ' err.identifier];...
0140             ['Error details are being written to the following file: '];...
0141             errLogFileName},...
0142             'VMT Status: Unexpected Error',...
0143             'error');
0144         fid = fopen(errLogFileName,'W');
0145         fwrite(fid,err.getReport('extended','hyperlinks','off'));
0146         fclose(fid);
0147         rethrow(err)
0148     else
0149         msgbox(['An unexpected error occurred. Error code: ' err.identifier],...
0150             'VMT Status: Unexpected Error',...
0151             'error');
0152         rethrow(err);
0153     end    
0154 end

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