Entry and exit points / trajectories through a rectangle (orbits)

Hi, I have a set of (x,y) points that are continuous orbit trajectories crossing a certain rectangle at specific times. The (x,y) points are always consecutive, meaning they are not scattered points but they follow a certain path with time. At some point they will enter the rectangle and eventually will leave it through another side. This will occur several times. I assume I need to find the (x,y) points that are inside the rectangle and then compute a diff to identify the boundaries. But I need to know the 'entry' and 'exit' points, and I don't always manage. If I identify them wrongly (entry and exit are switched) then the error is very large. Any idea on how to do that?

Attached a mat file with the (x,y) points = (satlon, satlat) that define the trajectories for which I want to find entry and exit points to the rectangle. Currently I'm using this code but it does not always work ok, only sometimes. Script for using the code also attached.

Thanks!

Answer 1

[Image Analyst]

Basically if you have the x and y of the trajectory, and the x and y of the sides, then you simply need to use interp1 to find where it crosses the rectangle lines.

[Albert Zurita]

But I need to identify if an intersection is an 'entry' or an 'exit' point.

[Image Analyst]

You can use inpolygon to identify which points in your x and y vectors are inside the rectangle. When it goes from outside to inside, those are the two points on either side of the rectangle wall, and the interpolated point in between them would be an entry point. Similar for an exit point where it goes from inside to outside.

[Albert Zurita]

In the end I want to find two vectors, an entry_index, and an exit_index. I don't mind if the entry and exit exactly intersect the rectangle, but at least they are entry and exit points to it. They are the indices of the original set of (x,y) points that define the key entry/exit points

[Image Analyst]

I guess I'm not sure why you think inpolygon won't give that to you. Did you even try it?

[Albert Zurita]

Yes yes, I have updated the original question with the sample points and the code I'm currently using considering inpolygon

[Image Analyst]

Make it easy for me to help you in a few minutes. Give the code to read in that mat file into your variables. I'll check back later.

[Albert Zurita]

no problem, I have updated the mat file and provided a m script that directy reads the file, computes the entry-exit points and plots two cases, one where it does fine and another where it doesn't.thanks!

Answer 2

[Star Strider]

This should be a relatively straightforward interpolation problem in that the only calculations necessary should be to first find the approximate indices of the respective latitude and longitude crossings and then use those to interpolate to get more precise values for the crossings. That part works well enough and produces what may be appropriate times (I created a time vector for this, use the correct one). However to check the calculations, I also interpolated the crossing longitudes for the target latitudes and the crossing latitudes for the crossing longitudes. That did not go as well as I would have liked, and they are not plotting correctly.

Getting the correct latitudes and longitudes for the interpolatioon also required checking to avoid ‘end effects’ where the coordinates abruptly reversed sign, for example from -179 to +179, since my code interpreted that as crossing the target latitudes or longitudes (correct, however returning spurious data), so a correction for that had to be included.

The problems may be in my understanding of your data and of the coordinate system (I am not used to dealing with geographical coordinates) rather than an inherent problem in my code itself, so I’m posting this anyway.

type('entry_exit.m')

cycle = load('full_test.mat')   ;

in = inpolygon(cycle.satlon,cycle.satlat,[70 70 90 90],[-20 10 10 -20]);
in_ix = find(in==1);
din_ix = find(gradient(in_ix)>5);

din_on = [in_ix(1);in_ix(din_ix(2:2:end))];
din_off = [in_ix(din_ix(1:2:end));in_ix(end)];

figure;
plot(cycle.satlon(in),cycle.satlat(in),'.');hold on;
k=25;
plot(cycle.satlon(din_on(k)),cycle.satlat(din_on(k)),'r*');hold on;
plot(cycle.satlon(din_off(k)),cycle.satlat(din_off(k)),'g*');hold on;
title('OK');

figure;
plot(cycle.satlon(in),cycle.satlat(in),'.');hold on;
k=32;
plot(cycle.satlon(din_on(k)),cycle.satlat(din_on(k)),'r*');hold on;
plot(cycle.satlon(din_off(k)),cycle.satlat(din_off(k)),'g*');hold on;
title('KO');

cycle = load('full_test.mat')

cycle = struct with fields:
    satlat: [233280×1 double]
    satlon: [233280×1 double]

slat = cycle.satlat

slat = 233280×1

81.4180 81.3900 81.3214 81.2131 81.0665 80.8836 80.6665 80.4174 80.1388 79.8332

slon = cycle.satlon

slon = 233280×1

-0.4600 -4.4854 -8.4667 -12.3692 -16.1623 -19.8207 -23.3252 -26.6624 -29.8247 -32.8093

L = numel(slat);                                                                % Length Of Data Array
Fs = 2;                                                                         % Sampling Events / Time Unit
tv = linspace(0, numel(slat), numel(slat)).'/Fs;                                % Sampling Times

latlim = [-20 10];                                                              % Get Times & Longitudes For Each Latitude

for  k = 1:numel(latlim)
    zx = find(diff(sign(slat - latlim(k))));
    latzx{k} = zx(abs(slat(zx) - latlim(k)) < 5);
end

latzx

latzx = 1×2 cell array
    {768×1 double}    {769×1 double}

tlatx = NaN(numel(latlim), max(numel(latzx{:})));
lonlatx = tlatx;

