Clipping trajectories to eurofirs

[obbe79]

Hello,
I would like to clip so6 trajectories to specific FIRs.
For instance, I have collected all the trajectories crossing a specific FIR:

bdx_so6_LIMM = so6_m1.inside_bbox(eurofirs['LIMM'])

and now I would like to clip them.
Unfortunately the clip command does not work with the eurofirs objects (I don't have AIRAC files):

bdx_so6_LIMM['BCS361'].clip(eurofirs['LIMM'].flatten())

I get the following error:
NotImplementedError: Multi-part geometries do not provide a coordinate sequence

I'm probably missing something. Is there a better way to do what I'm trying to do?

Thank you,
Rob

[xoolive]

Thank you for your report @obbe79

I was able to reproduce the bug with the so6 example on the repository

so6_LIMM = so6.inside_bbox(eurofirs['LIMM'])
so6_LIMM['BAW33CE'].clip(eurofirs['LIMM'].flatten())

I think it is more a problem on the SO6 side than on the FIR side. Let me investigate further when I have more time (not this week); I have fixed a similar issue with regular ADS-B flights and considered a different way to address the question. I can probably correct the SO6 in a similar manner.

[obbe79]

Hi,
thank you for your prompt answer.
I noticed that the inside_bbox trajectories contain also trajectories that do not cross the FIR but cross the bounding box.
So I tried to check if the problem were still present if I extracted the intersecting trajectories:

inter_so6_LIMM = bdx_so6_LIMM.intersects(eurofirs['LIMM'])

Unfortunately I still get the abovementioned error on some trajectories.
However I noticed that the error appears when processing a trajectory that enters and exits multiple times the area due to the jagged boundary of the LIMM FIR.

(the red line is the trajectory and the blue line is the FIR boundary)
Could this be the problem?
Thank you,
Rob

[xoolive]

This is indeed the problem I suspected :)

I will fix it the same way it is done now with Flight classes, i.e. you will get the trajectory between the first point of entry and the last point of exit of the trajectory. (In practice, this is what makes the most sense from an ATC point of view)

What you noticed with inside_bbox is correct according to the specification. The thing is that since this first check is very fast to compute, it makes sense to apply it before intersects (and not try to check for an intersection if it doesn't cross the bounding box)

[obbe79]

Is it possible to join multiple intersections? For example, I would like to extract all the trajectories that go through the LIMM FIR and LIRR FIR. Is there an easy way to do this?

I was able to do it creating a polygon that joins the italian FIRs, and then creating a new ExtrudedPolygon() object from which I created a new Airspace() object.

import traffic.core.airspace as asp

FIR_IT = cascaded_union([eurofirs['LIMM'].flatten(),eurofirs['LIRR'].flatten(),eurofirs['LIBB'].flatten()]) 
fir_sector = asp.ExtrudedPolygon(FIR_IT, lower=0, upper=40000) 
FIR_IT_airspace = asp.Airspace('ita_fir',[fir_sector])

I'm not sure how lower and upper must be specified (is it an altitude? in feet? meters?)

Thank you.

[xoolive]

lower and upper are altitudes in flight levels (100*ft). The airspace is an extruded polygon between those two altitudes. For FIR it is usually not very complex: GND (ground) to UNL (unlimited) or FIR covers the lower space and UIR the upper space.

Note to myself, implement SO6.__add__ and so6.Flight.__add__ s.t. this is possible (as it is with ADSB data):

sum(f for f in so6 if f.intersects(eurofirs['LIMM']) and f.intersects(eurofirs['LIRR']))

[xoolive]

@obbe79 I started to have a look but my faulty flight actually did not intersect the sector (so clip would fail, it will return None in the corrected version).

Could you please try bdx_so6_LIMM['BCS361'].to_pickle('bcs361.pkl') and send the pickled file to my personal email so that I can be sure I debug this particular case?

Thank you for your support

[xoolive]

@obbe79 also please be aware that most flights will actually fly above LIMM FIR. (stops at FL295, but UIR covers FL295 and above) Current clip method only looks at the ground footprint of the airspace while intersects look at the altitude extent of the airspaces.

Now that I look at it again, I think a future version of clip could "easily" work on airspaces as well, but that's not the case yet!

By the way, I forgot I had this, but you can do eurofirs['LIMM'] + eurofirs['LIRR']? 😄

[obbe79]

Thank you for your answer.

Yes, summing the FIRs works! 😄

I also sent you the trajectory you asked for.
As for UIRs, is there a different dict for those?

[xoolive]

@obbe79 Could you please try now?

Also be careful if you use pickled version of your so6, I made a change in the library: timestamps are now explicitly marked with UTC timezone. This may cause issues if you cached data before... (just rewrite your cached files)

If you parse directly so6 files every time you run your code (slower), you will not see any difference.

About your specific flight:

import matplotlib.pyplot as plt

from traffic.data import SO6, eurofirs
from traffic.drawing import EuroPP, countries

ex = SO6.from_file('bcs361.pkl')

with plt.style.context('traffic'):
    
    fig, ax = plt.subplots(subplot_kw=dict(projection=EuroPP()))
    
    ax.add_feature(countries())
    
    ex['BCS361'].plot(ax)
    ex['BCS361'].clip(eurofirs['LIMM'].flatten()).plot(ax)
    eurofirs['LIMM'].plot(ax)
    
    ax.set_extent(eurofirs['LIMM'])