最近，微信中去过哪些省份的小程序十分火爆，用户通过选择去过哪些省份和地级市， 提交之后后台就自动生成地图。其实该地图要绘制起来并不困难。 这里给出绘制去过哪些省份的R代码。供参考

代码和R脚本 map of china.zip

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")

参考
## https://zhuanlan.zhihu.com/p/25231546 
## http://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html 
## http://www.statsoft.org/wp-content/uploads/2016/09/Lecture6_HKMapVis.html