extracting overlapping portions of 2 similar tracks

Question

Sarvesh 2024-9-12，12:06

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152265-extracting-overlapping-portions-of-2-similar-tracks

评论： Mathieu NOE 2024-9-16，8:53

sampleTracks.mat

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2024-9-12，12:51

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152265-extracting-overlapping-portions-of-2-similar-tracks#answer_1515520

编辑：Mathieu NOE 2024-9-13，6:59

在 MATLAB Online 中打开

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Sarvesh 2024-9-14，0:34

@Mathieu NOE yes i think interpolation should be fine.

Mathieu NOE 2024-9-16，8:53

在 MATLAB Online 中打开

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