extracting overlapping portions of 2 similar tracks

hi and thanks in advance

I have cyclone tracks extracted from 2 different sources labeled as trackBom and trackRe (stored in sampleTracks.mat). I am trying to bias-correct some variables along these tracks. To do that, I need to extract the overlapping portions of the tracks. The issue is that these tracks do not exactly align in space or time. I have attached a figure showing the original tracks and the overlapping portions that I want to extract. What i want is to extract the overlapping (or almost overlapping) portions of the tracks as shown in the picture and store them as trackBomLap and trackReLap to impliment my bias correction metods.

your help is much appreciated.

Answer 1

[Mathieu NOE]

hello

you could do something like this ... I have manually checked which track is the longest and the code works here because trackBom is shorter than trackRe. If another situation arise , the logic have probably to be revisited according to which one is the shortes (easy) - compare lR1 with lB1

tried manually on the 4 datasets (ok), but the code is shown below for the first data set - also easy here to adapt to the other data sets.

code :

load('sampleTracks.mat')

%   Name            Size            Bytes  Class        Attributes
% 
%   trackBom1      43x2              1945  timetable              
%   trackBom2      40x2              1897  timetable              
%   trackBom3      24x2              1641  timetable              
%   trackBom4      17x2              1529  timetable              
%   trackRe1       59x2              2201  timetable              
%   trackRe2       91x2              2713  timetable              
%   trackRe3       27x2              1689  timetable              
%   trackRe4       37x2              1849  timetable 

% trackRe1
[xR1,yR1,lR1] = getXYfromTT(trackRe1);

% trackBom1
[xB1,yB1,lB1] = getXYfromTT(trackBom1);

% find nearest point (start zone)
[d1,i1] = min((xR1 - xB1(1)).^2+(yR1 - yB1(1)).^2);

% find nearest point (end zone)
[d2,i2] = min((xR1 - xB1(end)).^2+(yR1 - yB1(end)).^2);

% extract portion of segment of trackRe
xR1e = xR1(i1:i2);
yR1e = yR1(i1:i2);

plot(xR1,yR1,'k',xB1,yB1,'r',xR1e,yR1e,'c',xR1(i1),yR1(i1),'db',xR1(i2),yR1(i2),'dg','markersize',15);
xlabel('lon');
ylabel('lat');
legend('trackRe','trackBom','trackRe extract', 'Location', 'south')

%%%%%%%%%%%%%%%%%%
function [x,y,ll] = getXYfromTT(TT)
    x = TT.lon;
    y = TT.lat;
    ll = sum(sqrt(diff(x).^2 + diff(y).^2)); % length of track
end

[Image Analyst]

Or you could place the legend better using the 'Location' input:

legend('trackRe','trackBom','trackRe extract', 'Location', 'south')

[Mathieu NOE]

@Image Analyst you're right, on my PC the legend was not hidding some part of the plot

will update my answer accordingly

have a great day

[Sarvesh]

@Mathieu NOE. Thank you for your response and the code. I tried implimenting the code on my data and I have an issue (which i know has very simple solution but i am unable to figure out) that i was hoping you could help me with. I want to output 2 arrays that are the overlapping portions from both tracks. atm the code outputs xR1e and yR1e which are from RE data only. How can I output the x and y from BOM data corrosponding to the overlaps?

Also, manualling checking 200+ tracks is going to take some time ( i have additional tracks from other sources which would mean close to 1000+ tracks), is there anyway to automate the track length checking part and extract the track based on this automation?

PS I am not very fluent in coding

Update: I have figured out how to extrack the data from both Re and Bom tracks. but i am still unable automate the code.

[Mathieu NOE]

hello again @Sarvesh

automate the code shouldn't be too difficult , but it depends how all your data are stored (like all data is in one bid file or , as I would prefer, one file contains one data (or Bom and Re data but only for one record) .... please clarify and if possible share some sample files

also , if you have updated my code , pls share the modified version so I am on the same page as you

all the best

[Sarvesh]

Hi @Mathieu NOE,

I have attached the data and the code (along with a fuction that is required fo the main code) that i modified. I have commented my issues and some clarifications in the beginning of the code (might be a bit messy but i tried my best).

much appreciated.

[Mathieu NOE]

I try to find a way to tackle this problem :

i am only extracting over lapping tracks of bom and re if and only if their heights are same. This is because when i will correct the bias, i need the sizes of both final wind data to be same. I need to better code this portion so that i can use more tracks to correct the bias. any suggestions?

the way we proceed today does not ensure that both lapping tracks ( bom and re) have same dimensions, as the spatial increment may be quite different between the two data sets

somehow , if you really need the results to be of same size we need to interpolate both tables to the same length, so again use retime

is it what you want to do ?

[Sarvesh]

@Mathieu NOE yes i think interpolation should be fine.

[Mathieu NOE]

hello again

so the idea here is we will change the sampling rate of the shortest timetable to match the dimension of the taller one

this is done in these lines :

    % to get same timetable height, resample the shortest timetable
    
    if height(tempBomLap) < height(tempReLap) % resample tempBomLap
        dt = 1/(6*60*60)*height(tempReLap)/height(tempBomLap); % the term 1/(6*60*60) is related to your 6 hours interval resampling (above in your code)
        tempBomLap = retime(tempBomLap,'regular','linear','SampleRate',dt);
 
    else  % resample tempReLap
        dt = 1/(6*60*60)*height(tempBomLap)/height(tempReLap); % the term 1/(6*60*60) is related to your 6 hours interval resampling (above in your code)
        tempReLap = retime(tempReLap,'regular','linear','SampleRate',dt);
    end
 

