欢迎关注天善智能，我们是专注于商业智能BI，人工智能AI，大数据分析与挖掘领域的垂直社区，学习，问答、求职一站式搞定！
对商业智能BI、大数据分析挖掘、机器学习，python，R等数据领域感兴趣的同学加tstoutiao，邀请你进入数据爱好者交流群，数据爱好者们都在这儿。
作者：李誉辉
四川大学在读研究生
1.简介
此次教程是关于使用ggplot2包进行地理信息可视化的，
笔者将从头到尾，条分缕析的梳理ggplot2地图可视化的知识点。
让读者对ggplot2地图可视化的过程有更深入的理解。
并可以作为知识库进行查询。
ggplot2支持2种地理数据模型：
* sp, sp对象，全称SpatialPolygonsDataFrame, 其数据类型为数据框，
使用函数geom_map()/geom_polygon() coord_map()绘制。
* sf，sf对象，全称Simple feature,
使用函数geom_sf()/stat_sf() coord_sf()绘制。
若要增加文本注释，则使用geom_sf_label(),geom_sf_text()。
2.sp数据类型
2.1 geom_map()地图对象
geom_map(mapping = NULL, data = NULL, stat = "identity", ..., map, na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
关键参数：
- mapping, 表示要在地图上映射的美学特征，用aes()指定。
- data, 表示要展示在地图上的美学特征的数据集。为数据框格式，无指定则从前面继承。
- stat, 表示对美学特征的统计变换。
- 更多统计变换
- ..., 表示其它美学特征参数。如colour = "red",或size = 3指定border参数。
- map, 表示地图数据集(构成地图行政区域多边形边界的坐标数据),
- 为数据框格式。必须包含至少3个变量:
- x坐标或long(经度)，y坐标或lat(纬度)，region或id(行政区域多边形编号)。
- na.rm, 表示是否移除缺失值，默认FALSE则出现warning，TRUE则直接移除不显示warning。
- show.lengend,表示是否显示美学特征的图例，默认NA表示有美学特征则显示，
- FALSE表示不显示，TRUE表示一直显示。
- 如果有多种美学特征，可以用逻辑向量来设定要显示的几种的美学特征的图例。
- inherit.aes, 表示是否继承前面的美学特征，FALSE表示不继承；
- TRUE则表示将前面的几何特征与现在的几何特征结合起来，一起显示。
美学特征：
geom_map()支持多种美学特征，相比其它geom_xxx()对象，增加了map-id特征。
常见美学特征： map-id, alpha, colour, fill, group, linetype, size。
几何特征详细使用方法
：
R_ggplot2基础（一）
2.1.1 数据集构成
library(ggplot2) # 编造数据集 ids <- factor(c("1.1", "2.1", "1.2", "2.2", "1.3", "2.3")) positions <- data.frame( id = rep(ids, each = 4), # each=4表示每个多边形由4个点构成，即为四边形 x = c(2, 1, 1.1, 2.2, 1, 0, 0.3, 1.1, 2.2, 1.1, 1.2, 2.5, 1.1, 0.3, 0.5, 1.2, 2.5, 1.2, 1.3, 2.7, 1.2, 0.5, 0.6, 1.3), y = c(-0.5, 0, 1, 0.5, 0, 0.5, 1.5, 1, 0.5, 1, 2.1, 1.7, 1, 1.5, 2.2, 2.1, 1.7, 2.1, 3.2, 2.8, 2.1, 2.2, 3.3, 3.2) ) values <- data.frame( id = ids, value = c(3, 3.1, 3.1, 3.2, 3.15, 3.5) )
2.1.2 地图可视化初步
library(ggplot2) ggplot(values, aes(fill = value)) # 按每个id给多边形区域填充颜色，色值与value变量匹配 geom_map(aes(map_id = id), map = positions, colour = "magenta", size = 1) expand_limits(positions) # 扩展映射参数的显示范围，不能缩小，只能扩展
调整坐标轴显示范围
library(ggplot2) ggplot(values, aes(fill = value)) geom_map(aes(map_id = id), map = positions, colour = "magenta", size = 1) expand_limits(positions) ylim(0, 3)
2.1.3 map参数
通过map参数指定数据集，画国家地图，世界地图。
library(ggplot2) data("USArrests") # 美国犯罪率数据集 crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests) crimesm <- reshape2::melt(crimes, id = 1) # 宽转长 # 绘制地图 states_map <- map_data("state") # ggplot2内置州级地图数据集 USA_map <- ggplot(crimes, aes(map_id = state)) geom_map(aes(fill=Murder), map=states_map,colour="cyan",size = 0.5) # 指定map参数 scale_fill_distiller(palette = "RdPu",direction= 1) # 修改颜色标度 expand_limits(x = states_map$long, y = states_map$lat) # 扩大显示，显示所有在经纬度内的 USA_map
地图分面:
library(ggplot2) ggplot(crimesm, aes(map_id = state)) geom_map(aes(fill = value), map = states_map,colour = "cyan", size=0.5) scale_fill_distiller(palette = "RdPu", direction= 1) expand_limits(x = states_map$long, y = states_map$lat) facet_wrap(vars(variable)) # 地图分面，按变量variable封装分面
2.2 coord_map()地图投影变换
注意：coord_map()现在很多时候已经被geom_sf coord_sf() 替代了。
而且后者支持更多的地图投影类型。但是地图投影原理是一样的。
地图投影变换是指将球面地图投影在平面上，而geom_map本质上是笛卡尔坐标系下的多边形，
coord_map()变换非常消耗计算资源。变换时最好清空内存：rm(list = ls()); gc()。
语法：
coord_map(projection = "mercator", ..., parameters = NULL, orientation = NULL, xlim = NULL, ylim = NULL, clip = "on") coord_quickmap(xlim = NULL, ylim = NULL, expand = TRUE, clip = "on")
关键参数：
- projection, 表示指定地图投影类型，有好几十种，详情见?mapproj::mapproject()。
- ...,parameters,其它传递给mapproj::mapproject()的参数。
- ...表示已命名的参数；parameters表示未命名的参数，指定parameters参数则忽略...。 parameters = c(lat, lon, rot),分别表示投影中心点的经纬度和。
- 所有的标准地图投影都是基于本初子午线(the Prime Meridian)的。
- orientation, 表示指定投影方向，默认为c(90, 0, mean(range(x)))，
- 第1个数字为纬度，北纬为正，南纬为负。
- 第2个数字为经度，东经为正，西经为负。
- 第3个数字表示旋转角度，顺时针旋转为正。
- xlim/ylim, 表示手动指定x/y轴的边界(对应经度/纬度)。
- clip, 表示是否修剪图层以适应绘图面板的尺寸。
- 默认“on”表示修剪，“off”则表示否。绝大多数情况下，不需要修剪。
- “off”可能造成在绘图面板的margins上出现几何对象。
- expand, 为逻辑值，默认TRUE表示会增加一个扩展因子，使得坐标轴不与坐标轴重叠，
- 很多时候我们发现0刻度不在x轴或y轴，会有点偏移，就是这个原因。
对于任何地图投影，其经度方向的比例尺和纬度方向的比例尺都是变化的。
这也是运算量大的原因之一。
coord_quickmap()就是固定纵横比，即固定lat/long。这样相对墨卡托地图投影更快，
但对于其它地图投影，并不能代替。
2.2.1 projection类型
常用地图投影：
赤道投影：
赤道投影的中心位于本初子午线上，其纬线为水平的平行线(不一定是直线)。 赤道投影在赤道附近的精确度非常好。
- mercator,墨卡托投影
- 特征是：经纬线都是直线，分别相互平行，不能显示南北极极点。
- 此投影形状为等角投影。由于该投影维持局部角度关系不变，所以能很好地描绘微小形状。
- 大面积变形使得墨卡托投影不适用于常规地理世界地图。 常用于绘制航空旅行，风向，洋流等地图。
- sinusoidal, 正弦曲线投影
- 特征是：经线是基于正弦函数的曲线，纬线是直线。
- 用于世界地图时，无论是否存在等角变形，它都能保持等积。
- mollweide, 摩尔维特投影
- 特征： 此投影也称为巴比涅投影、椭圆投影或等面积投影。
- 所有纬线都是直线，所有经线都是等间距的椭圆弧。
- 唯一例外的是中央子午线，中央子午线是直线。极点是点。
- 适用于绘制世界专题或分布地图，经常采用不连续的形式。
极方位投影：
对于极方位投影，经线沿线的距离是精确的，但纬线圆沿线的变形从中心向外增大 极方位投影最适用于 30° 纬线半径范围内的区域，因为这些区域的变形程度最小。
极方位投影在南北极附近的精度特别好。
- azequidistant,等距方位投影，
- 特征： 这种投影会将所有经线和纬线划分为相等的部分，以保持等距离属性。 以某一城市为中心的斜轴投影也被普遍使用。
- 常用于空中导航和海上导航的路线，
- 这些地图将定位到一个重要的位置作为中心点并使用相应的投影方法。
- azequalarea,亚尔勃斯等积圆锥投影
- 特征：经线是相交于一个公共点的间距相等的直线。
- 极点表示为弧，而不表示为单个点。
- 纬线是间距不等的同心圆，距离极点越近，同心圆的间距越小。 这种投影最适合于东西方向分布的大陆板块或国家，而不适合南北方向分布的大陆板块。
- gnomonic,球心投影
- 特征：此投影是一种从地球中心进行观察的平面透视投影。可以从任意方位进行投影。
- 无论方位如何，所有的大圆都将投影为直线。
- 此投影适用于导航，因为大圆将突出显示距离最短的路线。
- orthographic, 正射投影
- 特征： 此透视投影从无穷远处观察地球。这样便可提供地球的三维图像。
- 面积比例尺从中心开始随距离的增加而减小。在半球边缘处面积比例尺为零。 此投影多用于美观的展示图而不是技术应用。
圆锥投影：
关于本初子午线对称的极地投影。
除bonne投影外，其它经线都是垂直与同心圆环。
- simpleconic,简单圆锥投影
- 特征：所有圆形纬线彼此等距，即在经线方向上间距相等。
- 无论使用一条纬线还是两条纬线，情况都是如此。
- 对主要沿东西方向延伸的中纬度地区进行区域制图。
- 普遍适用于小国家/地区的地图集地图。
- lambert,兰勃特等角圆锥投影
- 特征：此投影是最适用于中纬度的一种投影。其描绘形状比描绘面积更准确。
- 美国国家平面坐标系对所有具有较大东西范围的区域均使用此投影。
投影地图图片：
墨卡托投影
正弦曲线投影
高尔立体投影
摩尔维特投影
等距方位投影
亚尔勃斯等积圆锥投影
球心投影
简单圆锥投影
兰勃特等角圆锥投影
2.2.2 地图投影变换初步
library(ggplot2) USA_map coord_map() # 默认mercator投影 # USA_map coord_map(projection = "mercator") # 结果与上面一样 USA_map coord_quickmap() # 相比上面的，有少许失真
2.2.3 地图投影深入
接下来，我们将针对不同的地图投影，使用ggplot2包内置数据集，绘制投影地图。
ggplot2内置地图数据集列表：(“world”,“usa”,“state”,“county”,“nz”)
2.2.3.1 世界地图
rm(list = ls()); gc() # 清空内存 library(ggplot2) world_data <- map_data("world") # ggplot2内的数据集 # 使用geom_polygon结果一样 world_map <- ggplot(world_data) geom_polygon(aes(x = long, y = lat, group = group, fill = region), colour = "cyan", size = 0.5, show.legend = FALSE) world_map coord_map() # 默认墨卡托地图投影 world_map coord_map(projection = "mollweide") # 摩尔维特投影，这里不能增加parameters参数 world_map coord_map(projection = "sinusoidal") # 正弦曲线投影
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 795072 42.5 1628671 87 1628671 87 ## Vcells 2339195 17.9 8388608 64 8384404 64
2.2.3.2 局部区域地图
rm(list = ls()); gc() # 清空内存 library(ggplot2) county_data <- map_data("county") county_map <- ggplot(county_data) geom_polygon(aes(x = long, y = lat, group = group, fill = region), colour = "cyan", size = 0.5, show.legend = FALSE) county_map coord_map("gnomonic") # 球心投影 county_map coord_map("azequidistant") # 等距方位投影 county_map coord_map("azequalarea") # 亚尔勃斯等积圆锥投影 county_map coord_map("bonne", lat0 = 50) # 彭纳投影
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 790580 42.3 1698867 90.8 1698867 90.8 ## Vcells 2109747 16.1 10146329 77.5 10146329 77.5
2.2.3.3 半球地图
rm(list = ls()); gc() # 清空内存 library(ggplot2) # 世界地图， world <- map_data("world") worldmap <- ggplot(world, aes(x = long, y = lat, group = group)) geom_polygon(fill = "cyan",color = "magenta") scale_y_continuous(breaks = (-2:2) * 30) scale_x_continuous(breaks = (-4:4) * 45) # 正射投影 worldmap coord_map("ortho") # 默认北极为中心点。 # 南极为中心点 worldmap coord_map("ortho", orientation = c(-90, 0, 0)) worldmap coord_map("ortho", orientation = c(-90, 0, 30)) # 顺时针旋转30度 worldmap coord_map("ortho", orientation = c(41, -74, 0)) # 随便取一点
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 791366 42.3 1744096 93.2 1744096 93.2 ## Vcells 2042362 15.6 10146329 77.5 10146329 77.5
3.sf数据模型
3.1 geom_sf()以美国人口调查数据为例，数据来源（https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html）,需要科学上网。
rm(list = ls()); gc() # 清空内存 library(ggplot2) path1 <- "E:/R_input_output/data_input/cb_2013_us_state_20m/cb_2013_us_state_20m.shp" USA_cencus <- sf::st_read(dsn = path1) ggplot(USA_cencus) geom_sf(aes(fill = NAME), show.legend = FALSE) # 按州的名字填充颜色,图例过长不显示 coord_sf(xlim = c(-170, -60)) # 设定显示范围
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 791329 42.3 1744096 93.2 1744096 93.2 ## Vcells 2110470 16.2 12255594 93.6 12236244 93.4 ## Reading layer `cb_2013_us_state_20m' from data source `E:\R_input_output\data_input\cb_2013_us_state_20m\cb_2013_us_state_20m.shp' using driver `ESRI Shapefile' ## Simple feature collection with 52 features and 9 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -179.1473 ymin: 17.88481 xmax: 179.7785 ymax: 71.35256 ## epsg (SRID): 4269 ## proj4string: proj=longlat datum=NAD83 no_defs
rm(list = ls()); gc() # 清空内存 library(ggplot2) path1 <- "E:/R_input_output/data_input/cb_2013_us_state_20m/cb_2013_us_state_20m.shp" USA_cencus <- sf::st_read(dsn = path1) ggplot(USA_cencus) geom_sf(aes(fill = ALAND),color = "cyan", size = 0.5) # 按人口数据填充颜色, coord_sf(xlim = c(-170, -60)) # 设定显示范围 scale_fill_distiller(palette = "RdPu",direction = -1) # 调整颜色标度
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 900291 48.1 1744096 93.2 1744096 93.2 ## Vcells 1786205 13.7 9804475 74.9 12236244 93.4 ## Reading layer `cb_2013_us_state_20m' from data source `E:\R_input_output\data_input\cb_2013_us_state_20m\cb_2013_us_state_20m.shp' using driver `ESRI Shapefile' ## Simple feature collection with 52 features and 9 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -179.1473 ymin: 17.88481 xmax: 179.7785 ymax: 71.35256 ## epsg (SRID): 4269 ## proj4string: proj=longlat datum=NAD83 no_defs
3.2 coord_sf()
coord_sf()中有参数crs用于指定投影方式，
需使用proj4string格式的字符串。如墨卡托投影用crs = " proj=merc"。
某些地图投影其字符串中包含经纬度等参数，如crs = " proj=laea lat_0=35 lon_0=-100"。
其它地图投影见Projection methods（https://proj4.org/operations/projections/index.html）。
crs本身包含坐标数据的，就不能再增加xlim, ylim等参数。
rm(list = ls()); gc() # 清空内存 library(ggplot2) path1 <- "E:/R_input_output/data_input/cb_2013_us_state_20m/cb_2013_us_state_20m.shp" USA_cencus <- sf::st_read(dsn = path1) ggplot(USA_cencus) geom_sf(aes(fill = NAME), show.legend = FALSE) # 按州的名字填充颜色,图例过长不显示 coord_sf(crs = " proj=aea lat_1=25 lat_2=50 lon_0=-100") # 设定地图投影 ggtitle("Albers equal-area projection") # 指定标题
## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 900543 48.1 1744096 93.2 1744096 93.2 ## Vcells 1786455 13.7 9804475 74.9 12236244 93.4 ## Reading layer `cb_2013_us_state_20m' from data source `E:\R_input_output\data_input\cb_2013_us_state_20m\cb_2013_us_state_20m.shp' using driver `ESRI Shapefile' ## Simple feature collection with 52 features and 9 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -179.1473 ymin: 17.88481 xmax: 179.7785 ymax: 71.35256 ## epsg (SRID): 4269 ## proj4string: proj=longlat datum=NAD83 no_defs
往期精彩：
- R_插值_拟合_回归_样条
- R_circlize包_和弦图（一）
- R_ggplot2基础（一）
- R_3D图(一)
公众号后台回复关键字即可学习
回复 爬虫 爬虫三大案例实战
回复 Python 1小时破冰入门
回复 数据挖掘 R语言入门及数据挖掘
回复 人工智能 三个月入门人工智能
回复 数据分析师 数据分析师成长之路
回复 机器学习 机器学习的商业应用
回复 数据科学 数据科学实战
回复 常用算法 常用数据挖掘算法
万水千山总是情，点个 “好看” 行不行！！！
