Airplanes above Mainz

Quiet a while ago I bought a SDR-RTL USB dongle. Originally to capture the signal from my outdoor temperature sensor. However, no one had decipherd the signal for that device so far and I was not motivated enough to invest too much time in it. I almost forgot the dongle and remembered it recently, when I read about Flight radar and dump1090, which logs the ID, position and elevation signal of airplanes. Almost all airplanes transmit this unencrypted signal about themselves at 1090 MHz the ‘Automatic Dependent Surveillance-Broadcast System (ADS-B)’. And I was curious how I could log and visualize those flight tracks from Frankfurt/Main airport towards Mainz.

I’ve build an R package to process the data from dump1090 and present some function for flat and 3D visualization for it.

Get the raw data

This section is only relevant, if you intent to use your own recorded flight tracks. Otherwise you can jump down to section “Get it into R”. To capture the data, a SDR-RTL dongle is needed as well as the software dump1090. Compile it and run the following commands in a terminal shell. The first command captures the signal and converts it to a list, the second command collects the list via a network connection on port 30003 on localhost and appends it to a file.

./dump1090 --aggressive --quiet --net
nc -d 127.0.0.1 30003 >> flights_$(date +%Y%m%d).csv

 

My antenna is not the best and its position is not optimal, therefore I do not expect to capture all air traffic. However, it were a few MB per day.

Get it into R

Start R and load the package:

library(FlightControl)

 

Run `convert` with the recorded dates the path to where the csv files are (`source.dir`) and where to save the RData files (`target.dir`). This converts the raw CSV files above in much smaller RData files. However, already with a tiny bit of loss, since `convert` calls `airplanes` with its default parameters. You can override this by setting `raw=TRUE` or  changing `min.altitude`/`max.altitude` and/or `unit=”feet”`. Setting raw to TRUE returns the data as is. Changing only the min/max altitude matches the flight IDs and flight positions and converts feet to meters, if unit isn’t set to “feet”.

dates <- c(strptime("2016-10-28", format="%Y-%m-%d"), strptime("2016-10-30", format="%Y-%m-%d"))
convert(dates, "data", "data", max.altitude=3000)

 

Now import those .RData files and append them if desired by passing either one or several dates to `import`. And since several callSigns can appear a few times a day, e.g. when a plane lands and then starts again, those have to be split, which is done by `split.pause.callSign` after `pause.limit` seconds (default 900).

flights <- import(dates, source.dir = "data", verbose=FALSE)
flights <- breakCallSign(flights)

Or as done here use the built in dataset.

data(flights)
flights <- breakCallSign(flights)

 

Apply some filters to delete unrealistic coordinates and too short flight tracks (less than 10 nodes).

flights <- flightsPositionSubset(flights, west=0, east=15, south=45, north=55, length = 10)

 

 

Visualization

This is not part of the package, since it depends on a lot of other packages, which do not necessarily need to be installed. Therefore, the following code is not executed as part of the package, just listed as example code!

Maps

Below  is a function to create a map with the flight-tracks above an OpenStreetMap. The R package `OpenStreetMap` depends on rJava, which I experienced some trouble with in configuring. Therefore I pulled that dependecy out of my FlightControl package and list it here as example.

suppressPackageStartupMessages(require(ggplot2))
suppressPackageStartupMessages(require(OpenStreetMap))
suppressPackageStartupMessages(require(viridis))
suppressPackageStartupMessages(require(raster))
#' Map of recorded flight tracks
#'
#' Creates a map of flight tracks with an underlying OpenStreetMap layer.
#'
#' @param flights the flights data.frame
#' @return ggplot2 object
#' @export
#' @author Joerg Steinkamp \email{joergsteinkamp@@yahoo.de}
#' @import OpenStreetMap
#' @importFrom ggplot2 autoplot aes geom_path annotate scale_colour_gradientn xlab ylab
#' @importFrom viridis inferno
#' @importFrom raster extent
mapFlights <- function(flights) {
  lon=lat=altitude=callSign=NULL
  extent <- extent(x=c(min(flights$lon), max(flights$lon)), y=c(min(flights$lat), max(flights$lat)))
  extent@xmin <- extent@xmin - (extent@xmax-extent@xmin) * 0.05
  extent@xmax <- extent@xmax + (extent@xmax-extent@xmin) * 0.05
  extent@ymin <- extent@ymin - (extent@ymax-extent@ymin) * 0.05
  extent@ymax <- extent@ymax + (extent@ymax-extent@ymin) * 0.05

  ## retrieve the OpenSteetMap data and reproject it
  map <- openmap(c(extent@ymax, extent@xmin), c(extent@ymin, extent@xmax), type = "apple-iphoto", minNumTiles=10)
  mapLatLon <- openproj(map)

  p <- autoplot(mapLatLon)
  p <- p + geom_path(data=flights, aes(x=lon, y=lat, color=altitude, group=callSign), alpha=0.5, size=0.8)
  p <- p + scale_colour_gradientn(colours=inferno(255))
  p <- p + xlab("longitude") + ylab("latitude")
  return(p)
}

Then create a map of one day with ascending (2016-10-28) and one day with descending (2016-10-30) airplanes above Mainz with the above function.

flights.asc <- flightsDirectionSplit(subset(flights, date &amp;amp;lt; as.POSIXlt("2016-10-28 24:00:00")), direction="asc")
flights.dec <- flightsDirectionSplit(subset(flights, date &amp;amp;gt; as.POSIXlt("2016-10-30 00:00:00")), direction="dec")
map <- mapFlights(flights.asc)
map + labs(title="Ascending (2016-10-28)")
map <- mapFlights(flights.dec)
map + labs(title="Decending (2016-10-30)")

 

 

3D movie

However, those maps a a bit boring. Therefore, the package rgl is used, to create a rotating 3D movie of the flight tracks with two possibilities:

1. Flight tracks as lines in the “Atmosphere”.

For the terrain I prepared a digital elevation model (`dem`) based on the Shuttle Radar Topographic Mission (SRTM), which is needed as ground surface. If you want to use the package for another region, you need to get and prepare your own DEM data.

For the visualization a set of other libraries is necessary, which I make use of in the next function.

suppressPackageStartupMessages(require(rgl))
suppressPackageStartupMessages(require(raster))
suppressPackageStartupMessages(require(rgdal))
suppressPackageStartupMessages(require(sp))
suppressPackageStartupMessages(require(viridis))

