Created
January 26, 2014 17:49
-
-
Save epijim/8636599 to your computer and use it in GitHub Desktop.
greatcircles
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ############################################################################# | |
| ####### Prelim stages ########## | |
| ############################################################################# | |
| library(maps) | |
| library(geosphere) | |
| library(RCurl) | |
| library(RJSONIO) | |
| library(ggmap) | |
| library(qdap) | |
| library(mapdata) | |
| setwd("~/Dropbox/Personal/fun with code/migrations") | |
| map("world", col="#f2f2f2", | |
| fill=TRUE, | |
| bg="white", | |
| lwd=0.05) | |
| # Load up my facebook data (downloaded using RFacebook package) | |
| my_friends_info <- read.csv( | |
| "~/Dropbox/Personal/fun with code/migrations/my_friends_info.csv") | |
| places <- read.csv( | |
| "~/Dropbox/Personal/fun with code/migrations/places.csv") | |
| # Make dataframe of people with no missing hometown and current | |
| migrations <- na.omit(my_friends_info) | |
| ############################################################################# | |
| ####### Geocode the locations ########## | |
| ############################################################################# | |
| ## geocode locations using google. | |
| ## Each location goes through as one http request, so it's pretty slow | |
| ## Hence I commented it out and am loading it up for the analysis | |
| # places$location <- as.character(places$location) | |
| # locations <- cbind(places, t(sapply(places$location,geocode, USE.NAMES=F))) | |
| # locations$lon <- as.numeric(locations$lon) | |
| # locations$lat <- as.numeric(locations$lat) | |
| # locations$number <- NULL | |
| # write.table(locations, "locations.csv", col.names=TRUE, sep=",") | |
| locations <- read.csv("~/Dropbox/Personal/fun with code/migrations/locations.csv") | |
| ############################################################################# | |
| ####### Make the map ########## | |
| ############################################################################# | |
| ######################################## | |
| ####### 1st go | |
| ######################################## | |
| map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.05) | |
| for (j in 1:length(migrations$location)) { | |
| place1 <- locations[as.character(locations$location) == as.character(migrations[j,]$hometown),] | |
| place2 <- locations[as.character(locations$location) == as.character(migrations[j,]$location),] | |
| inter <- gcIntermediate(c(place1[1,]$lon, place1[1,]$lat), c(place2[1,]$lon, place2[1,]$lat), | |
| n=100, | |
| breakAtDateLine=F, | |
| addStartEnd=TRUE) | |
| lines(inter, lwd=0.4) | |
| } | |
| ## Ahhh.... Seems it draws a straight line accross when someone crossed th international | |
| ## date line. Need to sort though out. | |
| ## So for the next try I will reformat the data so I can split the plots into left and right if they cross | |
| ## the dateline. | |
| ######################################## | |
| ####### Make the map - take 2 | |
| ######################################## | |
| ## make a new dataframe to play with | |
| mapdata <- migrations | |
| ## do a lookup to pull the location data, to keep things simple below (suggested on the flowing data comments) | |
| mapdata$long1 <- lookup(mapdata$location, locations[, c("location","lon")]) | |
| mapdata$lat1 <- lookup(mapdata$location, locations[, c("location","lat")]) | |
| mapdata$long2 <- lookup(mapdata$hometown, locations[, c("location","lon")]) | |
| mapdata$lat2 <- lookup(mapdata$hometown, locations[, c("location","lat")]) | |
| ######################################## | |
| ####### Make the map | |
| ######################################## | |
| pdf("MAP.pdf") | |
| map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.000001) | |
| lines <- nrow(mapdata) # Get the number of lines | |
| this_width=.5 # Choose line width | |
| this_colour = "black" | |
| ## Start the loop | |
| for (this_one in 1:lines) { | |
| long1 <- mapdata$long1[this_one] | |
| lat1 <- mapdata$lat1[this_one] | |
| long2 <- mapdata$long2[this_one] | |
| lat2 <- mapdata$lat2[this_one] | |
| inter <- gcIntermediate(c(long1, lat1), c(long2, lat2), | |
| n=100, | |
| breakAtDateLine=T, | |
| addStartEnd=TRUE) | |
| ## Step one, draw the people that dont cross the dateline | |
| if(is.numeric(inter)==T) { | |
| lines(inter, col=this_colour, lwd=this_width) | |
| } | |
| ## Draw left and right side if the person crosses | |
| if(is.numeric(inter)==F) { | |
| #Draw the GC on one side of the dateline in red | |
| lines(inter[[1]], col=this_colour, lwd=this_width) | |
| #And the other side of the dateline in green | |
| lines(inter[[2]], col=this_colour, lwd=this_width) | |
| } | |
| } | |
| dev.off() | |
| ######################################## | |
| ####### Map with points | |
| ######################################## | |
| png("MAPpoints.png", width = 800) | |
| map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.000001) | |
| lines <- nrow(mapdata) # Get the number of lines | |
| this_width=.5 # Choose line width | |
| this_colour = "#00000064" | |
| ## Start the loop | |
| for (this_one in 1:lines) { | |
| long1 <- mapdata$long1[this_one] | |
| lat1 <- mapdata$lat1[this_one] | |
| long2 <- mapdata$long2[this_one] | |
| lat2 <- mapdata$lat2[this_one] | |
| inter <- gcIntermediate(c(long1, lat1), c(long2, lat2), | |
| n=100, | |
| breakAtDateLine=T, | |
| addStartEnd=TRUE) | |
| ## Step one, draw the people that dont cross the dateline | |
| if(is.numeric(inter)==T) { | |
| lines(inter, col=this_colour, lwd=this_width) | |
| } | |
| ## Draw left and right side if the person crosses | |
| if(is.numeric(inter)==F) { | |
| #Draw the GC on one side of the dateline in red | |
| lines(inter[[1]], col=this_colour, lwd=this_width) | |
| #And the other side of the dateline in green | |
| lines(inter[[2]], col=this_colour, lwd=this_width) | |
| } | |
| points(jitter(mapdata$long1,20), jitter(mapdata$lat1,20), | |
| col="#FF00000A", bg="#000000", pch=21, cex=0.03) | |
| points(jitter(mapdata$long2,20), jitter(mapdata$lat2,20), | |
| col="#0000FF0A", bg="#000000", pch=21, cex=0.03) | |
| } | |
| dev.off() | |
| ############################################################################# | |
| ####### What % of my friends have hometown and current data? ########## | |
| ############################################################################# | |
| 100*(nrow(mapdata) / nrow(my_friends_info)) | |
| # 49.6% have migration data | |
| # how many still in hometown? | |
| 54 / nrow(my_friends_info) | |
| # 11.0 % are still in their hometown | |
| ############################################################################# | |
| ####### Distance between two points ########## | |
| ############################################################################# | |
| # Convert degrees to radians | |
| deg2rad <- function(deg) return(deg*pi/180) | |
| mapdata$longrad1 <- deg2rad(mapdata$long1) | |
| mapdata$longrad2 <- deg2rad(mapdata$long2) | |
| mapdata$latrad1 <- deg2rad(mapdata$lat1) | |
| mapdata$latrad2 <- deg2rad(mapdata$lat2) | |
| # Calculates the geodesic distance between two points specified by radian latitude/longitude using the | |
| # Spherical Law of Cosines (slc) | |
| gcd.slc <- function(long1, lat1, long2, lat2) { | |
| R <- 6371 # Earth mean radius [km] | |
| d <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R | |
| return(d) # Distance in km | |
| } | |
| mapdata$distance <- gcd.slc(mapdata$longrad1,mapdata$latrad1,mapdata$longrad2,mapdata$latrad2) | |
| ## Histogram | |
| png("histogram.png",width=800) | |
| hist(mapdata$distance, freq=T, xlab="KM between current location and hometown", | |
| main="Number of friends by distance migrated", | |
| ylab="Number", | |
| col="red", | |
| labels=T, | |
| breaks=10, | |
| xlim=c(0,20000), ylim=c(0, 200)) | |
| dev.off() | |
| kiwi <- read.csv("~/Dropbox/Personal/fun with code/migrations/mapdata.csv") | |
| library(beeswarm) | |
| png("swarm.png", width=800) | |
| beeswarm(distance ~ kiwi, data = kiwi, | |
| method = 'hex', | |
| #pch = 1, | |
| col = "Red", | |
| cex=0.7, | |
| xlab = 'Distance current location is from home (KM)', ylab = 'Hometown in NZ?', | |
| horizontal=T, | |
| #corral="gutter", | |
| labels = c('No', 'Yes')) | |
| boxplot(distance ~ kiwi, data = kiwi, | |
| add = T, names = c("",""), | |
| horizontal=T, | |
| outline=F, | |
| col="#0000ff22") | |
| dev.off() | |
| ############################################################################# | |
| ####### New Zealand ########## | |
| ############################################################################# | |
| pdf("MAPnz.pdf") | |
| map('worldHires', | |
| c('New Zealand'), | |
| xlim=c(161,180), ylim=c(-50,-31)) | |
| lines <- nrow(mapdata) # Get the number of lines | |
| this_width=1 # Choose line width | |
| this_colour = col=rgb(0, 0, 0, 0.1) | |
| ## Start the loop | |
| for (this_one in 1:lines) { | |
| long1 <- mapdata$long1[this_one] | |
| lat1 <- mapdata$lat1[this_one] | |
| long2 <- mapdata$long2[this_one] | |
| lat2 <- mapdata$lat2[this_one] | |
| inter <- gcIntermediate(c(long1, lat1), c(long2, lat2), | |
| n=100, | |
| breakAtDateLine=T, | |
| addStartEnd=TRUE) | |
| ## Step one, draw the people that dont cross the dateline | |
| if(is.numeric(inter)==T) { | |
| lines(inter, col=this_colour, lwd=this_width) | |
| } | |
| ## Draw left and right side if the person crosses | |
| if(is.numeric(inter)==F) { | |
| #Draw the GC on one side of the dateline in red | |
| lines(inter[[1]], col=this_colour, lwd=this_width) | |
| #And the other side of the dateline in green | |
| lines(inter[[2]], col=this_colour, lwd=this_width) | |
| } | |
| points(mapdata$long1, mapdata$lat1, | |
| col=rgb(1, 0, 0, 0.5), bg="#000000", pch=2, cex=0.5) | |
| points(mapdata$long2, mapdata$lat2, | |
| col=rgb(0, 0, 1, 0.1), bg="#000000", pch=6, cex=0.5) | |
| } | |
| dev.off() | |
| ############################################################################# | |
| ####### Dot maps of home and now ########## | |
| ############################################################################# | |
| pdf("locationnow.pdf") | |
| ## Location now | |
| map("world", col="#191919", fill=TRUE, bg="#000000", lwd=0.000001) | |
| points(jitter(mapdata$long1,20), jitter(mapdata$lat1,20), | |
| col="White", bg="#000000", pch=21, cex=0.03) | |
| dev.off() | |
| ## Hometown | |
| pdf("hometown.pdf") | |
| map("world", col="#191919", fill=TRUE, bg="#000000", lwd=0.000001) | |
| points(jitter(mapdata$long2,20), jitter(mapdata$lat2,20), | |
| col="White", bg="#000000", pch=21, cex=0.03) | |
| dev.off() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment