自学R的新手一枚。
看到某大神的帖子后模仿改进了下，希望大家喜欢。
Mapping the World’s Biggest Airlines
参考资料：http://spatialanalysis.co.uk/2012/06/mapping-worlds-biggest-airlines/
看了这位大神的地图后，发现自己真的小瞧了R的强大。
正好最近有相关业务需求就照葫芦画瓢做了一个。效果图如下：
不会贴图，可以去我博客看：
http://blog.sina.com.cn/s/blog_3f16beff0102v38u.html
R语言飞机航线夜景图研究
比那位大神的图多了光晕效果，图片更绚丽。
使用R版本3.1.0
大致思路如下：
1、数据处理：将shp数据转化为RData数据
2、地图旋转，中国是世界中心，自己写的算法效率一般般啦。(坐标变换，图型闭合处理)
3、绘制图像
相关资源地址：
航线、飞机场坐标：
http://openflights.org/data.html
板块地图、都市地图：
http://www.naturalearthdata.com/downloads/
R语言代码如下：
<br />
#程序中差不多用了这么些包，也许有没用到的，懒得查了<br />
#install.packages("ggmap")<br />
library(ggmap)<br />
library(ggplot2)<br />
library(maps)<br />
#install.packages("rgeos")<br />
library(rgeos)<br />
library(maptools)<br />
library(shapefiles)<br />
library(geosphere)<br />
library(plyr)<br />
library(sp)</p>
<p>#读取都市地图文件<br />
#读取版图地图文件<br />
urbanareasin<-readShapePoly('ne_10m_urban_areas.shp')<br />
worldmapsin<-readShapePoly('ne_10m_admin_0_countries.shp')<br />
#以下为格式转化<br />
worldmap<-fortify(worldmapsin)<br />
urbanareas<-fortify(urbanareasin)</p>
<p>gpclibPermit()<br />
#开始抽取机场数据<br />
worldport<-airports[airports$V5!='',c('V3','V5','V7','V8','V9')]<br />
names(worldport)<- c('city','code','lan','lon','att')<br />
worldport$lan<-as.numeric(as.character(worldport$lan))<br />
worldport$lon<-as.numeric(as.character(worldport$lon))<br />
#找出所有航线有标识的机场（这里的data3.redu2s是我个人的航线数据，你们可以用上文提到的地址的航线数据代替）<br />
lineinworld<-(data3.redu2s$AIRPORT_FROM_CODE%in%worldport$code)&(data3.redu2s$AIRPORT_TO_CODE%in%worldport$code)</p>
<p>#有243条航线无标识<br />
table(lineinworld)<br />
#colnames(data3.upro1)<br />
worldline<-data3.redu2s[lineinworld,c("AIRPORT_FROM_CODE","AIRPORT_TO_CODE")]<br />
flights.ag<-ddply(worldline, c("AIRPORT_FROM_CODE","AIRPORT_TO_CODE"), function(x) count(x$AIRPORT_TO_CODE))<br />
#计算三字码映射到机场<br />
flights.ll <- merge(flights.ag, worldport, by.x="AIRPORT_TO_CODE", by.y="code", all.x=T)<br />
flights.ll<-flights.ll[order(flights.ll$AIRPORT_FROM_CODE,flights.ll$AIRPORT_TO_CODE),]<br />
flights.lf <- merge(flights.ag, worldport, by.x="AIRPORT_FROM_CODE", by.y="code", all.x=T)<br />
flights.lf<-flights.lf[order(flights.lf$AIRPORT_FROM_CODE,flights.lf$AIRPORT_TO_CODE),]</p>
<p>#beijing.ll <- c(worldport$lon[worldport$code=="PEK"],worldport$lan[worldport$code=="PEK"])<br />
rts <- gcIntermediate(flights.lf[,c('lon', 'lan')], flights.ll[,c('lon', 'lan')], 100, breakAtDateLine=FALSE, addStartEnd=TRUE, sp=TRUE)<br />
#rts.ff <- fortify.SpatialLinesDataFrame(rts)flights.lf[,c('lon', 'lan')]<br />
rts <- as(rts, "SpatialLinesDataFrame")<br />
#航线坐标数据<br />
rts.ff <- fortify(rts)<br />
flights.ll$id <-as.character(c(1:nrow(flights.ll))) # that rts.ff$id is a char<br />
#航线信息与航线坐标信息关联<br />
table(gcircles$freq)<br />
gcircles <- merge(rts.ff, flights.ll, all.x=T, by="id") # join attributes, we keep them all, just in case<br />
#设置地图中心<br />
center <- 115 # positive values only - US centered view is 260</p>
<p># shift coordinates to recenter great circles<br />
#航线坐标计算中心距离<br />
gcircles$long.recenter <- ifelse(gcircles$long < center - 180 , gcircles$long + 360, gcircles$long)</p>
<p># shift coordinates to recenter worldmap<br />
#worldmap <- map_data ("world")<br />
#地图坐标偏移<br />
worldmap$long.recenter <- ifelse(worldmap$long < center - 180 , worldmap$long + 360, worldmap$long)<br />
urbanareas$long.recenter <- ifelse(urbanareas$long < center - 180 , urbanareas$long + 360, urbanareas$long)<br />
RegroupElements <- function(df, longcol, idcol){<br />
g <- rep(1, length(df[,longcol]))<br />
if (diff(range(df[,longcol])) > 300) { # check if longitude within group differs more than 300 deg, ie if element was split<br />
d <- df[,longcol] > mean(range(df[,longcol])) # we use the mean to help us separate the extreme values<br />
g[!d] <- 1 # some marker for parts that stay in place (we cheat here a little, as we do not take into account concave polygons)<br />
g[d] <- 2 # parts that are moved<br />
}<br />
g <- paste(df[, idcol], g, sep=".") # attach to id to create unique group variable for the dataset<br />
df$group.regroup <- g<br />
df<br />
}</p>
<p>#重写上面的语句<br />
gcircles.rg <- ddply(gcircles, .(id), RegroupElements, "long.recenter", "id")<br />
#压缩地图<br />
#以下方法是原始算法，但是效率极低，故在下面重写<br />
#worldmap.rg <- ddply(worldmap, .(group), RegroupElements, "long.recenter", "group")<br />
#urbanareas.rg <- ddply(urbanareassim, .(group), RegroupElements, "long.recenter", "group")</p>
<p>#开始写原始算法替换函数<br />
worldmap.mean<-aggregate(x=worldmap[,c("long.recenter")],by=list(worldmap$group),FUN=mean)<br />
#worldmap.diff<-aggregate(x=worldmap[,c("long.recenter")],by=list(worldmap$group),FUN=function(x){diff(range(x))})<br />
worldmap.min<-aggregate(x=worldmap[,c("long.recenter")],by=list(worldmap$group),FUN=min)<br />
worldmap.max<-aggregate(x=worldmap[,c("long.recenter")],by=list(worldmap$group),FUN=max)<br />
#worldmap.md<-cbind(worldmap.mean,worldmap.diff[,2])<br />
worldmap.md<-cbind(worldmap.mean,worldmap.min[,2],worldmap.max[,2])<br />
#colnames(worldmap.md)<-c("group","mean","diff")<br />
colnames(worldmap.md)<-c("group","mean","min","max")<br />
worldmapt<-join(x=worldmap,y=worldmap.md,by = c("group"))<br />
worldmapt$group.regroup<-1<br />
worldmapt[(worldmapt$max>180)&(worldmapt$min<180)&(worldmapt$long.recenter>180),c("group.regroup")]<-2<br />
worldmapt$group.regroup<-paste(worldmapt$group,worldmapt$group.regroup,sep=".")<br />
worldmap.rg<-worldmapt</p>
<p>urbanareas.mean<-aggregate(x=urbanareas[,c("long.recenter")],by=list(urbanareas$group),FUN=mean)<br />
#urbanareas.diff<-aggregate(x=urbanareas[,c("long.recenter")],by=list(urbanareas$group),FUN=function(x){diff(range(x))})<br />
urbanareas.min<-aggregate(x=urbanareas[,c("long.recenter")],by=list(urbanareas$group),FUN=min)<br />
urbanareas.max<-aggregate(x=urbanareas[,c("long.recenter")],by=list(urbanareas$group),FUN=max)<br />
#urbanareas.md<-cbind(urbanareas.mean,urbanareas.diff[,2])<br />
urbanareas.md<-cbind(urbanareas.mean,urbanareas.min[,2],urbanareas.max[,2])<br />
#colnames(urbanareas.md)<-c("group","mean","diff")<br />
colnames(urbanareas.md)<-c("group","mean","min","max")<br />
urbanareast<-join(x=urbanareas,y=urbanareas.md,by = c("group"))<br />
urbanareast$group.regroup<-1<br />
urbanareast[(urbanareast$max>180)&(urbanareast$min<180)&(urbanareast$long.recenter>180),c("group.regroup")]<-2<br />
urbanareast$group.regroup<-paste(urbanareast$group,urbanareast$group.regroup,sep=".")<br />
urbanareas.rg<-urbanareast</p>
<p>#worldmap.cp <- ddply(worldmap.rg, .(group.regroup), ClosePolygons, "long.recenter", "order") # use the new grouping var<br />
#闭合曲线<br />
worldmap.rg<-worldmap.rg[order(worldmap.rg$group.regroup,worldmap.rg$order),]<br />
worldmap.begin<-worldmap.rg[!duplicated(worldmap.rg$group.regroup),]<br />
worldmap.end<-worldmap.rg[c(!duplicated(worldmap.rg$group.regroup)[-1],TRUE),]<br />
worldmap.flag<-(worldmap.begin$long.recenter==worldmap.end$long.recenter)&(worldmap.begin$lat==worldmap.end$lat)<br />
table(worldmap.flag)<br />
worldmap.plus<-worldmap.begin[!worldmap.flag,]<br />
worldmap.end[!worldmap.flag,]<br />
worldmap.plus$order<-worldmap.end$order[!worldmap.flag]+1<br />
worldmap.cp<-rbind(worldmap.rg,worldmap.plus)<br />
worldmap.cp<-worldmap.cp[order(worldmap.cp$group.regroup,worldmap.cp$order),]</p>
<p>urbanareas.rg<-urbanareas.rg[order(urbanareas.rg$group.regroup,urbanareas.rg$order),]<br />
urbanareas.begin<-urbanareas.rg[!duplicated(urbanareas.rg$group.regroup),]<br />
urbanareas.end<-urbanareas.rg[c(!duplicated(urbanareas.rg$group.regroup)[-1],TRUE),]<br />
urbanareas.flag<-(urbanareas.begin$long.recenter==urbanareas.end$long.recenter)&(urbanareas.begin$lat==urbanareas.end$lat)<br />
table(urbanareas.flag)<br />
urbanareas.plus<-urbanareas.begin[!urbanareas.flag,]<br />
urbanareas.end[!urbanareas.flag,]<br />
urbanareas.plus$order<-urbanareas.end$order[!urbanareas.flag]+1<br />
urbanareas.cp<-rbind(urbanareas.rg,urbanareas.plus)<br />
urbanareas.cp<-urbanareas.cp[order(urbanareas.cp$group.regroup,urbanareas.cp$order),]</p>
<p>range(worldmap.cp$long)<br />
#urbanareast<-join(x=urbanareas,y=urbanareas.md,by = c("group"))<br />
# close polys<br />
#worldmap.cp <- ddply(worldmap.rg, .(group.regroup), ClosePolygons, "long.recenter", "order") # use the new grouping var<br />
wrld<- c(geom_polygon(aes(long.recenter,lat,group=group.regroup),<br />
size = 0.1,<br />
colour= "#090D2A",<br />
fill="#090D2A",<br />
alpha=1,<br />
data=worldmap.cp))<br />
urb <- c(geom_polygon(aes(long.recenter,lat,group=group.regroup),<br />
size = 0.3,<br />
color = "#FDF5E6",<br />
fill = "#FDF5E6",<br />
alpha = 1,<br />
data = urbanareas.cp))</p>
<p>#放大系数<br />
bigmap<-1<br />
airline<- geom_line(aes(long.recenter,lat,group=group.regroup,alpha=max(freq)^0.6*freq^0.4,color=0.9*max(freq)^0.6*freq^0.4),<br />
size=0.2*bigmap,<br />
data= gcircles.rg)<br />
airlinep<- geom_line(aes(long.recenter,lat,group=group.regroup,alpha=0.04*max(freq)^0.6*freq^0.4),<br />
color="#FFFFFF",<br />
size=2*bigmap,<br />
#alpha=0.08,<br />
data= gcircles.rg)<br />
#table(gcircles.rg$freq)<br />
airlineppp<- geom_line(aes(long.recenter,lat,group=group.regroup,alpha=0.02*max(freq)^0.6*freq^0.4),<br />
color="#ECFFFF",<br />
size=4*bigmap,<br />
#alpha=0.015,<br />
data= gcircles.rg)<br />
airlineppp<- geom_line(aes(long.recenter,lat,group=group.regroup,alpha=0.010*max(freq)^0.6*freq^0.4),<br />
color="#ECFFFF",<br />
size=8*bigmap,<br />
#alpha=0.015,<br />
data= gcircles.rg)<br />
airlinepppp<- geom_line(aes(long.recenter,lat,group=group.regroup,alpha=0.005*max(freq)^0.6*freq^0.4),<br />
color="#BBFFFF",<br />
size=16*bigmap,<br />
#alpha=0.005,<br />
data= gcircles.rg)</p>
<p>gcircles.rg[gcircles.rg$group.regroup==1.1,]</p>
<p># plot画图到文件plot2.png<br />
png(6000*(bigmap^0.8), 2000*(bigmap^0.8), file="plot2.png", bg="white")<br />
ggplot() +<br />
wrld+<br />
urb+<br />
airline+<br />
scale_colour_gradient(low="#D9FFFF", high="#ECFFFF")+<br />
scale_alpha(range = c(0, 1))+<br />
airlineppp+<br />
#scale_alpha_discrete(range = c(0.001, 0.005))+<br />
airlinepp+<br />
#scale_alpha(range = c(0.005, 0.015))+<br />
airlinep+<br />
#scale_alpha(range = c(0.015, 0.08))+<br />
#geom_polygon(aes(long,lat,group=group), size = 0.2, fill="#f9f9f9", colour = "grey65", data=worldmap) +<br />
#geom_line(aes(long.recenter,lat,group=group.regroup, color=freq, alpha=freq), size=0.4, data= gcircles.rg) + # set transparency here # set color gradient here<br />
theme(panel.background = element_rect(fill = '#00001C', color = '#00001C'))+<br />
theme(panel.grid.minor = element_blank())+<br />
theme(panel.grid.major = element_blank())+<br />
theme(axis.ticks = element_blank())+<br />
theme(axis.title.x = element_blank())+<br />
theme(axis.title.y = element_blank())+<br />
theme(axis.text.x = element_blank())+<br />
theme(axis.text.y = element_blank())+<br />
theme(legend.position = "none") +<br />
ylim(-65, 75) +<br />
coord_equal()<br />
dev.off()<br />
#我是寂寞的数据分析师
</p>