#' 3D flight tracks
#'
#' @param flights the flights data.frame
#' @author Joerg Steinkamp \email{joergsteinkamp@@yahoo.de}
#' @export
#' @import rgl
#' @importFrom methods as
#' @importFrom raster raster projectRaster crop projection
#' @importFrom rgdal project
#' @importFrom sp CRS spTransform proj4string coordinates
#' @importFrom grDevices terrain.colors
#' @importFrom viridis inferno
flightTracks3D <- function(flights, fname="flights_3D.mp4") {
  callSign=NULL
  extent <- extent(data.frame(x=c(min(flights$lon), max(flights$lon)), y=c(min(flights$lat), max(flights$lat))))
  extent@xmin <- extent@xmin - (extent@xmax-extent@xmin)*0.05
  extent@xmax <- extent@xmax + (extent@xmax-extent@xmin)*0.05
  extent@ymin <- extent@ymin - (extent@ymax-extent@ymin)*0.05
  extent@ymax <- extent@ymax + (extent@ymax-extent@ymin)*0.05

## project to UTM zone 32U
  proj <- "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
  extent <- as(extent, "SpatialPolygons")
  proj4string(extent) <- CRS("+init=epsg:4326")
  extent < spTransform(extent, CRS(proj))
  extent <- extent(extent)

## reproject lon/lat of flights data.frame
  flights.proj <- project(as.matrix(flights[c("lon", "lat")]), proj)
  flights$lon <- flights.proj[,1]
  flights$lat <- flights.proj[,2]

## load a digital elevation model (DEM from SRTM data)
  full.srtm <- raster("srtm.tif")
  full.srtm <- projectRaster(full.srtm, crs=CRS(proj))
  dem <- crop(full.srtm, extent)

## extract the coordinates and z-values for the 3D surface from the DEM
  xy <- coordinates(dem)
  dem.x <- sort(unique(xy[,1]))
  dem.y <- sort(unique(xy[,2]), decreasing=TRUE)
  dem.z <- t(as.matrix(dem))
  dem.z[dem.z==0] = NA
  dem.z[is.na(dem.z)] = min(dem.z, na.rm=TRUE)

## color lookup table for DEM
  zlim <- range(dem.z)
  zlen <- ceiling(zlim[2] - zlim[1])

  dem.col <- terrain.colors(zlen)
  dem.col <- dem.col[ ceiling(dem.z - zlim[1]) ]

  zlim <- range(flights$altitude, na.rm=TRUE)
  zlen <- ceiling(zlim[2] - zlim[1])

  flight.alt.col <- inferno(zlen)

  open3d(windowRect=c(100, 100, 1280, 720), scale=c(1, 1, 5), FOV=10)
  rgl.bg(color = "black")
  surface3d(dem.x, dem.y, dem.z, color = dem.col, back = "cull")

## both methods (for and lapply) for all flight tracks are very slow!
  for (id in unique(flights$callSign)) {
    f = subset(flights, callSign == id)
    if (nrow(f) == 0)
      next
    lines3d(f$lon, f$lat, f$altitude, col=flight.alt.col[f$altitude], alpha=0.5, lwd=2)
  }

  rgl.viewpoint(0, -75, fov=10, zoom=0.55)

## play3d(spin3d(axis = c(0, 0, 1), rpm = 2))
  movie3d(spin3d(axis = c(0, 0, 1), rpm = 2), 30, dir=tempdir(), convert=FALSE, clean=FALSE, verbose=FALSE)
  rgl.close()
  if (length(system2("which", "ffmpeg", stdout=TRUE)) &amp;gt; 0) {
    system(paste0("ffmpeg -loglevel error -y -i ", file.path(tempdir(), "movie%03d.png"), " -s hd720 -r ntsc -vcodec libx264 -pix_fmt yuv420p ", fname))
    system(paste0("rm ", file.path(tempdir(), "movie*.png")))
  } else {
    message(paste0("'ffmpeg' not found. Files are in '", tempdir(), "'"))
  }
}

flightTracks3D(flights.asc, "ascending_flights_3D.mp4")
flightTracks3D(flights.dec, "decending_flights_3D.mp4")

 

 

 

 

2. Spatially interpolated minimum flight altitude as transparent surface.

The code for the second visualization is quiet similar to the first one. However, it needs much more input than just two days for reliable spatial interpolation.

suppressPackageStartupMessages(require(rgl))
suppressPackageStartupMessages(require(raster))
suppressPackageStartupMessages(require(rgdal))
suppressPackageStartupMessages(require(sp))
suppressPackageStartupMessages(require(viridis))

