Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Plot sun elevations along a flight path.
# 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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.