interactive HALO tracks

This notebook shows how to plot the flight tracks of the HALO aircraft. First, we’ll have only a quick look and later on we do some more work in order to get the data onto an interactive map.

Preparations

First, we’ll import pylab and the EUREC4A data catalog.

import eurec4a
cat = eurec4a.get_intake_catalog()

Inspecting the dataset

ds = cat.HALO.BAHAMAS.PositionAttitude['HALO-0122'].to_dask()
ds
<xarray.Dataset>
Dimensions:     (time: 331760)
Coordinates:
  * time        (time) datetime64[ns] 2020-01-22T14:57:36 ... 2020-01-23T00:1...
    lat         (time) float64 ...
    lon         (time) float64 ...
    alt         (time) float32 ...
Data variables:
    roll        (time) float32 ...
    pitch       (time) float32 ...
    heading     (time) float32 ...
    tas         (time) float32 ...
    trajectory  <U4 'HALO'
Attributes: (12/14)
    title:                           EUREC4A position and attitude data
    description:                     10Hz subset based on DLR BAHAMAS data, p...
    mission:                         EUREC4A
    platform:                        HALO
    instrument:                      BAHAMAS
    flight_id:                       HALO-0122
    ...                              ...
    featureType:                     trajectory
    Conventions:                     CF-1.7
    version:                         1.1
    history:                         acquired by HALO BAHAMAS, processed and ...
    file_created:                    File created by M. Klingebiel (email: ma...
    DODS_EXTRA.Unlimited_Dimension:  time

Just to have a better visual impression, we can create a quick overview plot of the data:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("./mplstyle/book")

plt.plot(ds.lon, ds.lat)
center_lat = float(ds.lat.mean())
plt.gca().set_aspect(1./np.cos(np.deg2rad(center_lat)))
plt.xlabel("longitude / °")
plt.ylabel("latitude / °")
plt.show()
_images/flight_tracks_leaflet_5_0.png

Reducing the size of the dataset

Later on, we want to plot all flights on an interactive map. Currently the dataset is rather large as the aircraft location has been recorded continuously at a high data rate. While this is good for quantitative analysis, this leads to poor interactive performance. So before going further, it is a good idea to reduce the amount of data while keeping the visual impression.

A possible idea to reduce the amount of required data is that plotting a line already does linear interpolation between two coordinate points. So if we would remove those points which are close to the linear interpolation between its neighboring points, the visual impression will stay almost the same. This idea has already been stated by Ramer Douglas and Peucker and is illustrated at Wikipedia. While the algorithm is not hard to write, it is difficult to do it efficiently in python. Thus, I’ve skipped it and use the simplification library in stead.

Note

Many algorithms for shape processing are contained in the beautiful shapely library. In general, I’d recommend using that library, but it requies the GEOS library which is a bit tricky to install.

from simplification.cutil import simplify_coords_idx
def simplify_dataset(ds, tolerance):
    indices_to_take = simplify_coords_idx(np.stack([ds.lat.values, ds.lon.values], axis=1), tolerance)
    return ds.isel(time=indices_to_take)

We can now use that algorithm to generate a simplified version of the dataset with a tolerance of \(10^{-5}\) degrees, which otherwise looks the same to the previous version.

dssimplified = simplify_dataset(ds, 1e-5)
dssimplified
<xarray.Dataset>
Dimensions:     (time: 8225)
Coordinates:
  * time        (time) datetime64[ns] 2020-01-22T14:57:36 ... 2020-01-23T00:1...
    lat         (time) float64 13.08 13.08 13.08 13.08 ... 13.07 13.07 13.07
    lon         (time) float64 -59.49 -59.48 -59.48 ... -59.52 -59.51 -59.5
    alt         (time) float32 7.656 60.04 130.9 194.4 ... 79.67 34.11 5.528
Data variables:
    roll        (time) float32 -0.2912 -2.435 -3.119 ... -1.214 2.606 1.582
    pitch       (time) float32 5.956 10.79 10.57 9.97 ... 0.4644 0.8037 2.903
    heading     (time) float32 75.55 75.75 75.53 75.02 ... 76.17 75.49 74.5 73.7
    tas         (time) float32 81.97 94.35 97.48 100.9 ... 67.07 65.64 59.13
    trajectory  <U4 'HALO'
Attributes: (12/14)
    title:                           EUREC4A position and attitude data
    description:                     10Hz subset based on DLR BAHAMAS data, p...
    mission:                         EUREC4A
    platform:                        HALO
    instrument:                      BAHAMAS
    flight_id:                       HALO-0122
    ...                              ...
    featureType:                     trajectory
    Conventions:                     CF-1.7
    version:                         1.1
    history:                         acquired by HALO BAHAMAS, processed and ...
    file_created:                    File created by M. Klingebiel (email: ma...
    DODS_EXTRA.Unlimited_Dimension:  time

We can now compare those two tracks side by side while keeping a look at the number of data points required.

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(ds.lon, ds.lat)
ax1.set_title(f"{len(ds.time)} data points")
ax2.plot(dssimplified.lon, dssimplified.lat)
ax2.set_title(f"{len(dssimplified.time)} data points")

for ax in (ax1, ax2):
    ax.set_aspect(1./np.cos(np.deg2rad(center_lat)))
    ax.set_xlabel("longitude / °")
    ax.set_ylabel("latitude / °")
    ax.label_outer()


plt.show()

ratio = len(dssimplified.time) / len(ds.time)
print(f"compression ratio: {ratio*100:.2f} %")
_images/flight_tracks_leaflet_11_0.png
compression ratio: 2.48 %

The dataset size has been substantially reduced while the visual impression stayed the same.

First interactive map

In order to show the map a little bit more interactively, we use the ipyleaflet library which creates a bridge between ipython and the Leaflet JavaScript library.

import ipyleaflet

As we will need to convert many tracks to ipyleaflet layers later on, the easiest is to create a little function for that purpose right away:

def track2layer(track, color="green", name=""):
    return ipyleaflet.Polyline(
        locations=np.stack([track.lat.values, track.lon.values], axis=1).tolist(), 
        color=color,
        fill=False,
        weight=2,
        name=name
    )

With the help of that little function, creating a map is now like a breeze:

testmap = ipyleaflet.Map(center=(13.3, -57), zoom=7)
testmap.add_layer(track2layer(dssimplified))
display(testmap)

All in one

Let’s see if we can add all the flights and provide a layer switcher so that we can have a look at all flights individually. We’ll start by loading and simplifying all out track data into a local dictionary. Requesting the datasets from the server is IO bound and thus can simply be accellerated using a ThreadPool:

from multiprocessing.pool import ThreadPool
pool = ThreadPool(20)

def get_dataset(flight_id):
    ds = cat.HALO.BAHAMAS.PositionAttitude[flight_id].to_dask()
    return flight_id, ds.load()

full_tracks = dict(pool.map(get_dataset, cat.HALO.BAHAMAS.PositionAttitude))

We still have to simplify the dataset, which is done here:

tracks = {flight_id: simplify_dataset(ds, 1e-5)
          for flight_id, ds in full_tracks.items()}

Let’s also quickly grab some colors from a matplotlib colorbar, such that we can show all tracks in individual colors:

import matplotlib
colors = [matplotlib.colors.to_hex(c)
          for c in plt.cm.inferno(np.linspace(0, 1, len(tracks)))]

We can now start with a new empty map. Let’s also have a different basemap.

m = ipyleaflet.Map(
    basemap=ipyleaflet.basemaps.Esri.NatGeoWorldMap,
    center=(13.3, -57), zoom=7
)

We’ll add all the tracks as individual layers to the map

for (flight_id, track), color in zip(tracks.items(), colors):
    m.add_layer(track2layer(track, color, flight_id))

and add a scale, a legend, layer controls and a full screen button to the map and show it. If you want to zoom in, you can for example shift-click and drag a rectangle over the area you want to zoom in more closely.

m.add_control(ipyleaflet.ScaleControl(position='bottomleft'))
m.add_control(ipyleaflet.LegendControl(dict(zip(tracks, colors)),
                                       name="Flights",
                                       position="bottomright"))
m.add_control(ipyleaflet.LayersControl(position='topright'))
m.add_control(ipyleaflet.FullScreenControl())
display(m)