trajectories

[espinielli]

Hi,
I am looking for a decent way to deal with 4D (3D + timestamp) / 7D (4D + speed vector) trajectory data in R since long.
I am super intrigued by the last section in the README where the 4D trajectory is mentioned:

The sphier system explicitly accommodates representations of temporal
dynamics. As the simplest level, a vector trajectory can be represented through
extension of the vertex component to include a t column. The maximum that
the sphier vertex table can contain is the four columns of (x, y, z,
t).

Could you elaborate more on it?

I have an intuition that silicate could help but terminology and R expertise fail me.
To give you a little bit of context my use case is conversion of ADS-B data,
i.e. position reports broadcast by aircraft, to trajectories and then analysis of trajectories
(intersection with airspaces, bundling in flows,...)

Any comments, ideas, clarification is welcom.

[mpadge]

Hi @espinielli, you're definitely in roughly the right place here, and @mdsumner is the man to answer your questions. This sphier repo was a place to work out some ideas on the way to silicate, which should be able to get to what you want, as your intuition rightly suggests. @mdsumner I'll be very keen to hear your answer on this, because i've not yet done any silicate trajectory stuff - what's the status? Workflow? Specific functions / tools?

[mdsumner]

I'd say that silicate can represent such a data set, and round-trip it from whatever format it's in to others - but there's really nothing built-in to work with it.

I'd be really interested to explore specific applications, why don't we get a shared example and start solving some problems?

The flight_tracks data set might give an idea of how I see silicate dealing with tracks, but I haven't pursued it since that was put in as an example data set. See how vertices exist in 4-space, 4D geometry. silicate knows about x_, y_, z_, and t_ but it only preserves those (and any others), the only thing that really does anything with them with anglr::plot3d

 flight_tracks
class       : PATH
type        : Sequential
vertices    : 18942 (4-space)
paths       : 144 (1-space)
coordinates : 18942
crs         : +proj=longlat +datum=WGS84 +no_defs
> plot(flight_tracks)
> SC0(flight_tracks)  ## convert to structural segment form 
class       : SC0
type        : Structural
vertices    : 18942 (4-space)
primitives  : 18798 (2-space)
crs         : +proj=longlat +datum=WGS84 +no_defs

anglr does something with a plot. The next step might be a plot4d() function that would animate this in time. I guess you want something like a discretization of airspace and a way to quantify time spent in 3D air-cells? Please feel free to expand on what you would want.

library(anglr)
 plot3d(flight_tracks); rgl::aspect3d(1, 1, 0.5)

[mdsumner]

oops and I see a print bug,

primitives  : 18798 (2-space)

segment primitives are 1-space because they are just lines, will fix

[mpadge]

@mdsumner Can you please transfer the issue over to the more appropriate repo of silicate, so we can continue the discussion there? I'm excited to see where this goes!

[mdsumner]

(good idea)

[mdsumner]

And my apologies @espinielli moving an issue is a pretty brutal move, I hope you don't mind! Happy to explore discussions

[espinielli]

@mdsumner and @mpadge thank you for your replies and no problem at all to move the issues where it best belong.

My use case is really a practical one: determine the portion of a flight trajectory in an airspace.
With Airspace I have one simple ground-projected polygon (no holes) with one altitude interval.
The example from parse_airspace_prisme in my WIP package {trrrj} gives you an idea:

> library(trrrj)
> library(magrittr)
> es <- system.file("extdata", "airspace_prisme_462.prisme", package = "trrrj")
> readr::read_lines(es) %>%
   parse_airspace_prisme()
Simple feature collection with 589 features and 4 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 350.6942 ymin: 173.0292 xmax: 364.0546 ymax: 187.7911
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
      unit airspace fl_min fl_max                       geometry
