% This is a script that shows how velocity data can be % transformed between BEAM coordinates and ENU coordinates. Beam % coordinates are defined as the velocity measured along the four /five % beams of the instrument. % ENU coordinates are defined in an earth coordinate system, where % E represents the East-West component, N represents the North-South % component and U represents the Up-Down component. % This script can be used for transforming data collected by Signature % instruments. The instruments orientation, its frequency and the existance % of an AHRS are taken into account. % This script is written with .mat files from SIGNATURE DEPLOYMENT % Variables names might vary if exported from Ocean Contour! % load .mat file clear all load('test_plan_1000.mat') %input name of your datafile % ----------------------------------------------------------------------------- % select deployment mode that you want to transform mode='avg' % !!set the declination at deployment site that needs to be corrected for!! declination = 0 % ----------------------------------------------------------------------------- switch mode case 'avg' dataModeWord = 'Average'; configModeWord = 'Average'; case 'avg1' dataModeWord = 'Alt_Average'; configModeWord = 'Alt_Average'; case 'burst1' dataModeWord = 'Alt_Burst'; configModeWord = 'Alt_Burst'; case 'burst' dataModeWord = 'Burst'; configModeWord = 'Burst'; end % define if you want to process a specific index or the entire timeseries index = 1:size(Data.([dataModeWord '_Time']), 1); if isfield(Data, [dataModeWord '_VelUp2']) || isfield(Data, [dataModeWord '_VelBeam4']) || isfield(Data, [dataModeWord '_VelZ2']) twoZs = 1; else twoZs = 0; end % make the assumption the beam mapping is the same for all measurements in a data file activeBeams = Data.( [dataModeWord '_BeamToChannelMapping'] )( 1, : )+1; activeBeams = activeBeams(find(activeBeams > 0)); numberOfBeams = length( activeBeams ); if numberOfBeams <= 2 disp( 'Transformations require at least 3 active beams.' ) T_beam2xyz = nan; return end %% load transformation matrix T beam2xyz_vec = Config.([ configModeWord '_Beam2xyz'] ); T_beam2xyz=beam2xyz_vec % If exported from Ocean Contour, uncomment and reshape T_beam2xyz as needed. %if numberOfBeams == 4, % T_beam2xyz = [beam2xyz_vec(1:4); beam2xyz_vec(5:8); beam2xyz_vec(9:12); beam2xyz_vec(13:16)]; %else % if numberOfBeams == 3, % T_beam2xyz = [beam2xyz_vec(1:3); beam2xyz_vec(4:6); beam2xyz_vec(7:9)]; % end %end % Special case of Signature55 if strcmp(Config.InstrumentName, 'Signature55kHz'), T_beam2xyz =[ 1.9492 -0.9746 -0.9746; ... 0 -1.6881 1.6881; ... 0.3547 0.3547 0.3547]; end % Fix bug in transformation matrix in version 2163 if Config.FWVersion <= 2163 && numberOfBeams == 4, T_beam2xyz(2,:) = -T_beam2xyz(2,:); end %set length of timeseries L=length( Data.( [ dataModeWord '_Time' ] )); % ----------------------------------------------------------------------------- % Start of transformation % ----------------------------------------------------------------------------- % BEAM -> XYZ % check if data is in BEAM if strcmpi( Config.( [ configModeWord '_CoordSystem' ] ), 'BEAM' ) disp( 'Velocity data is in Beam coordinate system.' ) xAllCells = zeros( L, Config.( [ configModeWord '_NCells' ] ) ); yAllCells = zeros( L, Config.( [ configModeWord '_NCells' ] ) ); zAllCells = zeros( L, Config.( [configModeWord '_NCells' ] ) ); if twoZs == 1 z2AllCells = zeros( L, Config.( [ configModeWord '_NCells' ] ) ); end xyz = zeros( size( T_beam2xyz, 2 ), L ); beam = zeros( size( T_beam2xyz, 2 ), L ); for nCell = 1:Config.( [configModeWord '_NCells' ] ) for i = 1:numberOfBeams beamChannel = activeBeams(i); fieldName = [dataModeWord '_VelBeam' num2str(beamChannel)]; beam(i, :) = Data.(fieldName)(index, nCell)'; end %transform BEAM to XYZ xyz = T_beam2xyz * beam; xAllCells( :, nCell ) = xyz( 1, : )'; yAllCells( :, nCell ) = xyz( 2, : )'; zAllCells( :, nCell ) = xyz( 3, : )'; if twoZs == 1 z2AllCells( :, nCell ) = xyz( 4, : )'; end end Config.( [configModeWord '_CoordSystem'] ) = 'XYZ'; Data.( [ dataModeWord '_VelX' ] )(index,1:nCell) = xAllCells; Data.( [ dataModeWord '_VelY' ] )(index,1:nCell) = yAllCells; if twoZs == 1 Data.( [ dataModeWord '_VelZ1' ] )(index,1:nCell) = zAllCells; Data.( [ dataModeWord '_VelZ2' ] )(index,1:nCell) = z2AllCells; else Data.( [ dataModeWord '_VelZ' ] )(index,1:nCell) = zAllCells; end disp( 'BEAM transformed to XYZ' ) end ahrsUsed = false; hprUsed = false; % ----------------------------------------------------------------------------- % check if data is already in ENU if strcmpi( Config.( [ configModeWord '_CoordSystem' ] ), 'ENU' ) disp( 'Velocity data is already in ENU coordinate system.' ) return else if strcmpi( Config.( [ configModeWord '_CoordSystem' ] ), 'XYZ' ) disp( 'Velocity data is in XYZ coordinate system.' ) end K = 3; Ncells = Config.( [configModeWord '_NCells' ] ); EAllCells = zeros( L , Ncells ); NAllCells = zeros( L , Ncells ); UAllCells = zeros( L , Ncells ); if twoZs == 1 U2AllCells = zeros( L, Config.( [ configModeWord '_NCells' ] ) ); K = 4; end Name = ['X','Y','Z']; ENU = zeros( K, Ncells); xyz = zeros( K, Ncells); % shifting T matrix depending on instrument orientation (up or down) for sampleIndex = index(1):index(end) orientation = bitand(bitshift(uint32(Data.( [dataModeWord '_Status' ])(sampleIndex)), -25),7); if (orientation == 5) signXYZ=[1 -1 -1 -1]; else signXYZ=[1 1 1 1]; end % loading AHRS transformation matrix if available if (orientation == 7) if ~ahrsUsed disp('AHRS matrix is used'); ahrsUsed = true; end % AHRS xyz2enuAHRS = transpose(reshape(Data.([dataModeWord '_AHRSRotationMatrix'])(sampleIndex,:),3,3)); % Make declination rotation matrix H = [cosd(declination) sind(declination) 0; -sind(declination) cosd(declination) 0; 0 0 1]; xyz2enuAHRSdecl = H*xyz2enuAHRS; % Update the transformation matrix Data.([dataModeWord '_AHRSRotationMatrix'])(sampleIndex,:)=reshape(transpose(xyz2enuAHRSdecl),1,9); % Update the heading value heading = 90-(atan2(xyz2enuAHRSdecl(2,1),xyz2enuAHRSdecl(1,1))*180/pi); if heading < 0, heading = heading + 360; end if heading >= 360, heading = heading - 360; end Data.([dataModeWord '_Heading'])(sampleIndex) = heading; xyz2enu = xyz2enuAHRSdecl; else if ~hprUsed disp('HPR matrix is used'); hprUsed = true; end % Apply declination heading = Data.([dataModeWord '_Heading'])(sampleIndex) + declination; if heading < 0, heading = heading + 360; end if heading >= 360, heading = heading - 360; end Data.([dataModeWord '_Heading'])(sampleIndex) = heading; % Regular compass hh = pi*(Data.([dataModeWord '_Heading'])(sampleIndex)-90)/180; pp = pi*Data.([dataModeWord '_Pitch'])(sampleIndex)/180; rr = pi*Data.([dataModeWord '_Roll'])(sampleIndex)/180; % Make heading matrix H = [cos(hh) sin(hh) 0; -sin(hh) cos(hh) 0; 0 0 1]; % Make tilt matrix (pitch and roll) P = [cos(pp) -sin(pp)*sin(rr) -cos(rr)*sin(pp);... 0 cos(rr) -sin(rr); ... sin(pp) sin(rr)*cos(pp) cos(pp)*cos(rr)]; % Make resulting transformation matrix xyz2enu = H*P; end if (twoZs == 1) xyz2enu(1,3) = xyz2enu(1,3)/2; xyz2enu(1,4) = xyz2enu(1,3); xyz2enu(2,3) = xyz2enu(2,3)/2; xyz2enu(2,4) = xyz2enu(2,3); xyz2enu(4,:) = xyz2enu(3,:); xyz2enu(3,4) = 0; xyz2enu(4,4) = xyz2enu(3,3); xyz2enu(4,3) = 0; end for i = 1:K if (twoZs == 1) && (i >= 3) axs = [ Name(3) num2str((i-2),1) ]; else axs = Name(i); end xyz( i, : ) = signXYZ(i) * Data.( [ dataModeWord '_Vel' axs] )( sampleIndex, 1:Ncells )'; end % transform XYZ to ENU ENU = xyz2enu * xyz; Data.( [ dataModeWord '_VelEast' ] )( sampleIndex, 1:Ncells ) = ENU( 1, : )'; Data.( [ dataModeWord '_VelNorth' ] )( sampleIndex, 1:Ncells ) = ENU( 2, : )'; if twoZs == 0 Data.( [ dataModeWord '_VelUp' ] )( sampleIndex, 1:Ncells ) = ENU( 3, : )'; else Data.( [ dataModeWord '_VelUp1' ] )( sampleIndex, 1:Ncells ) = ENU( 3, : )'; Data.( [ dataModeWord '_VelUp2' ] )( sampleIndex, 1:Ncells ) = ENU( 4, : )'; end end Config.( [ configModeWord '_CoordSystem' ] ) = 'ENU'; disp( 'XYZ transformed to ENU' ) end