extracting overlapping portions of 2 similar tracks

4 次查看（过去 30 天）
Sarvesh
Sarvesh 2024-9-12
评论： Mathieu NOE 2024-9-16
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.

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

回答（1 个）

Mathieu NOE
Mathieu NOE 2024-9-12
编辑：Mathieu NOE 2024-9-13
Ran in:
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更早的评论
Mathieu NOE
Mathieu NOE 2024-9-16
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by