Borrar filtros
Borrar filtros

extracting overlapping portions of 2 similar tracks

52 visualizaciones (últimos 30 días)
Sarvesh
Sarvesh el 12 de Sept. de 2024 a las 12:06
Comentada: Mathieu NOE el 16 de Sept. de 2024 a las 8:53
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.

Respuestas (1)

Mathieu NOE
Mathieu NOE el 12 de Sept. de 2024 a las 12:51
Editada: Mathieu NOE el 13 de Sept. de 2024 a las 6:59
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 comentarios
Sarvesh
Sarvesh el 14 de Sept. de 2024 a las 0:34
@Mathieu NOE yes i think interpolation should be fine.
Mathieu NOE
Mathieu NOE el 16 de Sept. de 2024 a las 8:53
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

Iniciar sesión para comentar.

Categorías

Más información sobre Line Plots en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by