A Mode S based wind field

ADS-B data broadcasted by aircraft contains information about groundspeed and true track angle of its trajectory. When proper requests are sent by a Secondary surveillance radar, aircraft also send more information with true airspeed and heading angle in specific BDS5,0 and BDS6,0 messages (see Junzi Sun’s website)

Groundspeeds and true track angles are derived from the GNSS positions whereas true airspeed is computed with traditional onboard instruments like Pitot tubes.

Here, aircraft behave as a distributed network of moving sensors and researchers [1] [2] have been recommending methods to derive wind fields from Mode S data. The following method is a very basic approach to compute wind.

Please note I am not a meteorology specialist, so feel free to improve this page where it should.

The example above is wind averaged between 25°W and 55°E and between 32°N and 65°N, from FL200 and above on February 23th 2019, between 14:00 and 16:30 UTC.

from traffic.core import Traffic

t = Traffic.from_file("<fill here>")

t_extended = (
    traffic
    # download and decode EHS messages (DF 20/21)
    .query_ehs()
    # resample/interpolate one sample per second
    .resample("1s")
    # median filters
    .filter(altitude=23, track=53, heading=53, groundspeed=53, TAS=53)
    # wind triangle computation
    .compute_wind()
    # median filter
    .filter(wind_u=53, wind_v=53)
    # resample one sample per minute
    .resample("1 min")
    # do not use multiprocessing to avoid denial of service
    .eval(desc="preprocessing")
)

# t_extended.to_pickle("wind_backup.pkl")

The result of this computation is a set of trajectories: each aircraft yields one point per minute with a 4D-position (timestamp, latitude, longitude, altitude) and a wind vector decomposed along a zonal speed (wind_u) and a meridional speed (wind_v).

We then use ipyleaflet to display the wind field. The Velocity widget requires two 2D matrices of zonal and meridional components of the wind. The following method rounds lat/lon coordinates to the closest integer and average wind in each resulting cell.

import xarray as xr

def compute_grid(traffic: Traffic) -> xr.Dataset:

    avg = (
        traffic
        # remove NaN values, just in case
        .query("wind_u == wind_u")
        # prepare coordinates for the 4d-grid, also remove NaN in wind
        .assign(
            # round coordinates to the closest .33 latitude/longitude
            lat_=lambda df: (3 * df.latitude.round(0)) / 3,
            lon_=lambda df: (3 * df.longitude.round(0)) / 3,
            # This basic version averages on all altitudes/timeranges
            # but it is easy to use the following fields to display
            # wind fields in particular time ranges and altitude levels.
            alt_=lambda df: (3e-3 * df.altitude).round(0) / 3e-3,
            hour=lambda df: df.timestamp.dt.round("h"),
        )
        # compute the average wind
        .data[["wind_u", "wind_v", "lat_", "lon_"]]
        .groupby(["lat_", "lon_"])
        .mean()
    )

    # Unstack then fill the holes where possible (2D interpolation)
    u = avg[["wind_u"]].unstack().interpolate().values
    v = avg[["wind_v"]].unstack().interpolate().values

    return xr.Dataset(
        data_vars={
            "u_wind": xr.DataArray(u, coords=avg.index.levels),
            "v_wind": xr.DataArray(v, coords=avg.index.levels),
        }
    )

The following is a basic rendering delegated to ipyleaflet library.

from ipyleaflet import Map, basemaps
from ipyleaflet.velocity import Velocity

# t_extended = Traffic.from_file("wind_backup.pkl")

map_ = Map(
    center=(52, 15),
    zoom=4,
    interpolation="nearest",
    basemap=basemaps.CartoDB.DarkMatter,
)

wind = Velocity(
    data=compute_grid(t_extended),
    zonal_speed="u_wind",
    meridional_speed="v_wind",
    latitude_dimension="lat_",
    longitude_dimension="lon_",
    velocity_scale=0.002,
    max_velocity=150,
)

map_.add_layer(wind)

map_