Plot sun elevations along a flight path.

sun_flight_map.R

# sun_flight_map.R
#
# Plot sun elevations along a flight path. 
#
# Airport information can be sourced from http://openflights.org/data.html.
# Note that the time zone offsets provided in their airports.dat never count
# daylight savings time.
#
# This code is for fun only, so it comes with no guarantee of accuracy and
# is not to be used for any serious purpose including making any decisions 
# at all ever.
#
# Author: Tim Raupach
# Blog post: http://www.cutflat.net/words/2014/05/19/flight-paths/

require(data.table)
require(maps)
require(geosphere)
require(ggplot2)
require(maptools)
require(scales)

flightmap = function(from, to, startTime, endTime,
                     airports, buffer=5, textSize=12, ptInterval=20) {
  # Map a flight path with an indication of sun elevations along the way.
  # 
  # Args:
  #   start, end: IATA airport codes for to/from airports.
  #   startTime, endTimes: start and end local times, in a format 
  #                        understood by as.POSIXct.
  #   airports: data.table of airport information.
  #   buffer: Number of degrees to add around the route on the map 
  #           (NA for the whole world, default: 5).
  #   textSize: Text size on the map (default: 18).
  #   ptInterval: Plot a light point every ptInterval minutes (default: 20 m).
  #
  # Returns: A ggplot2 plot.
  
  # Find airport information.
  fromAirport = airports[IATA == from]
  toAirport = airports[IATA == to]
  stopifnot(dim(fromAirport)[1] == 1 & dim(toAirport)[1] == 1 )
  
  # Convert start and end times into UTC times.
  startUTC = fromAirport[, as.POSIXct(startTime, tz="UTC") - 3600*tzOffset]
  endUTC = toAirport[, as.POSIXct(endTime, tz="UTC") - 3600*tzOffset]
  
  # Get coordinates for the start/end points.
  startCoord = fromAirport[, list(long, lat)]
  endCoord = toAirport[, list(long, lat)]
  
  # Calculate distance, flight time, and average speed.
  flightDist = distVincentyEllipsoid(startCoord, endCoord) / 1000    # [km].
  flightTime = as.numeric(difftime(endUTC, startUTC, units="hours")) # [h].
  averageSpeed = flightDist / flightTime                             # [km/h].
  
  # Get points along the flight path; one point per minute (why not?).
  minutes = flightTime * 60
  flightPath = data.table(min=seq(1, minutes), 
                          gcIntermediate(startCoord, endCoord, n=minutes))
  setnames(flightPath, "lon", "long")
  
  # Get the UTC time for each point along the journey.
  flightPath[, UTCtime := startUTC+(60*seq(1, length(lat)))]
  
  # Get the sun elevation for each point.
  flightPath[, sunElevation := solarpos(cbind(long, lat), UTCtime)[2], by=min]
  
  # Set the map limits. Use 178 to avoid problems with flight paths
  # that cross the -180/180 longitude point. 
  if(!is.na(buffer)) {
    longLimits = flightPath[, c(max(min(long)-buffer, -178), 
                                min(max(long)+buffer, 178))]
    latLimits = flightPath[, c(min(lat)-buffer, max(lat)+buffer)]
  } else {
    longLimits = c(-178, 178)
    latLimits = c(-90, 90)
  }
  
  # Draw a map OF THE WORLD!
  worldMap = fortify(map_data("world"))
  map = ggplot(data=worldMap, aes(x=long, y=lat)) + 
    geom_map(aes(map_id=region), map=worldMap, fill="grey", 
             colour="black", size=.4) + theme_bw(textSize) +
    scale_x_continuous(limits=longLimits) +
    scale_y_continuous(limits=latLimits) +
    labs(x=parse(text="Longitude~group('[',{}^o~E,']')"),
         y=parse(text="Latitude~group('[',{}^o~N,']')"),
         title=sprintf("%s (%s) to %s (%s), %.1f hours, %.0f km/h", 
                       from, startTime, to, endTime, flightTime, 
                       averageSpeed)) + coord_fixed()
    
  # Set up colours. When the sun is above 12 degrees it's day time, 
  # below -12 degrees it's night. In between it's spooky twilight.
  flightPath[, colour := "#FF0000"] 
  flightPath[sunElevation <= -12, colour := "#00008B"]
  gradient = colorRampPalette(c("darkblue", "red"))(nCols <- 100)
  flightPath[sunElevation > -12 & sunElevation < 12,
             colour := gradient[round(rescale(sunElevation, c(1, nCols)), 0)]]
  
  # Quantify how bright it will be on a scale from zero to ten.
  flightPath[, sunSize := 0]
  flightPath[sunElevation > -12, 
             sunSize := rescale(sunElevation, c(0, 10), c(-12, 90))]
  
  # What colour are the first and last points on the line?
  startColour = flightPath[order(UTCtime), colour[1]]
  endColour = flightPath[order(UTCtime, decreasing=TRUE), colour[1]]
    
  # Draw on the start and end points and the flight path.
  map = map +
    geom_path(data=flightPath, aes(group=NA, colour=colour), size=2) +
    geom_point(data=startCoord, colour=startColour, size=10) +
    geom_point(data=endCoord, colour=endColour, size=10) +
    geom_point(data=rbind(startCoord, endCoord), colour="white", size=7) +
    scale_colour_identity()

  # Add points indicating the brightness, every ptInterval minutes.
  sunPoints = flightPath[seq(ptInterval, dim(flightPath)[1]-ptInterval, 
                             by=ptInterval)][sunSize > 0]
  sunPoints[, offset := 0.1*max(sunSize)]
  map = map + 
    geom_point(data=sunPoints, aes(size=sunSize), colour="orange") +
    geom_point(data=sunPoints, aes(size=sunSize-offset), colour="white") +
    scale_size_continuous(guide=FALSE, trans="identity", range=c(1, 10))
  
  return(map) 
}

# # Example usage:
# 
# # Load airport information.
# airports = read.csv("airports.dat", header=FALSE)
# names(airports) = c("id", "name", "city", "country", "IATA", 
#                     "ICAO", "lat", "long", "alt", "tzOffset")
# airports = data.table(airports)
# 
# # Fix LAX info for daylight savings time.
# airports[IATA == "LAX", tzOffset := tzOffset + 1] 
# 
# # Plot Dubai to LAX.
# print(flightmap(from="DXB", to="LAX", 
#                 startTime="2014-05-20 8:55",
#                 endTime="2014-05-20 14:15", 
#                 airports=airports, 
#                 ptInterval=60, buffer=10))