1  LSAZCTA  LSAZESL    125    245 POLYGON ((352.2102 187.0826...
2  LSAZCTA  LSAZESL      0    245 POLYGON ((410.0669 156.5964...
3  LSAZCTA  LSAZESL    105    245 POLYGON ((384.1205 193.2281...
4  LSAZCTA  LSAZESL      0    245 POLYGON ((402.759 165.3511,...
5  LSAZCTA  LSAZESL    125    245 POLYGON ((363.1386 181.2812...
6  LSAZCTA  LSAZWSL    115    245 POLYGON ((304.3561 184.9918...
7  LSAZCTA  LSAZWSL    125    245 POLYGON ((329.2045 178.9664...
8  LSAZCTA  LSAZWSL    125    245 POLYGON ((325.3235 145.3995...
9  LSGGTMA  LSGGSUD      0    195 POLYGON ((247.8045 80.88635...
10 LSAZCTA  LSAZESL    125    245 POLYGON ((362.0952 158.4151...

So I would be interested to find the interceptions of the trajectory with an airspace.
A trajectory then being a bunch of 4D points (in the simple case where I do not have speed vector).

I am thinking of taking segments (consecutive pairs of 4D points) and keep the ones completely internal to the airspace plus finding intersection with airspace using the bisection algorithm for the one segments that have one (and only one) of the 4D points outside.

Any thoughts, again, super welcome.

[espinielli]

I played with your code

library(anglr)
library(rgl)
library(silicate)

plot3d(flight_tracks)
rgl::aspect3d(1, 1, 0.5)
axes3d( edges=c("x--", "y--", "z") )

wit the addition of runways (maybe from osmdata) it would make it further awesome...!

But I am also/more interested in the analysis of the trajectories; once you have the parts belonging
to an airspace you can calculate the flight hours, the number of flights, bundle flows together, determine the level flight portions (which in the terminal area contribute to the fuel burn inefficiency due to not flying continuous descent operations), ...

...

[mdsumner]

I hear you, it's a good example and I hope to spend some time on it - thanks!

[mdsumner]

Hi @espinielli, do you know how to make sense of coordinates like 352.21019662192;187.082627470795?

If you have a reprex with reading an airspace and comparing it to a flight track that would be good. I'll delve into it over time, but those coordinates don't make sense atm.

[mdsumner]

Here's a version of flight tracks that puts on background image fwiw

library(anglr)
library(rgl)
library(silicate)

xy <- sc_coord(flight_tracks)  ## I happen to know this is longlat
im <- ceramic::cc_location(raster::extent(range(xy$x_), range(xy$y_)) + 0.5)
dem <- ceramic::cc_elevation(im)

mesh <- as.mesh3d(dem, image_texture = im)
tracks <- reproj::reproj(flight_tracks, 
                            raster::projection(im))

plot3d(tracks, lit = F)
plot3d(mesh, add = TRUE)
rgl::aspect3d(1, 1, 0.1)
axes3d( edges=c("x--", "y--", "z") )

(ceramic needs a mapbox key but you can replace with any raster-dem + raster-image)

[espinielli]

You found a bug in my code: I am not correctly converting from the internal format (.prisme) to a degree decimal...
Those coordinates are in minutes decimal.
Fixing in {trrrj} right now
Thanks!

[mdsumner]

nice one - I'd very much appreciate a reprex that gets some airspace polygons and some flights ready to go - but another question, how is the altitude component of the airspace/s defined?

Really nothing you can tell me "too simple", I have no idea about these data :)

[espinielli]

In aviation (for airspaces and above a certain altitude) you use flight levels, https://en.wikipedia.org/wiki/Flight_level

[mdsumner]

Ah, cool - that's exactly what I needed ;)

Tell me about the lifetimes of these, or can we treat them as essentially permanent (for our purposes)?

[espinielli]

I would start with static airspaces.
In fact an airspace can exist or disappear during the days or from one AIRAC to another.

Just to give an idea, you could have 3 ATC (air traffic control) positions controlling 3 sectors during, let's say, the days and 1 single collapse sector (made of the 3) during the night.
(Note that in this scenario anyway the 3 separate sectors and the collapsed one have different IDs so they exist all the time as separate entities only they are active in specific period of time.)

But the building blocks are the so called Elementary Sectors, which you can use to build all the rest once you know the composition definition (and the time it applies to)

[mdsumner]

Ok, reasonably complex. In one sense this is easy, define the regions and figure out which volume each track point is in (it's reasonable I think to interpolate down to a equal-time interval, forget about intersecting line segments and volumes). And the easiest way to do that is build a 3D raster (maybe abstractly with pure base logic) and find out which cell every point is in. Then you have to give lifetimes - and updated identity - to each cell whenever spaces change.

I've done a lot of stuff like this, but I don't know anything like it currently in R that makes it easiesh. And as described above it's not particularly silicate-ish, But we did a lot of this for elephant seal tracks in 4D ocean models, so I'm interested to apply that scenario to this.

There is some volumetric intersection stuff in R (in Rvcg maybe), and that might make things simpler in one sense.

[mdsumner]

You do have a massive problem with the spaces wrapping around the globe though, so unless you are easily confined to parts of the globe (it probably is true that a given airspace is easy to work with locally) then https://github.com/r-spatial/libs2 - and the stars package, might be a better place to bring this up in the long run, but I don't think you'll find any easy solutions along those lines for a good while.

[espinielli]

Re: trajectory/airspace intersection, yes I was thinking of upsampling the trajectory to my desired resolution (application dependent) and simply keep what is inside the volume.

Maybe I am optimizing prematurely, but I could be dealing with 30k trajectories (in wide Europe) per day covering few years with (at most) 1 position report per second...and it feels like it could be slow.
But surely worth it. I am especially interested in trying using it as you proposed via raster but I am a noob, so any direction as where to start is very much appreciated.

Re: date line wrapping/poles: I am mainly dealing with (wide) Europe: from Atlantic Ocean till Armenia/Georgia (longitude) and from Morocco up to Svalbard (latitude), see EUROCONTROL Member States, so it is not a problem now.
I am sure we will be bitten by it one day or another!