From 815ff64b08291432f968bdcecbd07a1e75054bdc Mon Sep 17 00:00:00 2001
From: oboebel <oboebel@bkli03p121.dmawi.de>
Date: Wed, 11 Oct 2023 15:44:41 +0200
Subject: [PATCH] New tracking method Singel Source plus Topography tracking.
---
artoa.ini | 21 ++--
plugins/tracking/SiSoTopo.m | 242 +++++++++++++++++++++++++++++++-----
2 files changed, 226 insertions(+), 37 deletions(-)
diff --git a/artoa.ini b/artoa.ini
index acb0d3e..7f63f70 100644
--- a/artoa.ini
+++ b/artoa.ini
@@ -1,16 +1,20 @@
% ARTOA ini file for dimes first three floats
[files]
-%-soundsourcefile C:\MyData\RAFOS\HAFOS.soso
--soundsourcefile C:\MyData\RAFOS\TESTwithOFFSETS.soso
--floatfile ..\HAFOS.flo % not used in artoa4 yet
--bathyfile ..\HAFOS.prj % not used in artoa4 yet
--etopo C:\MyData\TOPO\etopo1_ice_c.flt
+-soundsourcefile C:\MyData\RAFOS\HAFOS.soso
+%-soundsourcefile C:\MyData\RAFOS\TESTwithOFFSETS.soso
+-floatfile ..\HAFOS.flo % not used in artoa4 yet
+-bathyfile ..\HAFOS.prj % not used in artoa4 yet
+-etopo C:\MyData\TOPO\etopo1_ice_c.flt
+-floatProfiles \\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat
+-trackingContours C:\MyData\ARTOA\contours_Grid5000_Step50.mat
+
[directory]
--interim C:\MyData\RAFOS\
--rfb C:\MyData\RAFOS\
--soundsource C:\MyData\RAFOS\
+-interim C:\MyData\RAFOS\
+-rfb C:\MyData\RAFOS\
+-soundsource C:\MyData\RAFOS\
+-ancillary C:\MyData\ARTOA\
% the following parameters are not yet used in artoa4
-argos ..
-rft ..
@@ -23,6 +27,7 @@
-surface
-ric C:\MyData\RAFOS\
+
[filemask]
-interim *.itm
-rfb *.rfb
diff --git a/plugins/tracking/SiSoTopo.m b/plugins/tracking/SiSoTopo.m
index 6096318..8c86d9d 100644
--- a/plugins/tracking/SiSoTopo.m
+++ b/plugins/tracking/SiSoTopo.m
@@ -1,7 +1,7 @@
-function [positions, clockError] = SiSoTopo(pData, pSoundsourcePositions, pCombinationDetails, pStartEndPosition)
+function [positions, clockError] = SiSoTopo_v2(pData, pSoundsourcePositions, pCombinationDetails, pStartEndPosition)
%MYTRACKINGMETHOD A sample function to write your own tracking algorithms
%
-% Parameters:
+% Parameters:plo
% pData (struct):
% Data structured by soundsource. The toa is already corrected
% (e.g. doppler correction etc.)
@@ -21,6 +21,8 @@ function [positions, clockError] = SiSoTopo(pData, pSoundsourcePositions, pCombi
% clockError (double):
% The estimated clock error.
+global artoaConfig;
+global artoaDataInput;
%% Get reference position and all soundsource names that are involved
referencePosition = cellfun(@str2double, strsplit(pCombinationDetails.referencePosition{1}));
soundsourceNames = strsplit(pCombinationDetails.soundsources{1});
@@ -31,40 +33,215 @@ pCombinationDetails
pSoundsourcePositions
%% Return NaN because this is a sample
-positions = NaN(size(pData.(soundsourceNames{1}).date, 1), 2);
-clockError = [];
+positions = NaN(size(pData.(soundsourceNames{1}).date, 1), 2);
+clockError = NaN(size(pData.(soundsourceNames{1}).date, 1),1);
+
+%% Additional data needed
+%%% addpath('C:\Users\oboebel\Documents\MATLAB\FloatTools\lib\utils\seawater_library') % for pressure to depth conversion
+%%% wurde nach %%%%% C:\Users\oboebel\Documents\MATLAB\artoa4\lib\vendor\seawater_library verschoben
+
+% the topography file with contour lines
+% % %
+% load('C:\Users\oboebel\Documents\MATLAB\artoa4development\tools\contours_Grid5000_Step50.mat')
+% replaced by next line
+load(artoaConfig.files.trackingContours)
+
+% get the corresponding profile depth. Profile depths are spaced by 10
+% days. To match daily rafos fixes, we use linear interpolation in-between
+% to obtain daily pressure readings. NOte: the is intended for un-grounded
+% floats. Grounded floats will actually have daily bottom pressure
+% readings. This might be queried by comparing the depth readings in the
+% RAFOS data block with the BOTTOM block pressure readings.
+% global artoaDataInput
+% profileInfo = artoaDataInput.rfb.PROFILE_DEPTH
+% profilePressure = profileInfo(:,8)
+% profileID = profileInfo(:,1)
+% now the pressure date and the distance date need to be matched.
+%--> interpolate profile pressures based on distance dates
+
+% this needs to be recoded in line with artoa4 structure, rfb structure,
+% and apex to profile data, currently hard coded to 8882
+%load('C:\Users\oboebel\Documents\MATLAB\artoa4development\floatdata.mat')
+% % % load('\\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat')
+floatProfilesAll = load(artoaConfig.files.floatProfiles);
+% added link from local nemo_data.mat to \\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat'
+
+[floatProfileDataSingle, index2Float] = nemo.float.find(floatProfilesAll.nemo_data,artoaDataInput.rfb.FLOAT.floatname);
+profileNumbers = 1:size(floatProfileDataSingle,2)
+launchposition = artoaDataInput.rfb.FLOAT.launchpos;
+% % % %%% 8882
+% % % index2Float = 4;
+% % % profileNumbers = 1:77;
+% % % launchposition = [-74.5704, -29.1891]; % plot the launch position % hardcoded, to be adjusted
+% % % %%% 8885
+% % % index2Float = 6;
+% % % profileNumbers = 1:36;
+% % % launchposition = [-74.5704, -29.1891]; % dummy, achtung, falsche werte
+% % % %%% 9223
+% % % index2Float = 22;
+% % % profileNumbers = 1:36;
+% % % launchposition = [-70, -8.]; % plot the launch position % hardcoded, to be adjusted
+cutOff = 1995; % dbar, if maxPressure less, float is considered grounded (black points), otherwise plot red points
+%%% Future: add pop-up dialog to modify parameters if wanted.
+
+
+% all data is loaded, here we go:
% control plots %%
-figure
-lati = [-78 -60]
-loni = [-70 0]
-h = axesm('MapProjection','stereo','MapLatLimit',lati,'MapLonLimit',loni,'Frame','on','Grid','on')
-MapHandle = worldmap(lati,loni);
-
-% plot the map
-load('C:\Users\oboebel\Documents\MATLAB\CPTCMAP\WeddellContours.mat')
-helper = fliplr(cArray{3}');
-latMap = helper(:,1);
-lonMap = helper(:,2);
-plotm(latMap,lonMap) % index 3 = 2000m , then upwards in 50m increments
-hold on
-
-% program uses only first selected sound source
-centerSoSo = pSoundsourcePositions.W18a
-distance = [pData.W18a.distance{:}]';
-distance = distance/60/1.852 %--> distence in arch length
-
-for i = 1:length(distance)
- [lat,lon] = scircle1(centerSoSo(1),centerSoSo(2),distance(i))
- %plotm(lat,lon,'r');
- [xi,yi] = polyxpoly(latMap,lonMap,lat,lon); % hier muss die zum Druck passende latMap lonMap verwendung finden.
- plotm(xi,yi,'xk')
- drawnow
+figure(2)
+lati = [-75 -72]; % Markus Floats --> move to popup
+loni = [-42 -20];
+
+lati = [-75 -60]; % Neumayer Floats
+loni = [-35 -5];
+h = axesm('MapProjection','stereo','MapLatLimit',lati,'MapLonLimit',loni,'Frame','on','Grid','on');
+MapHandle = worldmap(lati,loni); hold on;
+
+% % prepare the map
+% helper = fliplr(cArray{3}');
+% latMap = helper(:,1);
+% lonMap = helper(:,2);
+% hold on
+
+sourceName = pCombinationDetails.soundsources;
+sourceName = sourceName{1}
+[token,remain] = strtok(sourceName,' ');
+sourceName = token;
+title(sourceName)
+
+plotm(launchposition(1),launchposition(2),'ko') % plot the launch position % hardcoded, to be adjusted
+
+eval(['centerSoSo = pSoundsourcePositions.',sourceName])
+plotm(centerSoSo(1),centerSoSo(2),'ob','MarkerSize',10)
+eval(['distanceVar = [pData.',sourceName,'.distance{:}]']); % these are the distances from the TOA, calculated on the fly by artoa4
+distanceVar = distanceVar';
+distanceVar = distanceVar/60/1.852; %--> distence in arch length
+eval(['pongDate = [pData.',sourceName,'.date{:}]']);
+pongDate = fliplr(rd2dmy(pongDate'));
+pongDate = datetime((pongDate)); % this is the array of pong dates
+startDate = pongDate(1);
+
+for i = 1:length(profileNumbers)
+% % % profileData = nemo_data{index2Float}{i}.PROFILE_DATA;
+ profileData = floatProfileDataSingle{i}.PROFILE_DATA;
+ pressureMax(i,1) = profileData(1,7);
+ depthMax(i,1) = sw_dpth(pressureMax(i,1),centerSoSo(1));
+ % temperatureMax(i,1) = profileData(1,5);
+ profileTech = floatProfileDataSingle{i}.PROFILE_TECHNICAL_DATA;
+ pressureMaxTime(i,1) = datetime(profileTech.xmit_ascent_start_time); % this is when the pofile started
+ profileNumber(i,1) = profileTech.xmit_profile_number;
+end
+
+if ~pCombinationDetails.backwardTracking;
+ currentRefPosition = pCombinationDetails.referencePosition;
+elseif pCombinationDetails.backwardTracking
+ currentRefPosition = pCombinationDetails.referencePositionEnd;
+end
+currentRefPosition = cell2mat(currentRefPosition);
+[token,remain] = strtok(currentRefPosition,' ');
+currentRefPosition = [str2num(token), str2num(remain)]
+hp_select = plotm(currentRefPosition(1),currentRefPosition(2),'k*');
+textm(currentRefPosition(1),currentRefPosition(2),' refPos');
+
+x = datenum(pressureMaxTime) ; % the profile timestamps
+v = depthMax; % the profile depths at the profile timestamps
+xq = datenum(pongDate); % the RAFOS timestamps
+vq = interp1(x,v,xq); % the interpolated profile depths at the RAFOS timestamps
+pongDepth = vq;
+
+%figure(3) % pongDepth only covers the period of this sound source
+x = datenum(pressureMaxTime) ; % the profile timestamps
+v = pressureMax; % the profile depths at the profile timestamps
+xq = datenum(pongDate); % the RAFOS timestamps
+vq = interp1(x,v,xq); % the interpolated profile depths at the RAFOS timestamps
+% plot(x,v,'o'); hold on
+% plot(xq,vq,'-')
+% ylabel('pressure [dbar]')
+% plot( get(gca,'Xlim')', [cutOff cutOff]')
+
+for i = 1:length(distanceVar) % loop over all distances / received pongs
+ [lat,lon] = scircle1(centerSoSo(1),centerSoSo(2),distanceVar(i));
+ figure(2)
+% plotm(lat,lon,'r');
+ % find the contour fitting to this depth - latMap, lon Map
+ currentDepth = pongDepth(i);
+ currentDate = pongDate(i);
+
+ if ~isnan(currentDepth)
+% currentDepth = sw_dpth(currentDepth,centerSoSo(1));
+
+ helper = abs(abs(contourLevels)-currentDepth);
+ index2contour = find( min(abs(helper)) == helper) ; % this is the contour closes to the current position
+ helper = fliplr(cArray{index2contour}');
+ latMap = helper(:,1);
+ lonMap = helper(:,2);
+
+ index = (latMap > -68) | (latMap < -74.6) | (lonMap < -60);
+ latMap(index) = NaN;
+ lonMap(index) = NaN;
+
+ figure(2)
+ plotm(latMap(2:end),lonMap(2:end),'k:','MarkerSize',1) % plot the contour line
+
+ [xi,yi] = polyxpoly(latMap,lonMap,lat,lon);
+% P1 = polyshape(latMap,lonMap)
+% P2 = polyshape(lat,lon)
+% C = intersect(P1,P2)
+% xi = C(1)
+% yi = C(2)
+
+ if ~isempty(xi) % continue/plot only if intersection exists
+
+ hp_temp = plotm(xi,yi,'.','MarkerSize',10,'Color',[.6 .6 .6]); % all potential positions in grey
+ [ARCLEN, AZ] = distance(currentRefPosition(1),currentRefPosition(2),xi,yi); % chose nearest Neighbor
+ indexNearby = find(ARCLEN == min(ARCLEN)) ;
+ xi = xi(indexNearby);
+ yi = yi(indexNearby);
+ % currentRefPosition = [xi,yi];
+ hp_select = plotm(xi,yi,'.','MarkerSize',10,'Color',[.6 .6 .6]); % overplot the closest position
+
+ [i length(distanceVar) currentDepth xi yi];
+ if currentDepth > sw_dpth(cutOff,centerSoSo(1)) % float touches bottom; cutoff[dbar]=1995, cutoff[m]=1965
+ set(hp_select,'color','r'); % float did touch bottom.
+ else
+ if ~pCombinationDetails.backwardTracking; % ergo forward tracking
+ % if currentDate < (startDate+360) % parking depth nominally 2000m
+ if currentDate < (startDate) % parking depth nominally 2000m
+ set(hp_select,'color','k');
+ else % parking depth nominally 800m
+ set(hp_select,'color','g');
+ end
+ elseif pCombinationDetails.backwardTracking; % ergo backward tracking
+ if currentDate >= (startDate) %
+ set(hp_select,'color','g');
+ else % parking depth nominally 800m % der Teil hier ist noch unlogisch und muss überarbeitet werden
+ set(hp_select,'color','g');
+ end
+ end
+ end
+ positions(i,:) = [xi, yi]; % see line 34
+ clockError(i,:) = 0;
+
+ delete(hp_temp)
+
+ % check if profile label shall be added
+ helper = find(abs(datenum(currentDate)-x) <= 2 ); %0.75 current Pong date is within a day of profile date
+ if ~isempty(helper)
+ textm(xi,yi,num2str(profileNumber(helper)),'HorizontalAlignment','left','VerticalAlignment','bottom')
+ end
+
+ end
+ drawnow
+ pause
+ end
end
+if 1 == 2
+ exportgraphics(gcf,'C:\MyData\APEX\8882.pdf','Append',true)
+end
@@ -72,3 +249,10 @@ end
end
+
+
+
+
+
+
+
--
GitLab