Making a map of flights

gistfile1.r

# Load Airports
# This comes from: http://openflights.org/data.html
airports <- read.csv('airports.dat', header=F);
colnames(airports) <- c('id', 'name', 'city', 'county', 'faa', 'icao', 'lat', 'lng', 'alt', 'timezone', 'dst')
airports$faa <- as.character(airports$faa)

# Load trips
# a file with: "source-airport,destination-airport,date"
# e.g., "MSP,SFO,2014-01-01"
trips <- read.csv('trips.csv', header=F);
colnames(trips) <- c('source', 'dest', 'date')
trips$source <- as.character(trips$source)
trips$dest <- as.character(trips$dest)
trips$cnt <- 1

# Remove directionality:
trips.2 <- rbind(cbind(trips$source, trips$dest, trips$date, trips$cnt), cbind(trips$dest, trips$source, trips$date, trips$cnt))
trips.2 <- as.data.frame(trips.2)
colnames(trips.2) <- c('source', 'dest', 'date', 'cnt')
trips.2$source <- as.character(trips.2$source)
trips.2$dest <- as.character(trips.2$dest)
trips.2$cnt <- 1
trips.2 <- trips.2[which(trips.2$source > trips.2$dest),]

trips <- aggregate(trips$cnt, by=list(trips$source,trips$dest), FUN=sum)
colnames(trips) <- c('source', 'dest', 'count')

trips.dest <- cbind(trips$source, trips$dest)
trips.dest <- table(trips.dest)
trips.dest <- data.frame(cnt=trips.dest)
trips.dest$cnt.trips.dest <- as.character(trips.dest$cnt.trips.dest)

# Merge (for source)
trips <- merge(trips, airports, by.x='source', by.y='faa', all.x=T);
trips <- trips[,c('source','dest','count','lat','lng')]
colnames(trips)[4:5] <- c('source.lat', 'source.lng')

# Merge (for dest)
trips <- merge(trips, airports, by.x='dest', by.y='faa', all.x=T);
trips <- trips[,c('source', 'dest','count', 'source.lat', 'source.lng', 'lat', 'lng')]
colnames(trips)[6:7] <- c('dest.lat', 'dest.lng')

# Projections
trips.source.project <- mapproject(list(y=trips$source.lat,x=trips$source.lng))
trips.dest.project <- mapproject(list(y=trips$dest.lat,x=trips$dest.lng))
trips$source.lat.p <- trips.source.project$y
trips$source.lng.p <- trips.source.project$x
trips$dest.lat.p <- trips.dest.project$y
trips$dest.lng.p <- trips.dest.project$x
trips$grp <- 1:nrow(trips)

# Airport counts
trips <- merge(trips, trips.dest, by.x='source', by.y='cnt.trips.dest', all.x=T)
colnames(trips)[ncol(trips)] <- 'source.cnt'
trips <- merge(trips, trips.dest, by.x='dest', by.y='cnt.trips.dest', all.x=T)
colnames(trips)[ncol(trips)] <- 'dest.cnt'

# Minor formatting
t <- rbind(cbind(trips$source.lng.p, trips$source.lat.p, trips$grp, trips$count, trips$source.cnt),
cbind(trips$dest.lng.p, trips$dest.lat.p, trips$grp, trips$count, trips$dest.cnt))
t <- as.data.frame(t)
colnames(t) <- c('lng.p', 'lat.p', 'grp', 'count', 'freq.cnt')

# Build the actual map
states <- map_data("usa",project="albers",par=c(39,45))
p <- ggplot(t, aes(lng.p, lat.p, group=grp))
p <- p + geom_polygon(data=states, aes(x=long, y=lat, group=group), colour="#333333", fill="#111111", size=0.2)
p <- p + theme(	panel.grid.major = element_line(colour=NA),
				panel.grid.minor = element_line(colour=NA),
				panel.background = element_rect(fill="black", colour=NA),
				axis.ticks = element_line(colour=NA),
				axis.title.x = element_blank(),
				axis.title.y = element_blank(),
				axis.text.x = element_blank(),
				axis.text.y = element_blank());
p <- p + geom_point(aes(size=sqrt(freq.cnt)), colour='#008c9e', alpha=0.3) + geom_line(aes(colour=count), size=0.4)
plot(p)