Visualizing Flight Paths on a 3-D Interactive Earth Visual
Ankit Agarwal
While going through my blog on Vislauzaing Plane Crash Data , one of my friends asked me if to try and plot the Flight Routes of Crashed Planes on a 3-D view of Earth.
As interesting as it sounds, this was certainly an excellent suggestion and I started exploring web to find out some R package which could help me achieve this feat. Finally I found the package threejs which could do this easily.
Let us see how it can be done .
#
# As always, install and load the requisite packages first
#
#
#
list.of.packages <-
c("geosphere", # For spatial methods
"threejs", # threejs is used for 3-D interactive Earth Visualization
"rworldmap", # For creating earth map
"leaflet", # Leaflet for R provides functions to control and integrate Leaflet, a JavaScript library for interactive maps, within R.
"rgeos", # Provides functions for handling operations on topologies.
"raster", # For raster image
"DT", # For creating interactive tables
"ggplot2",
"sp" , # For Spatial processing of data
"ggmap", # To reverse geocode Long/Lat
"knitr", # TO enable 3-D visualization embedding in the HTML page
"rglwidget",
"rgl"
)
#
#You should be able to simply reuse the following lines of code as is
#
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#
if(length(new.packages)) suppressMessages(suppressWarnings(install.packages(new.packages)))
#
# By now we have installed the requisite packages. Time to load them .
#
lapply(list.of.packages,function(x){suppressMessages(suppressWarnings(library(x,character.only=TRUE)))})## [[1]]
## [1] "geosphere" "sp" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[2]]
## [1] "threejs" "geosphere" "sp" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "rworldmap" "threejs" "geosphere" "sp" "stats"
## [6] "graphics" "grDevices" "utils" "datasets" "methods"
## [11] "base"
##
## [[4]]
## [1] "leaflet" "rworldmap" "threejs" "geosphere" "sp"
## [6] "stats" "graphics" "grDevices" "utils" "datasets"
## [11] "methods" "base"
##
## [[5]]
## [1] "rgeos" "leaflet" "rworldmap" "threejs" "geosphere"
## [6] "sp" "stats" "graphics" "grDevices" "utils"
## [11] "datasets" "methods" "base"
##
## [[6]]
## [1] "raster" "rgeos" "leaflet" "rworldmap" "threejs"
## [6] "geosphere" "sp" "stats" "graphics" "grDevices"
## [11] "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "DT" "raster" "rgeos" "leaflet" "rworldmap"
## [6] "threejs" "geosphere" "sp" "stats" "graphics"
## [11] "grDevices" "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "ggplot2" "DT" "raster" "rgeos" "leaflet"
## [6] "rworldmap" "threejs" "geosphere" "sp" "stats"
## [11] "graphics" "grDevices" "utils" "datasets" "methods"
## [16] "base"
##
## [[9]]
## [1] "ggplot2" "DT" "raster" "rgeos" "leaflet"
## [6] "rworldmap" "threejs" "geosphere" "sp" "stats"
## [11] "graphics" "grDevices" "utils" "datasets" "methods"
## [16] "base"
##
## [[10]]
## [1] "ggmap" "ggplot2" "DT" "raster" "rgeos"
## [6] "leaflet" "rworldmap" "threejs" "geosphere" "sp"
## [11] "stats" "graphics" "grDevices" "utils" "datasets"
## [16] "methods" "base"
##
## [[11]]
## [1] "knitr" "ggmap" "ggplot2" "DT" "raster"
## [6] "rgeos" "leaflet" "rworldmap" "threejs" "geosphere"
## [11] "sp" "stats" "graphics" "grDevices" "utils"
## [16] "datasets" "methods" "base"
##
## [[12]]
## [1] "rglwidget" "knitr" "ggmap" "ggplot2" "DT"
## [6] "raster" "rgeos" "leaflet" "rworldmap" "threejs"
## [11] "geosphere" "sp" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "rgl" "rglwidget" "knitr" "ggmap" "ggplot2"
## [6] "DT" "raster" "rgeos" "leaflet" "rworldmap"
## [11] "threejs" "geosphere" "sp" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"#
# Set this property to enable 3-D visualization
#
knit_hooks$set(webgl = hook_webgl)Based on the source and destination of a flight from the flights data (Refer my BLOG for details around plane crash data), I gecoded the Source and Destination addresses using Google Geocoding API.
I will load the geocoded addresses for source and destination (limiting it to 1000 addresses to avoid cluttering of paths) and will visualize the paths and source and destination location using various spatial methods inclusing the 3-D interactive visualization of paths.
#
# Read the file
#
input_data<-read.csv("C:\\PB Backup\\R\\Great_Circles.csv", header = TRUE, stringsAsFactors = FALSE)
#
str(input_data)## 'data.frame': 999 obs. of 4 variables:
## $ SourceLong: num -0.236 5.93 1.025 -78.429 2.341 ...
## $ SourceLat : num 44 43.1 51.1 44.5 48.9 ...
## $ DestLong : num 2.34 24.06 4.48 -90.2 14.43 ...
## $ DestLat : num 48.9 -26.6 51.9 38.6 50.1 ...#
#
head(input_data)## SourceLong SourceLat DestLong DestLat
## 1 -0.235721 43.98616 2.34120 48.85693
## 2 5.930290 43.12534 24.06189 -26.61898
## 3 1.025000 51.07486 4.47846 51.92285
## 4 -78.429441 44.53110 -90.19956 38.62775
## 5 2.341200 48.85693 14.43322 50.07908
## 6 -75.610702 42.93709 2.34120 48.85693#
#
# Let us use the function datatable that creates an HTML widget to display rectangular data (a matrix or data frame) using JavaScript using package DT
#
datatable(input_data, rownames = FALSE)#
#
#Converting data to Spatial Format
While the purpose of the blog is to create a 3-D visualization of Flight routes, let us do some basic Spatial Data processing to get some understanding of Spatial Analytics in R
#
# Converting Data to spatial objects
#
source_df<-data.frame(SourceLong=input_data$SourceLong,SourceLat=input_data$SourceLat)
#
# Create object of class SpatialPoints using SpatialPoints function from sp package
#
source_sp<-SpatialPoints(source_df, proj4string=CRS("+proj=longlat"))
#
#
str(source_sp)## Formal class 'SpatialPoints' [package "sp"] with 3 slots
## ..@ coords : num [1:999, 1:2] -0.236 5.93 1.025 -78.429 2.341 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "SourceLong" "SourceLat"
## ..@ bbox : num [1:2, 1:2] -162.7 -54.8 175 69.7
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "SourceLong" "SourceLat"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat"#
head(source_sp)## class : SpatialPoints
## features : 1
## extent : -0.235721, -0.235721, 43.98616, 43.98616 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat#
# Convert to Spatial Dataframe
#
source_spdf <- SpatialPointsDataFrame(source_sp, data = source_df)
#
str(source_spdf)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 999 obs. of 2 variables:
## .. ..$ SourceLong: num [1:999] -0.236 5.93 1.025 -78.429 2.341 ...
## .. ..$ SourceLat : num [1:999] 44 43.1 51.1 44.5 48.9 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:999, 1:2] -0.236 5.93 1.025 -78.429 2.341 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "SourceLong" "SourceLat"
## ..@ bbox : num [1:2, 1:2] -162.7 -54.8 175 69.7
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "SourceLong" "SourceLat"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat"#
head(source_spdf)## SourceLong SourceLat
## 1 -0.235721 43.98616
## 2 5.930290 43.12534
## 3 1.025000 51.07486
## 4 -78.429441 44.53110
## 5 2.341200 48.85693
## 6 -75.610702 42.93709#
# Similarly for the Destination Locations
# Converting Data to spatial objects
#
dest_df<-data.frame(DestLong=input_data$DestLong,DestLat=input_data$DestLat)
#
# Create object of class SpatialPoints using SpatialPoints function from sp package
#
dest_sp<-SpatialPoints(dest_df, proj4string=CRS("+proj=longlat"))
#
# Convert to Spatial Dataframe
#
dest_spdf <- SpatialPointsDataFrame(dest_sp, data = dest_df)Calculating DIstance between the Source and Destination and the bearing between the Source and Destination
I’ll create the distance and bearing of each source location from the respective destination location using the distHaversine and bearing functions from geosphere package.
Bearing: A bearing is the angular difference away from a north or south baseline, ranging from 0° to 90°. Determining your bearing begins with a north or south baseline, whichever is closer to your direction. The angle is then measured east or west from the baseline.
#
# Create a Combined data frame for display purpose only
#
comb_df<-data.frame(input_data)
#
# Calculate distance between Source and Destination
#
comb_df$distance<-distHaversine(source_sp,dest_sp)
#
comb_df$bearing<-bearing(dest_sp,source_sp)
#
# Display the combined dataframe using interactive table
#
datatable(comb_df, rownames = FALSE)Let us prepare a simple plot of Source and Destination location
#
# Plot source points
#
plot(comb_df[1:2],col= "red", pch = 3,xlab="",ylab="")
#
# Add Destination points
#
points(comb_df[3:4],col="blue")
USA (primarily East Coast) and Europe are seemingly worst hit by plane crashes as appear from the graph.
A better visual can be created using rworldmap package and ggplot2 package
#
(worldMap <- getMap() ) # rworldmap function to get Long/Lat Polygon data for all countries ## class : SpatialPolygonsDataFrame
## features : 243
## extent : -180, 180, -90, 83.64513 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
## variables : 49
## names : ScaleRank, LabelRank, FeatureCla, SOVEREIGNT, SOV_A3, ADM0_DIF, LEVEL, TYPE, ADMIN, ADM0_A3, GEOU_DIF, GEOUNIT, GU_A3, SU_DIF, SUBUNIT, ...
## min values : 1, 1, Adm-0 country, Afghanistan, AFG, 0, 2, Country, Afghanistan, ABW, 0, Afghanistan, ABW, 0, Afghanistan, ...
## max values : 6, 6, Admin-0 countries, Zimbabwe, ZWE, 1, 2, Sovereign country, Zimbabwe, ZWE, 1, Zimbabwe, ZWE, 0, Zimbabwe, ...#
world.points <- fortify(worldMap) # Convert data into dataframe using fortify from ggplot2## Regions defined for each Polygons#
#
head(world.points)## long lat order hole piece id group
## 1 61.21082 35.65007 1 FALSE 1 Afghanistan Afghanistan.1
## 2 62.23065 35.27066 2 FALSE 1 Afghanistan Afghanistan.1
## 3 62.98466 35.40404 3 FALSE 1 Afghanistan Afghanistan.1
## 4 63.19354 35.85717 4 FALSE 1 Afghanistan Afghanistan.1
## 5 63.98290 36.00796 5 FALSE 1 Afghanistan Afghanistan.1
## 6 64.54648 36.31207 6 FALSE 1 Afghanistan Afghanistan.1#
world.points$region <- world.points$id
#
world.df <- world.points[,c("long","lat","group", "region")]
#
worldmap <- ggplot() +
geom_polygon(data = world.df, aes(x = long, y = lat, group = group)) +
geom_point(aes(x=comb_df[,1], y=comb_df[,2]),color="yellow", size=1) + # Plot Source Locations
geom_point(aes(x=comb_df[,3], y=comb_df[,4]),color="green", size=1) + # Plot Dest Location
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45) +
theme_bw() +
theme(axis.line = element_line(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
#
worldmap
By adding the coord_map line to the ggplot call we can transform the map projection to a globe.
worldmap +coord_map("ortho", orientation=c(40, -40, -10))+theme(panel.background = element_rect(fill = 'gray', colour = 'white'), panel.grid.major = element_line(color = "white"),
panel.grid.minor = element_line(color = "white"))
By changing the parameters of the coord_map function we can change the map projection. Lots of information about map projections can be found here.
The azimuthal equidistant projection is of particular interest to this task as it can include the whole world and correctly shows distances and directions. Great circles (shortest paths between two points on a globe) will apear as straight lines on this projection.
Here I plot the sample points on an azimuthal equidistant projection with the North Pole at the centre of the map. (To avoid a polygon rendering difficulty with this particular projection I’ve removed the Antartica polygon)
worldmap <- ggplot() +
geom_polygon(data = world.df[world.df$region != 'Antarctica',], aes(x = long, y = lat, group = group)) +
geom_point(aes(x=comb_df[,1], y=comb_df[,2]),color="yellow", size=1) + # Plot Source Locations
geom_point(aes(x=comb_df[,3], y=comb_df[,4]),color="green", size=1) + # Plot Dest Location
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45) +
coord_map("azequidist",orientation=c(90, 0, 0))+
theme(panel.background = element_rect(fill = 'azure3'))
#
worldmap
Reverse Geocoding
Let us reverse geocode the Lat/Long information of these 50 most frequent source locations to get the name of the cities using revgeocode function from ggmap package that calls Google API for reverse-geocoding
# Let us also plot top 50 such Source locations from where the flights crahsed frequently
#
# Approximate locations as factors rounding off Long / Lat to 2 decimal places
#
source_da <- factor(sprintf("%.2f:%.2f",comb_df[,2], comb_df[,1]))
#
str(source_da)## Factor w/ 490 levels "-0.10:34.77",..: 343 328 413 350 382 326 25 427 54 213 ...#
head(source_da)## [1] 43.99:-0.24 43.13:5.93 51.07:1.02 44.53:-78.43 48.86:2.34
## [6] 42.94:-75.61
## 490 Levels: -0.10:34.77 -0.23:-78.52 -0.76:122.98 ... 9.74:-84.08#
# A table of Source frequencies
#
freq <- sort(table(source_da), decreasing=TRUE)
#
# The most frequent source airports in these data, possibly airports like NYC, London, Paris ?
#
frequent_destinations <- names(freq)[1:50]
#
# Subset the flight data by source frequency
#
idx <- source_da %in% frequent_destinations
#
# Get the Long Lat data for source airports having most crash frequency
#
LongLat <- unique(comb_df[idx,1:2])
#
frequent_flights <-comb_df[idx,]
#
str(LongLat)## 'data.frame': 50 obs. of 2 variables:
## $ SourceLong: num 2.34 -75.61 13.38 -110.98 -1.26 ...
## $ SourceLat : num 48.9 42.9 52.5 29.1 54.5 ...#
head(LongLat)## SourceLong SourceLat
## 5 2.341200 48.85693
## 6 -75.610702 42.93709
## 8 13.376980 52.51607
## 10 -110.984001 29.12122
## 12 -1.264001 54.54761
## 16 -118.243679 34.05223#
# Call Reverse Geocoding function
#
result <- suppressMessages(suppressWarnings(do.call(rbind,
lapply(1:nrow(LongLat),
function(i)revgeocode(as.numeric(LongLat[i,1:2])))) ))
#
# Let us see the output
#
datatable(subset(result,result != 'NA'))#Creating 3-D Interactive Globe
Now we will have some fun with a 3d interactive globe of Flight Route locations using the threejs library to render a 3d globe.
We will Plot frequent source airports (top 50 calculated above) as bars, and their flights routes as arcs. Adjust arc width and color by frequency.
#
#
(earth <- system.file("images/world.jpg", package="threejs")) # Get the image of globe## [1] "C:/Users/ankitagarwal5/Documents/R/win-library/3.2/threejs/images/world.jpg"#
# 3-D visual Dataframe
#
test_df <- data.frame(origin_lat = frequent_flights[,2], origin_long = frequent_flights[,1], dest_lat = frequent_flights[,4], dest_long = frequent_flights[,3])
#
globejs(img=earth, lat=LongLat[,2], long=LongLat[,1], arcs=test_df,
arcsHeight=0.3, arcsLwd=2, arcsColor="#ffff00", arcsOpacity=0.15,
atmosphere=TRUE, height=800, width = 800)You can click on the Globe and drag it to rotate it. Try it out.
We can plot paths of all the Flight Routes (1000 that we have) as well. Let us see how it looks
#
#
(earth <- system.file("images/world.jpg", package="threejs")) # Get the image of globe## [1] "C:/Users/ankitagarwal5/Documents/R/win-library/3.2/threejs/images/world.jpg"#
# 3-D visual Dataframe
#
test_df <- data.frame(origin_lat = comb_df[,2], origin_long = comb_df[,1], dest_lat = comb_df[,4], dest_long = comb_df[,3])
#
globejs(img=earth, lat=LongLat[,2], long=LongLat[,1], arcs=test_df,
arcsHeight=0.3, arcsLwd=2, arcsColor="red", arcsOpacity=0.15,
atmosphere=TRUE,bg="white", height = 800 , width = 800)You can click on the Globe and drag it to rotate it. Try it out.
Hope you have enjoyed this 3-D visualization of Globe.