knitr::opts_chunk$set(echo = TRUE, eval = TRUE, tidy = TRUE, highlight = TRUE, warning = FALSE, error = FALSE, message = FALSE, include = TRUE,fig.width = 9.4,fig.height = 4)
1、安装、加载相关包
# install.packages('ggmap')
library(ggmap)
library(ggplot2)
library(maps)
# install.packages('rgeos')
library(rgeos)
# install.packages('maptools')
library(maptools)
# install.packages('shapefiles')
library(shapefiles)
library(geosphere)
library(plyr)
library(sp)
# install.packages('rgdal') install.packages('sf')
library(rgdal)
2.1 机场信息
airports <- read.table("/resources/rstudio/airport/airports.dat", sep = ",",
header = F)
colnames(airports) <- c("Airport_ID", "Name", "City", "Country", "IATA_FAA",
"ICAO", "Latitude", "Longitude", "Altitude", "Timezone", "DST", "Tz_database_time_zone")
str(airports)
## 'data.frame': 8107 obs. of 12 variables:
## $ Airport_ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Name : Factor w/ 7908 levels "\\\\'t Harde",..: 2582 4235 4770 4854 5633 7611 4913 5114 6588 7058 ...
## $ City : Factor w/ 6977 levels "108 Mile Ranch",..: 2197 3641 4118 4197 4933 6719 4258 2171 5812 6206 ...
## $ Country : Factor w/ 240 levels "Afghanistan",..: 169 169 169 169 169 169 84 84 84 84 ...
## $ IATA_FAA : Factor w/ 5879 levels "","%u0","04G",..: 1709 2950 1906 2644 3914 5330 4987 1748 4378 4774 ...
## $ ICAO : Factor w/ 6784 levels "","\\N","%u04",..: 54 62 63 68 69 75 81 86 95 97 ...
## $ Latitude : num -6.08 -5.21 -5.83 -6.57 -9.44 ...
## $ Longitude : num 145 146 144 147 147 ...
## $ Altitude : int 5282 20 5388 239 146 19 112 283 165 251 ...
## $ Timezone : num 10 10 10 10 10 10 -3 -3 -3 -4 ...
## $ DST : Factor w/ 7 levels "A","E","N","O",..: 6 6 6 6 6 6 2 2 2 2 ...
## $ Tz_database_time_zone: Factor w/ 294 levels "\\N","Africa/Abidjan",..: 286 286 286 286 286 286 79 79 79 120 ...
knitr::kable(head(airports))
| 1 |
Goroka |
Goroka |
Papua New Guinea |
GKA |
AYGA |
-6.081689 |
145.3919 |
5282 |
10 |
U |
Pacific/Port_Moresby |
| 2 |
Madang |
Madang |
Papua New Guinea |
MAG |
AYMD |
-5.207083 |
145.7887 |
20 |
10 |
U |
Pacific/Port_Moresby |
| 3 |
Mount Hagen |
Mount Hagen |
Papua New Guinea |
HGU |
AYMH |
-5.826789 |
144.2959 |
5388 |
10 |
U |
Pacific/Port_Moresby |
| 4 |
Nadzab |
Nadzab |
Papua New Guinea |
LAE |
AYNZ |
-6.569828 |
146.7262 |
239 |
10 |
U |
Pacific/Port_Moresby |
| 5 |
Port Moresby Jacksons Intl |
Port Moresby |
Papua New Guinea |
POM |
AYPY |
-9.443383 |
147.2200 |
146 |
10 |
U |
Pacific/Port_Moresby |
| 6 |
Wewak Intl |
Wewak |
Papua New Guinea |
WWK |
AYWK |
-3.583828 |
143.6692 |
19 |
10 |
U |
Pacific/Port_Moresby |
cn_airports <- airports[airports$Country == "China", ] # 中国机场
2.2 航空公司信息
airlines <- read.table("/resources/rstudio/airport/airlines.dat", sep = ",",
header = F)
colnames(airlines) <- c("Airline_ID", "Name", "Alias", "IATA", "ICAO", "Callsign",
"Country", "Active")
str(airlines)
## 'data.frame': 6048 obs. of 8 variables:
## $ Airline_ID: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Name : Factor w/ 5959 levels "1-2-go","12 North",..: 4356 3 4 5 6 7 8 9 11 12 ...
## $ Alias : Factor w/ 140 levels "","\\N","Aero Spike",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ IATA : Factor w/ 1091 levels "","-","--","-.",..: 2 1 45 1 1 1 1 1 1 774 ...
## $ ICAO : Factor w/ 5768 levels "","---","--+",..: 3400 2064 4232 5647 4959 996 5164 5200 4449 3231 ...
## $ Callsign : Factor w/ 5216 levels ""," "," Aerovias de Integracion Regional",..: 1 2160 3329 1 1 1412 1276 1450 4061 3170 ...
## $ Country : Factor w/ 277 levels ""," Boonville Stage Line",..: 1 267 238 265 217 217 217 265 267 267 ...
## $ Active : Factor w/ 3 levels "n","N","Y": 3 2 3 2 2 2 2 2 2 3 ...
2.3 航线信息
routes <- read.table("/resources/rstudio/airport/routes.dat", sep = ",", header = F)
colnames(routes) <- c("Airline", "Airline_ID", "Source_airport", "Source_airport_ID",
"Destination_airport", "Destination_airport_ID", "Codeshare", "Stops", "Equipment")
str(routes)
## 'data.frame': 67663 obs. of 9 variables:
## $ Airline : Factor w/ 568 levels "2B","2G","2I",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Airline_ID : Factor w/ 541 levels "\\N","10","1006",..: 374 374 374 374 374 374 374 374 374 374 ...
## $ Source_airport : Factor w/ 3409 levels "AAE","AAL","AAN",..: 50 152 152 472 472 699 699 699 699 771 ...
## $ Source_airport_ID : Factor w/ 3271 levels "\\N","1","10",..: 1091 1092 1092 1094 1094 1724 1724 1724 1724 2497 ...
## $ Destination_airport : Factor w/ 3418 levels "AAE","AAL","AAN",..: 1557 1557 1901 1557 2154 1557 1979 2785 3001 1403 ...
## $ Destination_airport_ID: Factor w/ 3276 levels "\\N","1","10",..: 1113 1113 1090 1113 1763 1113 2850 1 2501 1081 ...
## $ Codeshare : Factor w/ 2 levels "","Y": 1 1 1 1 1 1 1 1 1 1 ...
## $ Stops : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Equipment : Factor w/ 3946 levels ""," 73W 733 73C",..: 2919 2919 2919 2919 2919 2919 2919 2919 2919 2919 ...
knitr::kable(head(routes))
| 2B |
410 |
AER |
2965 |
KZN |
2990 |
|
0 |
CR2 |
| 2B |
410 |
ASF |
2966 |
KZN |
2990 |
|
0 |
CR2 |
| 2B |
410 |
ASF |
2966 |
MRV |
2962 |
|
0 |
CR2 |
| 2B |
410 |
CEK |
2968 |
KZN |
2990 |
|
0 |
CR2 |
| 2B |
410 |
CEK |
2968 |
OVB |
4078 |
|
0 |
CR2 |
| 2B |
410 |
DME |
4029 |
KZN |
2990 |
|
0 |
CR2 |
3、读取地图数据
3.1 读取都市地图数据文件：
#### http://blog.csdn.net/lichangzai/article/details/52025248
#### file.exists('/resources/rstudio/airport/ne_10m_urban_areas.shp')
#### 检测数据文件是否存在
#### urbanareasin<-readShapePoly('/resources/rstudio/airport/ne_10m_urban_areas.shp')
#### 函数readShapePoly已废弃
#### readOGR()函数读取shp数据，需要把同名的shp、shx、dbf文件放到同一个文件夹里，否则读入数据就会报错
urbanareasin <- readOGR("/resources/rstudio/airport/ne_10m_urban_areas.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/resources/rstudio/airport/ne_10m_urban_areas.shp", layer: "ne_10m_urban_areas"
## with 11878 features
## It has 3 fields
## Integer64 fields read as strings: scalerank
urbanareas <- fortify(urbanareasin) # 数据转换
knitr::kable(head(urbanareas))
| -1.657182 |
37.38898 |
1 |
FALSE |
1 |
0 |
0.1 |
| -1.657182 |
37.38258 |
2 |
FALSE |
1 |
0 |
0.1 |
| -1.656224 |
37.38313 |
3 |
FALSE |
1 |
0 |
0.1 |
| -1.656221 |
37.38313 |
4 |
FALSE |
1 |
0 |
0.1 |
| -1.656124 |
37.38318 |
5 |
FALSE |
1 |
0 |
0.1 |
| -1.656120 |
37.38318 |
6 |
FALSE |
1 |
0 |
0.1 |
3.2 读取板块地图数据文件：
worldmapsin <- readOGR("/resources/rstudio/airport/ne_110m_admin_0_countries.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/resources/rstudio/airport/ne_110m_admin_0_countries.shp", layer: "ne_110m_admin_0_countries"
## with 177 features
## It has 63 fields
worldmap <- fortify(worldmapsin) # 数据转换
knitr::kable(head(worldmap))
| 61.21082 |
35.65007 |
1 |
FALSE |
1 |
0 |
0.1 |
| 62.23065 |
35.27066 |
2 |
FALSE |
1 |
0 |
0.1 |
| 62.98466 |
35.40404 |
3 |
FALSE |
1 |
0 |
0.1 |
| 63.19354 |
35.85717 |
4 |
FALSE |
1 |
0 |
0.1 |
| 63.98290 |
36.00796 |
5 |
FALSE |
1 |
0 |
0.1 |
| 64.54648 |
36.31207 |
6 |
FALSE |
1 |
0 |
0.1 |
3.3 绘制城市分布图
ggplot(worldmap, aes(x = long, y = lat, group = group)) + geom_polygon(fill = "#090D2A",
colour = "#090D2A") + geom_point(aes(x = long, y = lat, group = 1), data = urbanareas,
color = "white", shape = ".", size = 0.01, alpha = 0.3) + theme(panel.grid = element_blank(),
panel.background = element_rect(fill = "black"), axis.text = element_blank(),
axis.ticks = element_blank(), axis.title = element_blank())
3.4 绘制世界板块地图
ggplot(worldmap, aes(x = long, y = lat, group = group)) + geom_polygon(fill = "#090D2A",
colour = "grey") + theme(panel.grid = element_blank(), panel.background = element_rect(fill = "#090D2A"),
axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
3.5 绘制机场分布图
ggplot(worldmap, aes(x = long, y = lat, group = group)) + geom_polygon(fill = "#090D2A",
colour = "#090D2A") + geom_point(aes(x = Longitude, y = Latitude, group = 1),
data = airports, color = "white", shape = ".", size = 0.01, alpha = 0.7) +
theme(panel.grid = element_blank(), panel.background = element_rect(fill = "black"),
axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
4、数据预处理
4.1 处理机场数据，将IATA_FAA为空的移除，重命名变量,保留变量：城市，国家，代码，维度，经度，海拔高度
worldport <- airports[airports$IATA_FAA != "", c(3, 4, 5, 7, 8, 9)]
names(worldport) <- c("city", "country", "code", "lan", "lon", "alt")
head(worldport)
## city country code lan lon alt
## 1 Goroka Papua New Guinea GKA -6.081689 145.3919 5282
## 2 Madang Papua New Guinea MAG -5.207083 145.7887 20
## 3 Mount Hagen Papua New Guinea HGU -5.826789 144.2959 5388
## 4 Nadzab Papua New Guinea LAE -6.569828 146.7262 239
## 5 Port Moresby Papua New Guinea POM -9.443383 147.2200 146
## 6 Wewak Papua New Guinea WWK -3.583828 143.6692 19
4.2 计算航线的频数
library(dplyr)
routes_count <- routes %>% group_by(Source_airport, Destination_airport) %>%
summarise(cnt = n())
dim(routes_count)
## [1] 37595 3
knitr::kable(head(routes_count))
| AAE |
ALG |
1 |
| AAE |
CDG |
1 |
| AAE |
IST |
1 |
| AAE |
LYS |
1 |
| AAE |
MRS |
2 |
| AAE |
ORN |
1 |
4.3 在航线数据中增加出发机场、到达机场经纬度
routes_all <- routes_count %>% left_join(worldport, by = c(Source_airport = "code")) %>%
left_join(worldport, by = c(Destination_airport = "code")) %>% na.omit()
# routes_china <- routes_all[routes_all$country.x == 'China'|
# routes_all$country.y == 'China',]
4.4 提取北京机场的航线数据
routes_pek <- routes_all[routes_all$Source_airport == "PEK" | routes_all$Destination_airport ==
"PEK", ]
routes_pek$id <- row.names(routes_pek)
4.5 将北京航线进行曲线化
rts <- gcIntermediate(routes_pek[, c("lon.x", "lan.x")], routes_pek[, c("lon.y",
"lan.y")], 100, breakAtDateLine = FALSE, addStartEnd = TRUE, sp = TRUE)
rts <- as(rts, "SpatialLinesDataFrame")
# 航线坐标数据
rts.ff <- fortify(rts)
knitr::kable(head(rts.ff))
| 38.79932 |
8.977889 |
1 |
1 |
1 |
1.1 |
| 39.38165 |
9.446543 |
2 |
1 |
1 |
1.1 |
| 39.96557 |
9.914235 |
3 |
1 |
1 |
1.1 |
| 40.55116 |
10.380916 |
4 |
1 |
1 |
1.1 |
| 41.13849 |
10.846533 |
5 |
1 |
1 |
1.1 |
| 41.72766 |
11.311033 |
6 |
1 |
1 |
1.1 |
rts.ff1 <- left_join(rts.ff, routes_pek[, c("id", "cnt")]) # 增加航线频数
knitr::kable(head(rts.ff1))
| 38.79932 |
8.977889 |
1 |
1 |
1 |
1.1 |
2 |
| 39.38165 |
9.446543 |
2 |
1 |
1 |
1.1 |
2 |
| 39.96557 |
9.914235 |
3 |
1 |
1 |
1.1 |
2 |
| 40.55116 |
10.380916 |
4 |
1 |
1 |
1.1 |
2 |
| 41.13849 |
10.846533 |
5 |
1 |
1 |
1.1 |
2 |
| 41.72766 |
11.311033 |
6 |
1 |
1 |
1.1 |
2 |
4.6 绘制航线图
做出的航线图不够美观
ggplot(worldmap, aes(x = long, y = lat, group = group)) + geom_polygon(fill = "#090D2A",
colour = "#090D2A") + geom_point(aes(x = Longitude, y = Latitude, group = 1),
data = airports, color = "white", shape = ".", size = 0.01, alpha = 0.3) +
geom_line(data = rts.ff, aes(x = long, y = lat, group = group), color = "gold",
alpha = 0.95, size = 0.05) + theme(panel.grid = element_blank(), panel.background = element_rect(fill = "black"),
axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
5、地图旋转
center <- 115
# 航线坐标计算中心距离
gcircles <- rts.ff1
gcircles$long.recenter <- ifelse(gcircles$long < center - 180, gcircles$long +
360, gcircles$long)
gcircles.mean <- aggregate(x = gcircles[, c("long.recenter")], by = list(gcircles$id),
FUN = mean)
gcircles.min <- aggregate(x = gcircles[, c("long.recenter")], by = list(gcircles$id),
FUN = min)
gcircles.max <- aggregate(x = gcircles[, c("long.recenter")], by = list(gcircles$id),
FUN = max)
gcircles.md <- cbind(gcircles.mean, gcircles.min[, 2], gcircles.max[, 2])
colnames(gcircles.md) <- c("id", "mean", "min", "max")
gcirclest <- join(x = gcircles, y = gcircles.md, by = c("id"))
gcirclest$group.regroup <- 1
gcirclest[(gcirclest$max > 180) & (gcirclest$min < 180) & (gcirclest$long.recenter >
180), c("group.regroup")] <- 2
gcirclest$group.regroup <- paste(gcirclest$id, gcirclest$id.regroup, sep = ".")
gcircles.rg <- gcirclest
# 地图坐标偏移
worldmap$long.recenter <- ifelse(worldmap$long < center - 180, worldmap$long +
360, worldmap$long)
urbanareas$long.recenter <- ifelse(urbanareas$long < center - 180, urbanareas$long +
360, urbanareas$long)
## 世界地图重分组
worldmap.mean <- aggregate(x = worldmap[, c("long.recenter")], by = list(worldmap$group),
FUN = mean)
worldmap.min <- aggregate(x = worldmap[, c("long.recenter")], by = list(worldmap$group),
FUN = min)
worldmap.max <- aggregate(x = worldmap[, c("long.recenter")], by = list(worldmap$group),
FUN = max)
worldmap.md <- cbind(worldmap.mean, worldmap.min[, 2], worldmap.max[, 2])
colnames(worldmap.md) <- c("group", "mean", "min", "max")
worldmapt <- join(x = worldmap, y = worldmap.md, by = c("group"))
worldmapt$group.regroup <- 1
worldmapt[(worldmapt$max > 180) & (worldmapt$min < 180) & (worldmapt$long.recenter >
180), c("group.regroup")] <- 2
worldmapt$group.regroup <- paste(worldmapt$group, worldmapt$group.regroup, sep = ".")
worldmap.rg <- worldmapt
# 都市地图重分组
urbanareas.mean <- aggregate(x = urbanareas[, c("long.recenter")], by = list(urbanareas$group),
FUN = mean)
urbanareas.min <- aggregate(x = urbanareas[, c("long.recenter")], by = list(urbanareas$group),
FUN = min)
urbanareas.max <- aggregate(x = urbanareas[, c("long.recenter")], by = list(urbanareas$group),
FUN = max)
urbanareas.md <- cbind(urbanareas.mean, urbanareas.min[, 2], urbanareas.max[,
2])
colnames(urbanareas.md) <- c("group", "mean", "min", "max")
urbanareast <- join(x = urbanareas, y = urbanareas.md, by = c("group"))
urbanareast$group.regroup <- 1
urbanareast[(urbanareast$max > 180) & (urbanareast$min < 180) & (urbanareast$long.recenter >
180), c("group.regroup")] <- 2
urbanareast$group.regroup <- paste(urbanareast$group, urbanareast$group.regroup,
sep = ".")
urbanareas.rg <- urbanareast
# 闭合曲线
worldmap.rg <- worldmap.rg[order(worldmap.rg$group.regroup, worldmap.rg$order),
]
worldmap.begin <- worldmap.rg[!duplicated(worldmap.rg$group.regroup), ]
worldmap.end <- worldmap.rg[c(!duplicated(worldmap.rg$group.regroup)[-1], TRUE),
]
worldmap.flag <- (worldmap.begin$long.recenter == worldmap.end$long.recenter) &
(worldmap.begin$lat == worldmap.end$lat)
table(worldmap.flag)
## worldmap.flag
## FALSE TRUE
## 13 287
worldmap.plus <- worldmap.begin[!worldmap.flag, ]
worldmap.end[!worldmap.flag, ]
## long lat order hole piece id group long.recenter
## 7962 180.00000 64.979709 56 FALSE 1 135 135.1 180.00000
## 8502 178.72530 71.098800 596 FALSE 13 135 135.13 178.72530
## 10259 -64.61101 1.328731 57 FALSE 1 170 170.1 -64.61101
## 1615 -65.40228 -11.566270 31 FALSE 1 21 21.1 294.59772
## 1734 -65.35471 1.095282 90 FALSE 1 22 22.1 294.64529
## 2088 -64.42549 45.292040 118 FALSE 1 27 27.1 -64.42549
## 2273 -64.66941 63.392927 303 FALSE 2 27 27.2 -64.66941
## 2321 -64.33400 81.927750 351 FALSE 3 27 27.3 -64.33400
## 297 -66.27334 -21.832310 109 FALSE 1 4 4.1 293.72666
## 4131 180.00000 -16.555217 16 FALSE 2 53 53.2 180.00000
## 871 -65.37133 -65.896390 542 FALSE 1 6 6.1 294.62867
## 941 -66.29003 -80.255773 612 FALSE 4 6 6.4 293.70997
## 4642 -67.15129 80.515820 119 FALSE 1 65 65.1 292.84871
## mean min max group.regroup
## 7962 91.280910 27.28818 180.0000 135.1.2
## 8502 179.525745 178.72530 180.0000 135.13.1
## 10259 152.163643 -64.89045 294.6453 170.1.1
## 1615 97.981298 -64.96489 294.6781 21.1.2
## 1734 28.176028 -64.81606 294.6781 22.1.2
## 2088 228.346420 -64.79854 294.9437 27.1.1
## 2273 255.599857 -64.86231 294.9862 27.2.1
## 2321 258.255392 -64.33400 294.5197 27.3.1
## 297 134.484462 -64.97856 294.8820 4.1.2
## 4131 179.320038 178.59684 180.0000 53.2.2
## 871 119.894717 -64.88169 294.6875 6.1.2
## 941 20.664877 -64.48813 294.2583 6.4.2
## 4642 -2.849638 -63.68925 294.6761 65.1.2
worldmap.plus$order <- worldmap.end$order[!worldmap.flag] + 1
worldmap.cp <- rbind(worldmap.rg, worldmap.plus)
worldmap.cp <- worldmap.cp[order(worldmap.cp$group.regroup, worldmap.cp$order),
]
urbanareas.rg <- urbanareas.rg[order(urbanareas.rg$group.regroup, urbanareas.rg$order),
]
urbanareas.begin <- urbanareas.rg[!duplicated(urbanareas.rg$group.regroup),
]
urbanareas.end <- urbanareas.rg[c(!duplicated(urbanareas.rg$group.regroup)[-1],
TRUE), ]
urbanareas.flag <- (urbanareas.begin$long.recenter == urbanareas.end$long.recenter) &
(urbanareas.begin$lat == urbanareas.end$lat)
table(urbanareas.flag)
## urbanareas.flag
## FALSE TRUE
## 5 12106
urbanareas.plus <- urbanareas.begin[!urbanareas.flag, ]
urbanareas.end[!urbanareas.flag, ]
## long lat order hole piece id group long.recenter
## 97438 -64.99591 -31.72799 55 FALSE 1 1022 1022.1 -64.99591
## 97619 -64.99524 -32.35196 46 FALSE 1 1025 1025.1 -64.99524
## 542025 -65.00258 -40.73363 19 FALSE 1 5386 5386.1 294.99742
## 542541 -64.98944 -42.79195 53 FALSE 1 5390 5390.1 -64.98944
## 93532 -64.98884 -21.22354 25 FALSE 1 971 971.1 -64.98884
## mean min max group.regroup
## 97438 162.35840 -64.99591 294.9884 1022.1.1
## 97619 172.58527 -64.99524 294.9857 1025.1.1
## 542025 19.29929 -64.99136 294.9981 5386.1.2
## 542541 192.11383 -64.99714 294.9963 5390.1.1
## 93532 165.62015 -64.99024 294.9979 971.1.1
urbanareas.plus$order <- urbanareas.end$order[!urbanareas.flag] + 1
urbanareas.cp <- rbind(urbanareas.rg, urbanareas.plus)
urbanareas.cp <- urbanareas.cp[order(urbanareas.cp$group.regroup, urbanareas.cp$order),
]
# 绘制背景
wrld <- geom_polygon(aes(long.recenter, lat, group = group.regroup), size = 0.1,
colour = "#090D2A", fill = "#090D2A", alpha = 1, data = worldmap.cp)
urb <- geom_polygon(aes(long.recenter, lat, group = group.regroup), size = 0.3,
color = "#FDF5E6", fill = "#FDF5E6", alpha = 1, data = urbanareas.cp)
# 绘制航线 放大系数
head(gcircles.rg)
## long lat order piece id group cnt long.recenter mean
## 1 38.79932 8.977889 1 1 1 1.1 2 38.79932 73.95761
## 2 39.38165 9.446543 2 1 1 1.1 2 39.38165 73.95761
## 3 39.96557 9.914235 3 1 1 1.1 2 39.96557 73.95761
## 4 40.55116 10.380916 4 1 1 1.1 2 40.55116 73.95761
## 5 41.13849 10.846532 5 1 1 1.1 2 41.13849 73.95761
## 6 41.72766 11.311033 6 1 1 1.1 2 41.72766 73.95761
## min max group.regroup
## 1 38.79932 116.5846 1.
## 2 38.79932 116.5846 1.
## 3 38.79932 116.5846 1.
## 4 38.79932 116.5846 1.
## 5 38.79932 116.5846 1.
## 6 38.79932 116.5846 1.
bigmap <- 0.5
airline <- geom_line(aes(long.recenter, lat, group = group.regroup, alpha = max(cnt)^0.6 *
cnt^0.4, color = 0.9 * max(cnt)^0.6 * cnt^0.4), size = 0.2 * bigmap, data = gcircles.rg)
airlinep <- geom_line(aes(long.recenter, lat, group = group.regroup, alpha = 0.04 *
max(cnt)^0.6 * cnt^0.4), color = "#FFFFFF", size = 2 * bigmap, data = gcircles.rg)
airlinepp <- geom_line(aes(long.recenter, lat, group = group.regroup, alpha = 0.02 *
max(cnt)^0.6 * cnt^0.4), color = "#ECFFFF", size = 4 * bigmap, data = gcircles.rg)
airlineppp <- geom_line(aes(long.recenter, lat, group = group.regroup, alpha = 0.01 *
max(cnt)^0.6 * cnt^0.4), color = "#ECFFFF", size = 8 * bigmap, data = gcircles.rg)
airlinepppp <- geom_line(aes(long.recenter, lat, group = group.regroup, alpha = 0.005 *
max(cnt)^0.6 * cnt^0.4), color = "#BBFFFF", size = 16 * bigmap, data = gcircles.rg)
gcircles.rg[gcircles.rg$group.regroup == 1.1, ]
## [1] long lat order piece id
## [6] group cnt long.recenter mean min
## [11] max group.regroup
## <0 rows> (or 0-length row.names)
# 画图
ggplot() + wrld + urb + airline + scale_colour_gradient(low = "#D9FFFF", high = "#ECFFFF") +
scale_alpha(range = c(0, 1)) + airlineppp + airlinepp + airlinep + theme(panel.background = element_rect(fill = "#00001C",
color = "#00001C"), panel.grid = element_blank(), axis.ticks = element_blank(),
axis.title = element_blank(), axis.text = element_blank(), legend.position = "none") +
ylim(-65, 75) + coord_equal()
绘制北京机场的国内航线图
# 提取从北京机场到全国各地的航线数据
routes_pek1 <- routes_all[routes_all$Source_airport == "PEK" & routes_all$country.y ==
"China", ]
routes_pek1$id <- row.names(routes_pek1)
# 航线进行曲线化
rts1 <- gcIntermediate(routes_pek1[, c("lon.x", "lan.x")], routes_pek1[, c("lon.y",
"lan.y")], 100, breakAtDateLine = FALSE, addStartEnd = TRUE, sp = TRUE)
rts1 <- as(rts1, "SpatialLinesDataFrame")
# 航线坐标数据
rts.pek <- fortify(rts1)
knitr::kable(head(rts.pek))
| 116.5846 |
40.08011 |
1 |
1 |
1 |
1.1 |
| 116.5898 |
39.98608 |
2 |
1 |
1 |
1.1 |
| 116.5949 |
39.89205 |
3 |
1 |
1 |
1.1 |
| 116.6001 |
39.79802 |
4 |
1 |
1 |
1.1 |
| 116.6053 |
39.70399 |
5 |
1 |
1 |
1.1 |
| 116.6104 |
39.60996 |
6 |
1 |
1 |
1.1 |
rts.ff2 <- left_join(rts.pek, routes_pek1[, c("id", "cnt")]) # 增加航线频数
knitr::kable(head(rts.ff2))
| 116.5846 |
40.08011 |
1 |
1 |
1 |
1.1 |
1 |
| 116.5898 |
39.98608 |
2 |
1 |
1 |
1.1 |
1 |
| 116.5949 |
39.89205 |
3 |
1 |
1 |
1.1 |
1 |
| 116.6001 |
39.79802 |
4 |
1 |
1 |
1.1 |
1 |
| 116.6053 |
39.70399 |
5 |
1 |
1 |
1.1 |
1 |
| 116.6104 |
39.60996 |
6 |
1 |
1 |
1.1 |
1 |
bou2_4l <- readOGR("/resources/rstudio/airport/bou2_4l.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/resources/rstudio/airport/bou2_4l.shp", layer: "bou2_4l"
## with 1785 features
## It has 8 fields
## Integer64 fields read as strings: FNODE_ TNODE_ LPOLY_ RPOLY_ BOU2_4M_ BOU2_4M_ID
bou2_4l_ff <- fortify(bou2_4l) # 数据转换
knitr::kable(head(bou2_4l_ff))
| 121.4884 |
53.33265 |
1 |
1 |
0 |
0.1 |
| 121.4974 |
53.32104 |
2 |
1 |
0 |
0.1 |
| 121.5004 |
53.31389 |
3 |
1 |
0 |
0.1 |
| 121.5316 |
53.30765 |
4 |
1 |
0 |
0.1 |
| 121.6013 |
53.26009 |
5 |
1 |
0 |
0.1 |
| 121.6599 |
53.24680 |
6 |
1 |
0 |
0.1 |
bou2_4p <- readOGR("/resources/rstudio/airport/bou2_4p.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/resources/rstudio/airport/bou2_4p.shp", layer: "bou2_4p"
## with 924 features
## It has 7 fields
## Integer64 fields read as strings: BOU2_4M_ BOU2_4M_ID
bou2_4p_ff <- fortify(bou2_4p) # 数据转换
knitr::kable(head(bou2_4p_ff))
| 121.4884 |
53.33265 |
1 |
FALSE |
1 |
0 |
0.1 |
| 121.4995 |
53.33601 |
2 |
FALSE |
1 |
0 |
0.1 |
| 121.5184 |
53.33919 |
3 |
FALSE |
1 |
0 |
0.1 |
| 121.5391 |
53.34172 |
4 |
FALSE |
1 |
0 |
0.1 |
| 121.5738 |
53.34817 |
5 |
FALSE |
1 |
0 |
0.1 |
| 121.5840 |
53.34964 |
6 |
FALSE |
1 |
0 |
0.1 |
# 城市坐标
china_cities <- world.cities[world.cities$country.etc == "China" & world.cities$pop >
5e+05, ]
head(china_cities)
## name country.etc pop lat long capital
## 1575 Anshan China 1203894 41.12 122.95 0
## 1631 Anyang China 792081 36.08 114.35 0
## 3310 Baoding China 1027984 38.87 115.48 0
## 3319 Baotou China 1301768 40.60 110.05 0
## 3826 Beijing China 7602069 39.93 116.40 1
## 3996 Bengbu China 585037 32.95 117.33 0
ggplot(bou2_4p_ff, aes(x = long, y = lat, group = group)) + geom_polygon(colour = NA,
size = 0.2, fill = "grey70") + geom_point(aes(x = long, y = lat, group = 1),
data = china_cities, color = "red", fill = "red", shape = 21, size = 1,
alpha = 0.7) + geom_line(aes(x = long, y = lat, group = group, color = group),
data = rts.ff2, size = 0.5, alpha = 0.7) + theme(panel.background = element_rect(fill = "black"),
panel.grid = element_blank(), axis.ticks = element_blank(), axis.title = element_blank(),
axis.text = element_blank(), legend.position = "none")
