How to manipulate z-values for coordinate depending on the polygon they lie in?

Q&A
1

I have ~10k aircraft tracks (coordinates with altitude information) and I want to plot them in a way to show how far away they are from the lower boundary of a given airspace volume.

A quick primer on airspace structure

Airspaces are separated in controlled and uncontrolled volumes. In Europe everything above 10,000 ft altitude is controlled airspace and everything beneath it uncontrolled airspace. Simply speaking in controlled airspace there is an air traffic controller responsible for keeping aircrafts properly separated and in uncontrolled airspace this is the responsibility of the pilot. Most airports that see commercial traffic have therefore an extension of controlled airspace that reaches all the way to the ground. Its shape is similar to an inverted wedding cake:

weddingcake
weddingcake1414×2000 223 KB

What I tried so far

My idea was to use the airspace polygons to select all points of the tracks that lie inside that airspace with gmtselect and use gmtmath to subtract the lower altitude of said airspace from the altitude of the aircraft at that point.

  • If the result is negative the aircraft is below and therefore outside the airspace (see (1) in the illustration above)
  • If the result is positive but <1,000 ft then the aircraft is inside the airspace but too close to the floor (see (2) and (3) in the illustration above)

My problems

  • gmtselect only gives me a subset of my points inside the airspace. My track gets torn apart. This is not desired.
  • Needs to be done for every “cake floor” separately. It seems I can’t pass a list of polygons to gmtselect
  • Gets really tedious with complex airspace geometries at large airports

I haven’t found an elegant solution so far. How would you tackle the problem?

Thank you for your ideas and all the best,
Kristof

Example to play with

FRA-airspace
FRA-airspace2554×3283 342 KB

You can clearly see the complex airspace structure in magenta. The lines with the color gradient are aircraft tracks with their altitude color-coded.

coordinates for the runways
FRA-rwy-thr.txt (249 Bytes)

airspaces, the number after the “>” is the floor of the airspace in feet
airspace.txt (2.1 KB)

a number of sample tracks in the Frankfurt area
tracks.txt (286.6 KB)

#!/usr/bin/env bash

gmt makecpt -T0/39000 -Cseis > track.cpt
gmt convert tracks.txt -i7,6,3 -Fv -o0,1,3,4,2 > trackv.txt

gmt begin FRA-airspace
  gmt set MAP_FRAME_TYPE=plain
  gmt coast -R007:24:13.64/49:01:59.90/009:54:13.64/51:01:59.90r \
    -JL008:34:13.64/50:01:59.90/49:31:59.90/50:31:59.90/20c \
    -B+t"EDDF \057 FRA"+s"Frankfurt am Main, Germany" \
    -Ia/0.5p,lightblue -Saliceblue -W0.5p,lightblue \
    --FONT_TITLE=24p,Helvetica-Bold,black
  gmt plot airspace.txt -: -B -W5p,magenta2
  gmt plot FRA-rwy-thr.txt -W0.5p,black
  gmt plot trackv.txt -Sv1p+s -W2p+cl -Ctrack.cpt -i0,1,4,2,3
gmt end
2

Is the following worth to explore ?

  1. Create a circle
  2. Divide it in as many quadrant as needed (those are the different staircases)
  3. Define several inner circles where the closest to the center are negative values and the furthest are positive (say if your phase diagram has a 10cm radius, the lower boundary would have a radius of 5cm)

For a concrete example l, this is inspired by “mjo phase diagram”

3

I’m not sure I fully understand what you are proposing? Would you mind elaborating a bit further?

Basically what I’m trying to achieve is manipulating the z-value of a coordinate. Take position (2) in the first post. Say the aircraft altitude is 3800ft (green line) and the airspace floor is at 3500ft (magenta line). I want the “new” z-value to be the difference between the two, in this case 300ft.

4

I suggest you are trying to do vector data management and analysis for which GMT is not the best tool of choice. I make extensive use of Postgis for tasks like this, support for the DE-9IM model helps when determining vector spatial relationships, and my GMT scripts plot the data directly from a SQL query returning the data analysed as I require and formatted for GMT.

5

Hi Kristoff.

Interesting challenge. Could you add a figure of the expected result? It would help me to understand better the problem.

Thanks.

6

What I suggest is to look at the problem in a non-geographic phase portrait. All the following assume that aircrafts travel through airspace with incremental floor.

edited, see bottom

The circle would be divided in sectors (quadrant), and each sector represents an airspace. Say, **quadrant "1"** would be all airspaces with a floor at 1000ft; **quadrant "2"** all with a floor at 1500ft; *and so on* Let's assume the highest floor is 5000ft. This is the radius (+ a little more) of the circle. In quadrant 1, you draw an arc at 1500ft In quadrant 2, one at 2000ft *and so on* Now all you want is to plot the altitude of the aircraft within each quadrant as it moves toward the next (or previous) one. **Pro/Con** You don't care about the direction of the aircraft or its precise location, all that matters is the floor, which means that several "geographical sectors" can share the same "quadrant" Is that clearer ?

Edit : on second thought, a circle is not necessary at all. This can be linear (easier to plot). Left hand side is “departure airport” and right hand side is “arrival” (they can be the same). X is a pseudo distance through sectors and Y is altitude.

7

In others words, you are saying to create a new grid where it’s z value would be altitude of each airspace. Is it correct?

Once you created the grid, you could use grdtrack in your data to add a new column with the airspace altitude. Finally you will have to calculate the difference between the two columns.

8

I wonder if you need gmt spatial -N instead of gmt select.

9

Yup, I edited my post… A circle is too complicated a s the same result can be achieved on a linear frame

10

Thank you gentlemen for your ideas. Back to the drawing board. I’ll report back!