#' 3D surface plot
#'
#' @param flights the flights data.frame
#' @param fname the file name for the movie
#' @author Joerg Steinkamp \email{joergsteinkamp@@yahoo.de}
#' @export
#' @import rgl
#' @importFrom methods as
#' @importFrom raster raster projectRaster crop projection
#' @importFrom rgdal project
#' @importFrom sp CRS spTransform proj4string coordinates
#' @importFrom grDevices terrain.colors
#' @importFrom viridis inferno
flightsSurface3D &amp;lt;- function(flights, fname) {
  extent <- extent(data.frame(x=c(min(flights$lon), max(flights$lon)), y=c(min(flights$lat), max(flights$lat))))

  ## project to UTM zone 32U
  proj <- "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
  extent <- as(extent, "SpatialPolygons")
  proj4string(extent) <- CRS("+init=epsg:4326")
  extent <- spTransform(extent, CRS(proj))
  extent <- extent(extent)

## calculate a raster surface ("carpet") of the lowest flight altitude.
  min.altitude <- gridFlights(flights, min, proj.in="+init=epsg:4326", res.out=1000, proj.out=proj)
  carpet <- smoothFlightsGrid(min.altitude, "gam", k=23)
  carpet@data$altitude[is.na(min.altitude@data$altitude)] = NA
  carpet <- smoothFlightsGrid(carpet, "sph")

## load a digital elevation model (DEM from SRTM data)
  full.srtm <- raster("~/Cloud/Dropbox/WIP/FlightControl/data/srtm.tif")
  full.srtm <- projectRaster(full.srtm, crs=CRS(proj))
  dem <- crop(full.srtm, extent)

## extract the coordinates and z-values for the 3D surface from the DEM
  xy <- coordinates(dem)
  dem.x <- sort(unique(xy[,1]))
  dem.y <- sort(unique(xy[,2]), decreasing=TRUE)
  dem.z <- t(as.matrix(dem))
  dem.z[dem.z==0] = NA
  dem.z[is.na(dem.z)] = min(dem.z, na.rm=TRUE)

## color lookup table for DEM
  zlim <- range(dem.z)
  zlen <- ceiling(zlim[2] - zlim[1])

  dem.col <- terrain.colors(zlen)
  dem.col <- dem.col[ ceiling(dem.z - zlim[1]) ]

## eqivalent as above for altitude
  alt.x <- seq(carpet@bbox["x", "min"] + carpet@grid@cellsize["x"]/2, carpet@bbox["x", "max"] - carpet@grid@cellsize["x"]/2, carpet@grid@cellsize["x"])
  alt.y <- seq(carpet@bbox["y", "min"] + carpet@grid@cellsize["y"]/2, carpet@bbox["y", "max"] - carpet@grid@cellsize["y"]/2, carpet@grid@cellsize["x"])
  alt.z <- as.matrix(carpet)
  zlim <- range(alt.z, na.rm=TRUE)
zlen <- ceiling(zlim[2] - zlim[1])

  flight.alt.col <- inferno(zlen)

## flat look with "carpet"
  alt.col <- flight.alt.col[ ceiling(alt.z - zlim[1]) ]
  alt.alpha <- sqrt(max(alt.z) - alt.z) / sqrt(max(alt.z) - min(alt.z)) / 2

  open3d(windowRect=c(100, 100, 1280, 720), scale=c(1, 1, 5), FOV=10)
  rgl.bg(color = "black")
  surface3d(dem.x, dem.y, dem.z, color = dem.col, back = "cull")

  surface3d(alt.x, alt.y, alt.z, color = alt.col, alpha=0.7, back = "cull")
  rgl.viewpoint(0, -85, fov=10, zoom=0.55)
  movie3d(spin3d(axis = c(0, 0, 1), rpm = 2), 30, dir=tempdir(), convert=FALSE, clean=FALSE, verbose=FALSE)
  rgl.close()
  if (length(system2("which", "ffmpeg", stdout=TRUE)) &amp;gt; 0) {
    system(paste0("ffmpeg -loglevel error -y -i ", file.path(tempdir(), "movie%03d.png"), " -s hd720 -r ntsc -vcodec libx264 -pix_fmt yuv420p ", fname))
    system(paste0("rm ", file.path(tempdir(), "movie*.png")))
  } else {
    message(paste0("'ffmpeg' not found. Files are in '", tempdir(), "'"))
  }
}

First filter out those flights not to include in the analysis, e.g. remove those where the callSign starts with ‘DM’ (local flight tracks at low altitude).

require(plyr)
load("~/Path/to/large/file/20161015-20161215_split.RData")
flights <- subset(flights, !grepl("^DM", callSign))
flights <- flightsPositionSubset(flights, west=0, east=15, south=45, north=55, length = 10)
flights <- flightsDirectionSplit(flights)
total.dz <- ddply(flights[,c("callSign", "dz")], .(callSign), summarise, dz=sum(dz))
asc.flights <- subset(flights, callSign %in% total.dz$callSign[total.dz$dz &amp;gt;= 500])
dec.flights <- subset(flights, callSign %in% total.dz$callSign[total.dz$dz &amp;lt;= -500])

 

Further links related to tracking airplanes:

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s