1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
| setwd("C:\\Users\\jlzhang\\20180604")
#### 导入所需要的程序包library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(ggplot2)
library(rgdal)
## rgdal: version: 1.3-1, (SVN revision 747)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/jlzhang/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/jlzhang/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.2-7
rm(list = ls()) ## 删除所有对象
## 去过的地方
provinces_visited <- c("Beijing", "Hebei", "Heilongjiang",
"Liaoning", "Jilin", "Jiangxi", "Hunan",
"Zhejiang", "Hubei", "Chongqing", "Sichuan",
"Henan", "Guangxi", "Guangdong", "Fujian",
"Hainan", "Yunnan", "Guizhou")
provinces_visited_df <- data.frame(province_pinyin = provinces_visited,
visited_status = rep("Visited", length(provinces_visited)))
## 读取地图
country <- readOGR("bou1_4l.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\jlzhang\20180604\bou1_4l.shp", layer: "bou1_4l"
## with 1382 features
## It has 8 fields
## Integer64 fields read as strings: FNODE_ TNODE_ LPOLY_ RPOLY_ BOU1_4M_ BOU1_4M_ID
province <- readOGR("province_polygon.shp") ## 各省多边形
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\jlzhang\20180604\province_polygon.shp", layer: "province_polygon"
## with 34 features
## It has 10 fields
province$ID <- as.character(province$ID)
## 转换为ggplot2绘图用的data.frame
country_df <- fortify(country)
province_df <- fortify(province)
## Regions defined for each Polygons
## 提取省级数据中的信息
province_dat <- province@data
## 改正省级数据中关于香港和澳门的错误
index_HK <- which(province_dat$ID == "Xianggang")
province_dat$X[index_HK] <- 114.156121
province_dat$Y[index_HK] <- 22.37725
province_dat$ID[index_HK] <- "Hong Kong"
index_MC <- which(province_dat$ID == "Aomen")
province_dat$X[index_MC] <- 113.545681
province_dat$Y[index_MC] <- 20.197303
province_dat$ID[index_MC] <- "Macau"
## id与province_df数据表中的省出现的顺序相同， 这里id都从0开始， 所以需要减去1，以便匹配
province_dat$id <- as.character(1:nrow(province_dat) - 1)
## 增加省拼音
province_dat2 <- merge(province_dat, provinces_visited_df, by.x = "ID",
by.y = "province_pinyin", all.x = TRUE)
## 由于visited status 默认是 factor 类型， 这里需要先转换为字符串， 再处理
province_dat2$visited_status <- as.character(province_dat2$visited_status)
## 所有没有匹配上的省份都是没有去过的， 状态为 not yet
province_dat2$visited_status[is.na(province_dat2$visited_status)] <- "Not Yet"
province_df_merged <- merge(province_df, province_dat2, by = "id")
## 绘图
## 用geom_path绘制读取的polyline，因为国界线是polyline
## 多圆锥投影
## http://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/polyconic.htm
ggplot() + geom_path(data = country_df,
aes(x = long, y = lat, group=group), colour="darkblue") +
coord_map("polyconic") +
geom_polygon(data = province_df_merged,
aes(x = long, y = lat, group=group,
fill = visited_status), colour="grey") +
geom_text(mapping = aes(x = X, y = Y, label = ID),
data = province_dat, colour = 'Black')+
ggtitle("Provinces I have been to:") +
xlab("Longitude") +
ylab("Latitude") +
theme(legend.position = "bottom") +
scale_fill_discrete(name = "Province")
|