for k1 = 1:numel(latzx)
    for k2 = 1:numel(latzx{k1})
        idxrng = max(1, latzx{k1}(k2)-1) : min(latzx{k1}(k2)+1,L);
        % CheckData = [slat(idxrng) slon(idxrng) tv(idxrng)]
        tlatx(k1,k2) = interp1(slat(idxrng), tv(idxrng), latlim(k1));           % Crossing Time For This Latitude
        lonlatx(k1,k2) = interp1(slat(idxrng), slon(idxrng), latlim(k1));       % Crossing Longitude For This Latitude
        latlatx(k1,k2) = interp1(slon(idxrng), slat(idxrng), lonlatx(k1,k2));   % Crossing Latitude For This Latitude
    end
end
% latlatx.'
Latitude = table(tlatx(1,:).',lonlatx(1,:).',tlatx(2,:).',lonlatx(2,:).',  'VariableNames',["Time"+latlim(1), "Lon"+latlim(1), "Time"+latlim(2), "Lon"+latlim(2)])

Latitude = 769×4 table
    Time-20     Lon-20      Time10     Lon10 
    _______    _________    ______    _______

    92.615        -100.6    67.019    -93.791
    210.73        75.856    236.33     69.049
    396.15       -125.89    370.55    -119.09
    514.27        50.561    539.86     43.754
    699.68       -151.19    674.08    -144.38
    817.79        25.268    843.39      18.46
    1003.2       -176.48     977.6    -169.67
    1121.3     -0.025572    1146.9    -6.8341
    1306.7        158.22    1281.1     165.03
    1424.8        -25.32    1450.4    -32.128
    1610.2        132.93    1584.7     139.74
    1728.4       -50.615      1754    -57.423
    1913.8        107.63    1888.2     114.44
    2031.9        -75.91    2057.5    -82.716
    2217.3        82.342    2191.7     89.149
    2335.4        -101.2      2361    -108.01

lonlim = [70 90];                                                               % Get Times & Latitudes For Each Longitude

for  k = 1:numel(lonlim)
    zx = find(diff(sign(slon - lonlim(k))));
    lonzx{k} = zx(abs(slon(zx) - lonlim(k)) < 5);
end

lonzx

lonzx = 1×2 cell array
    {411×1 double}    {411×1 double}

tlonx = NaN(numel(lonzx), max(numel(lonzx{:})));
latlonx = tlonx;

for k1 = 1:numel(lonzx)
    for k2 = 1:numel(lonzx{k1})
        idxrng = max(1, lonzx{k1}(k2)-1) : min(lonzx{k1}(k2)+1,L);
        % CheckData = [slat(idxrng) slon(idxrng) tv(idxrng)]
        tlonx(k1,k2) = interp1(slon(idxrng), tv(idxrng), lonlim(k1));           % Crossing Time For This Longitude
        latlonx(k1,k2) = interp1(slon(idxrng), slat(idxrng), lonlim(k1));       % Crossing Latitude For This Longitude
        lonlonx(k1,k2) = interp1(slat(idxrng), slon(idxrng), latlonx(k1,k2));   % Crossing Longitude For This Longitude
    end
end
% lonlonx.'
Longitude = table(tlonx(1,:).',latlonx(1,:).',tlonx(2,:).',latlonx(2,:).',  'VariableNames',["Time"+lonlim(1), "Lat"+lonlim(1), "Time"+lonlim(2), "Lat"+lonlim(2)])

Longitude = 411×4 table
    Time70     Lat70     Time90     Lat90 
    ______    _______    ______    _______

    232.73     5.7719     177.1    -58.905
    474.77    -65.451    464.24    -76.344
    766.29    -77.646    762.38    -80.402
    1065.1    -80.779    1062.5     -81.41
    1365.3    -81.394    1362.7    -80.661
    1665.4    -80.248    1661.3     -77.24
    1963.2     -75.81    1951.5    -63.522
    2248.3     -55.95    2188.5     13.723
    2474.4     34.418    2444.1     69.199
    2744.6     72.419    2738.1     78.501
    3040.5     79.311    3037.4     81.019
    3340.2     81.224    3337.7     81.329
    3640.6     81.187    3637.6     79.807
    3940.2      79.14    3934.6     74.254
    4235.8     71.773    4213.4     46.492
    4502.9     30.114    4436.5    -47.419

figure
plot(slon, slat, '.-')
hold on
for k = 1:2
    plot(lonlatx(k,:), latlatx(k,:), 'pr', 'MarkerFaceColor','r')
    plot(lonlonx(k,:), latlonx(k,:), 'pg', 'MarkerFaceColor','g')
end
hold off
grid
% axis('equal')
xlim([70 90])
ylim([-20 10])
xlabel('Longitude')
ylabel('Latitude')

EDIT — (6 Nov 2023 at 05:24)

Revised code to plot the latitude and longitude intersections with the limits correctly. (It now works as I intend it to work.) The entry and exit times (when the actual time vector is supplied) should also now be verifiably correct. (Many of the intersections are close to the sampling points, so the sampling times nearest those intersections can be used if the exact intersection times themselves are not needed.)

.

[Albert Zurita]

Good!! I am however not that much interested in the exact intersection but the nearest point of the trajectory that is close to the edge. My code works ok but I discovered that it just doesn't work well whenever there is an isolated point of a given trajectory. This is the case for instance of the corners, where the code cannot identify an entry and a exit point because there only a single one...

[Star Strider]

Thank you!

If you only want the nearest indices to the entry and exit ponts, then use the ‘latzx’ and ‘lonzx’ values, respectively. Those should be the closest indices.