full code :

    %% 
ReBarraMerged = xReData;
observed_bom = xBomData;

%initialize final table variable
bomLap = [];
reLap = [];

for a = 1:height(dataMatch)
    % Re Track
    trackRe = ismemberx(dataMatch.idBarra(a), ReBarraMerged.id, ReBarraMerged); % this is a function i created. supplied
    time = trackRe.time;
    lon = trackRe.lon;
    lat = trackRe.lat;
    wind = trackRe.wind;
    trackRe = timetable(time, lon, lat, wind);
    trackRe = retime(trackRe, "regular", "linear", "TimeStep", hours(6));
    
    % Obs Track
    trackBom = ismemberx(dataMatch.idBom(a), observed_bom.id, observed_bom);
    time = trackBom.time;
    lon = trackBom.lon;
    lat = trackBom.lat;
    wind = double(trackBom.wind);
    trackBom = timetable(time, lon, lat, wind);
    trackBom = retime(trackBom, "regular", "linear", "TimeStep", hours(6));
    
    % find the shortest track between Re and Bom. compare lR1 and lR2
    % trackRe1
    [xRe1,yRe1,lengthRe1] = getXYfromTT(trackRe);
    % trackBom1
    [xBom1,yBom1,lengthBom1] = getXYfromTT(trackBom);
    
    % find the overlaping track between Bom and RE data. RE data will be output
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % find nearest start point
    [~,idx1] = min((xRe1 - xBom1(1)).^2+(yRe1 - yBom1(1)).^2);
    % find nearest end point
    [~,idx2] = min((xRe1 - xBom1(end)).^2+(yRe1 - yBom1(end)).^2);
    
    % extract portion of segment of trackRe
    xReLap = xRe1(idx1:idx2);
    yReLap = yRe1(idx1:idx2);
    tempReLap = trackRe(idx1:idx2, :);
    %tempReLap = table(time, xReLap, yReLap, 'VariableNames',{'time', 'lon', 'lat'});
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % for overlaping track between BOM and overlapping data. BOM data will be
    % output
    %%%YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
    % find nearest start point
    [~,idx3] = min((xBom1 - xReLap(1)).^2+(yBom1 - yReLap(1)).^2);
    % find nearest end point
    [~,idx4] = min((xBom1 - xReLap(end)).^2+(yBom1 - yReLap(end)).^2);
    
    xBomLap = xBom1(idx3:idx4);
    yBomLap = yBom1(idx3:idx4);
    tempBomLap = trackBom(idx3:idx4, :);
    %tempBomLap = table(time, xBomLap, yBomLap, 'VariableNames',{'time', 'lon', 'lat'});
    %%%YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
    
    % concate lapping data if height(tempBomLap) = height(tempReLap).
    % need to find better way to match the tracks so that final sizes of
    % both array are equal for bias correction
    %%%%%% ************

    % to get same timetable height, resample the shortest timetable
    
    if height(tempBomLap) < height(tempReLap) % resample tempBomLap
        dt = 1/(6*60*60)*height(tempReLap)/height(tempBomLap); % the term 1/(6*60*60) is related to your 6 hours interval resampling (above in your code)
        tempBomLap = retime(tempBomLap,'regular','linear','SampleRate',dt);
 
    else  % resample tempReLap
        dt = 1/(6*60*60)*height(tempBomLap)/height(tempReLap); % the term 1/(6*60*60) is related to your 6 hours interval resampling (above in your code)
        tempReLap = retime(tempReLap,'regular','linear','SampleRate',dt);
    end
        
      
    if height(tempBomLap) == height(tempReLap)
        disp(a)
        bomLap = [bomLap; tempBomLap];
        reLap = [reLap; tempReLap];

        f = figure;
        plot(xRe1,yRe1,'k',xBom1,yBom1,'r',xReLap,yReLap,'c',xRe1(idx1),yRe1(idx1),'dk',xRe1(idx2),yRe1(idx2),'dk','markersize',15);
        hold on
        plot(xBomLap, yBomLap, 'b', "Marker","+")
        xlabel('lon');
        ylabel('lat');
        %lgd = legend('trackRe','trackBom','trackRe extract')
        hold off
    end
    %%%%%% ************
    
end

[windObsKde, xValObs] = ksdensity(bomLap.wind, 'BoundaryCorrection','reflection');  % for OBS
[windReKde, xValRe] = ksdensity(reLap.wind, 'BoundaryCorrection','reflection');  % for RE

% figure;
% plot(xValObs, windObsKde, 'r', 'LineWidth', 2'); % plot Obs kern
% hold on
% plot(xValRe, windReKde, 'k', 'LineWidth', 2);    % plot RE kern 
% 
% lgd = legend('ATCD', 'Re');
% lgd.Box = 'off';
% lgd.Location = "northeast";
% lgd.Orientation = "vertical";
% 
% xlabel('Wind Speed (m s^-^1)');
% ylabel('Prob Density');
% grid on;
% box on
% hold off;

% clear unwanted vars to keep workspace clean
clear xRe1 yRe1 lengthRe1 xBom1 yBom1 lengthBom1 d1 idx1 d2 idx2 xReLap yReLap a...
       lgd trackRe time lon lat trackBom f time wind idx3 idx4 xBomLap...
       yBomLap tempBomLap tempReLap windObsKde xValObs windReKde xValRe