R绘制地图原理：
R使用一个个多边形（polygon）来表示每个区域，通过顺次连接GIS数据提供的每个区域多边形的坐标来逐点绘制这些多边形，所以理论上只要得到GIS数据就可绘制相应的地图。
地图绘制说明：
- 必须使用完整的中国地图（包含台湾及南海九段线）；
- R中大部分包对于中国地图支持不佳，论文发表中需手动添加清楚明确的南海九段线。
1.maps包
R中内置地图包。但是这个包中的中国地图缺乏台湾、西沙、南沙群岛，所以不要拿来做中国地图。
library(maps)
map('usa')
map('usa',regions = "Massachusetts")
map('state',interior = F)
map('state',boundary = F,lty=2,add = T)
map('state',fill=T,col = rainbow(200))
map('county','washington, san',names=T,plot = F)
map('state',xlim = range(ozone$x),ylim = range(ozone$y))
#经纬度范围
text(ozone$x,ozone$y,ozone$median)
box()
image.png
#标记每个区域失业率
data(unemp)
data(county.fips)
colors=c("#F1EEF6",'#D4B9DA','#C994C7','#DF65B0','#DD1C77','#980043')
unemp$colorBuckets <- as.numeric(cut(unemp$unemp,c(0,2,4,6,8,10,100)))
leg.txt <- c('<2%','2-4%','6-8%','8-10%','>10%')
colorsamtched <- unemp$colorBuckets[match(county.fips$fips,unemp$fips)]
map('county',col=colors[colorsamtched],fill=T,resolution = 0,lty=0)
#调用projection映射到专业地图坐标系
map('county',col=colors[colorsamtched],fill=T,resolution = 0,lty=0,
projection = 'polyconic')
#添加图例
leg.txt <- c('<2%','2-4%','6-8%','8-10%','>10%')
legend('topright',leg.txt,horiz = F,fill = colors)
title('2009年美国各市失业情况')
image.png
2. 从本地导入GIS地图
mapdata包和maptools包导入GIS地图数据。
library(mapdata)
library(maps)
map('china')
image.png
虽然这个包的中国地图版图完整，但提供的GIS信息比较旧，如重庆和四川没分开等。
maptools包可读入shapefile格式数据（最常见的地图文件格式）。
中国地图GIS数据可从GADM数据库下载（https://gadm.org/download_country_v3.html）如下图，选择shapefile或R格式：
image.png
或从国家基础地理信息中心（http://www.ngcc.cn/ngcc/
）下载shapefile文件（但我没找到）。
library(sp)
#sp包提供了处理空间数据的类和方法
library(maptools)
library(plyr)
library(ggplot2)
#读入地图数据
china.map <- readShapePoly('*.shp')
class(china.map)
plot(china.map) #仍在笛卡尔坐标系（扁平）
ggplot(china.map,aes(x=long,y=lat,group=group))+
geom_polygon(fill='grey')+
coord_map('polyconic') #投影到专业地图坐标系
3. 利用ggplot2绘制地图
准备地图信息。
#将地图数据转换为数据框（经纬度数据）
china.map1 <- fortify(china.map)
#提取行政区域信息
mymap <- china.map@data #@针对S4对象，类似$（S3对象）
#生成id用于匹配
mymap$id <- c(0:924)
#按id列匹配
china.map2 <- plyr::join(mymap,china.map1)
导入业务数据和绘制地图：
如中国各省份公共卫生机构分布数据以及某疾病发病率。
geom_polygon(color='grey')+scale_fill_gradient(low = 'pink',
high = 'red')+
coord_map('ployconic')
p
#给各省份添加标签，仍需获取省份的经纬度
#可通过geocode/getCoordinate从地图软件中获取
p+geom_text(aes(label=Group.1),size=3,data=province_data)+
#去掉灰色网格背景和不必要的经纬度信息
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank())
还可在地图上添加气泡图、饼图等。
除了绘制中国和世界地图，可从shape文件中获取市级地图GIS数据。
anhui <- subset(china.map2,NAME=="安徽省")
ggplot(anhui,aes(x=long,y=lat,group=group,fill=NAME))+
geom_polygon(color='skyblue',fill='skyblue')
要获取安徽省各个行政区的资料信息，可从GADM数据库中下载包含市级行政区的shapefile文件。
4. 从专业地图软件调用地图
R可调用谷歌、百度等地图数据，但每日API获取次数有限。
一般流程是：
- 获取经纬度
- 根据经纬度获得相应地图
- 按需处理数据和地图对象
