玩转数据可视化之R语言ggplot2：（八）ggplot2绘制空间地理数据图
玩转数据可视化之R语言ggplot2
- 🌸个人主页：JoJo的数据分析历险记
- 📝个人介绍：小编大四统计在读，目前保研到统计学top3高校继续攻读统计研究生
- 💌如果文章对你有帮助，欢迎关注、点赞、收藏、订阅专栏
本系列主要介绍R语言ggplot2的使用
参考资料:
ggplot2: Elegant Graphics for Data Analysis
文章目录
🌲8.ggplot2绘制空间地理数据图
绘制地理数据图在我们研究一些空间统计问题时很有效。例如绘制一张全国疫情情况热力图。现在用的比较多的是Arcgis软件绘制地理数据图，也可以使用python，python数据可视化未来我也会写相应的专栏。如果有精力的话，考虑再出一个Arcgis的专栏绘制空间地理数据图。言归正传，在ggplot2中，绘制空间地理数据图主要分为两个步骤：
- 第一个是使用相应地理数据先绘制地图
- 第二个是使用其他相关的数据在地图上添加信息
接下来我们来看看具体是如何实现这两个步骤的
# 导入基本库
library(ggplot2)
library(tidyverse)# 数据操作处理库，具体可以看我的R语言数据分析专栏
🌳8.1 使用多边形geom_polygon()绘制地图
在R中，绘制地图最简单的方法就是通过geom_polygon()绘制不同区域的边界。我们使用ggplot2中map_data()里的数据
我们在这里以美国密歇根州数据为例进行地图绘制，首先我们来看一下这个数据
mi_county <- map_data('county','michigan')#筛选出密西根州的地理数据
head(mi_county)
| long | lat | group | order | region | subregion | |
|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <int> | <chr> | <chr> | |
| 1 | -83.88675 | 44.85686 | 1 | 1 | michigan | alcona |
| 2 | -83.36536 | 44.86832 | 1 | 2 | michigan | alcona |
| 3 | -83.36536 | 44.86832 | 1 | 3 | michigan | alcona |
| 4 | -83.33098 | 44.83968 | 1 | 4 | michigan | alcona |
| 5 | -83.30806 | 44.80530 | 1 | 5 | michigan | alcona |
| 6 | -83.30233 | 44.77665 | 1 | 6 | michigan | alcona |
可以看出上述数据包含六列，long和lat分别表示经纬度，group表示子区域的分组，subregion表示所在的子区域名称
# 将subregion修改名为id
colnames(mi_county)[6] <- 'id'#修改列名
下面我们使用geom_point和geom_polygon()来更仔细的观察一下这个数据集
ggplot(mi_county, aes(long, lat)) +
geom_point(size = .25, show.legend = FALSE) +
coord_quickmap()
ggplot(mi_county, aes(long, lat, group = group)) +
geom_polygon(fill='white',col='grey50')+
coord_quickmap()
geom_point()和geom_polygon()在我们之前的内容中都介绍过，但是这里要注意的是，
在我们画多边形图的时候，要注意定义分组变量，否则多边形的连接会出现交叉混乱，这个很好理解。
还有我们在两个图中都使用了coord_quickmap()这是确保他们的经纬度一致。
虽然我们可以使用上述方法绘制地图，但是这个方法有一些局限性就是我们现实世界中，很少有这样的数据，比较多的往往是矢量数据。下面我们介绍如何使用矢量数据进行绘图
🌴8.2 sf包绘制地图
为了处理矢量地理数据，我们之前使用sp包进行分析，但是sp包的兼容性较差，因此我们现在使用sf包来绘制地图，
主要有geom_sf()和coord_sf()函数来配合sf包一起使用。
接下来我们使用ozmap包里面的矢量数据来介绍sf是如何工作的。
# 导入库
library(ozmaps)
library(sf)
oz_states <- ozmap_states
oz_states
Warning message:
"package 'ozmaps' was built under R version 4.1.3"
Warning message:
"package 'sf' was built under R version 4.0.5"
Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
Registered S3 method overwritten by 'geojsonsf':
method from
print.geojson geojson
| NAME | geometry | |
|---|---|---|
| <chr> | <MULTIPOLYGON [arc_degree]> | |
| 1 | New South Wales | MULTIPOLYGON (((150.7016 -3... |
| 2 | Victoria | MULTIPOLYGON (((146.6196 -3... |
| 3 | Queensland | MULTIPOLYGON (((148.8473 -2... |
| 4 | South Australia | MULTIPOLYGON (((137.3481 -3... |
| 5 | Western Australia | MULTIPOLYGON (((126.3868 -1... |
| 6 | Tasmania | MULTIPOLYGON (((147.8397 -4... |
| 7 | Northern Territory | MULTIPOLYGON (((136.3669 -1... |
| 8 | Australian Capital Territory | MULTIPOLYGON (((149.2317 -3... |
| 9 | Other Territories | MULTIPOLYGON (((167.9333 -2... |
ozmaps包包含了澳大利亚各州边界、地方政府区域、选举边界等数据。
上述得到了澳大利亚领土边界
这个输出显示了与数据相关的，并告诉我们数据本质上是一个包含9行2列的TIBLE。
sf数据的一个优势是显而易见的，我们可以很容易地看到数据的整体结构：澳大利亚由六个州和一些地区组成。有9个不同的地理单元，所以这个TIBLE中有9行。
最重要的列是geometry，表示每个州或地区的空间几何对象。其中每一个元素都是一个多元的多边形。从字面意思来看，这个元素包含了多边形多个顶点，对于这种数据，我们可以使用geom_sf()和coord_sf()直接绘制,甚至我们都不需要指定特定的美学特征aes()，具体如下：
ggplot(oz_states) +
geom_sf()+
coord_sf()
为了理解为什么这样是work的，我们来解释一下geom_sf()
- 1.对于上面这个数据集，我们不指定任何的aes,因为geom_sf()会自动去匹配
geometry - 2.如果数据是一个sf对象，geom_sf()会自动判断一个geometry列，即使列名不是geometry
- 3.我们也可以自己设置aes(geomtry = my_column)，当数据有很多geometry对象时
🖊️8.2.1 对地图分层
在某些情况下，我们可能希望将一张地图覆盖在另一张地图上。
ggplot2软件包支持允许将多个geom_sf()层添加到图形。
我将使用oz_state数据以不同的颜色绘制澳大利亚各个州，并将此图与澳大利亚选举区域的边界重叠。
为了节省绘图时间，这里使用rmapshaper::ms_simply包来简化原始数据。下面来看一下具体如何实现
# 选举边界
oz_votes <- rmapshaper::ms_simplify(abs_ced)
ggplot() +
#绘制州边界
geom_sf(data = oz_states, aes(fill = NAME), show.legend = FALSE) +
#绘制选举区域边界
geom_sf(data = oz_votes, fill = NA) +
coord_sf()
Registered S3 method overwritten by 'geojsonlint':
method from
print.location dplyr
我们使用fill来填充颜色，本例数据我们使用的是name的分类型变量，它并没有传递任何其他的信息，但是这种方法可以为我们提高一个信息就是我们可以使用fill来区分不同区域某一变量的情况，例如假设oz_state有人口这一数据，我们可以以人口来fill
🖌️8.2.2 对地图添加标签
上一章我们介绍了如何给图形添加注释，这一章我们来看一下如何给地图添加标签。我们可以使用geom_sf_label()
和geom_sf_text()，例如我们想要对悉尼地区个选举地域增加标签。我们可以使用geom_sf_label()显示各个选民区域的名字
# 筛选出悉尼大都会的选民区域
sydney_map <- ozmaps::abs_ced %>% filter(NAME %in% c(
"Sydney", "Wentworth", "Warringah", "Kingsford Smith", "Grayndler", "Lowe",
"North Sydney", "Barton", "Bradfield", "Banks", "Blaxland", "Reid",
"Watson", "Fowler", "Werriwa", "Prospect", "Parramatta", "Bennelong",
"Mackellar", "Greenway", "Mitchell", "Chifley", "McMahon"
))
# 绘制悉尼大都会选民区域图，并增加标签
ggplot(sydney_map) +
geom_sf(aes(fill = NAME), show.legend = FALSE) +
coord_sf(xlim = c(150.97, 151.3), ylim = c(-33.98, -33.79)) +
geom_sf_label(aes(label = NAME), label.padding = unit(1, "mm"))
Warning message in st_point_on_surface.sfc(sf::st_zm(x)):
"st_point_on_surface may not give correct results for longitude/latitude data"
🖍️8.2.3 在地图上绘制其他图形
虽然geom_sf()是一种特殊的绘图，但是和其他几何图形一样，允许通过层级的方法添加其他元素。
相当于geom_sf()后依然是一个geoms标准对象，我们可以做所有之前的处理。例如我们想要在某些地方标记点
oz_capitals <- tibble::tribble(
~city, ~lat, ~lon,
"Sydney", -33.8688, 151.2093,
"Melbourne", -37.8136, 144.9631,
"Brisbane", -27.4698, 153.0251,
"Adelaide", -34.9285, 138.6007,
"Perth", -31.9505, 115.8605,
"Hobart", -42.8821, 147.3272,
"Canberra", -35.2809, 149.1300,
"Darwin", -12.4634, 130.8456,
)
ggplot() +
geom_sf(data = oz_votes) +
geom_sf(data = oz_states, fill = NA) +
geom_point(data = oz_capitals, aes(x = lon, y = lat), colour = "red") +
coord_sf()
🌵8.3 地图投影
在上面的地图绘制中，我们默认的选择了笛卡尔平面，粗略来看，这是可以的，但如果我们关心准确性，这还不够好。这种方法有两个基本问题。
- 地球的形状
关于地球形状的一系列假设被称为大地基准面，虽然这对某些数据可视化可能并不重要，但对其他数据可视化来说则至关重要。
例如如果我们要研究北美区域，使用(NAD83),如果要研究全球区域，使用(WGS84)可能较好
- 地图的形状
地球是近似椭球的，但在大多数情况下，我们的空间地理数据需要绘制在二维平面上。要将椭球体的表面映射到一个平面，我们需要进行一些变形，这就是地图投影的工作。
总的来说，大地基准面（如WGS84）、地图投影类型和投影参数（如原点位置）指定了一个坐标参考系（CRS），这是一套完整的假设，用于将纬度和经度信息转换为二维地图。sf对象通常包括一个默认的CRS，我们可以使用st_crs()来访问
st_crs(oz_votes)
Coordinate Reference System:
User input: EPSG:4283
wkt:
GEOGCRS["GDA94",
DATUM["Geocentric Datum of Australia 1994",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Australia including Lord Howe Island, Macquarie Islands, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
BBOX[-60.56,93.41,-8.47,173.35]],
ID["EPSG",4283]]
上面返回了一个WKT字符串，可以看出oz_states默认的crs是4283，我们直接通过看EPSG代码来却数据的crs
st_crs(oz_votes) == st_crs(4283)
TRUE
在ggplot2中，CRS由coord_sf（）控制，以确保绘图中的每一层使用相同的投影。默认情况下，coord_sf（）使用与数据的geometry列关联的CRS。因为sf数据会自动提供一个合理的CRS选择，这个过程通常是不可见的，不需要用户干预。但是，如果需要自己设置CRS，可以通过将有效的用户输入传递给st_CRS（）来指定CRS参数。下面的示例说明了如何从默认的CRS切换到EPSG代码3112
ggplot(oz_votes) + geom_sf()
ggplot(oz_votes) + geom_sf() + coord_sf(crs = st_crs(3112))
🌾8.4 导入shp文件分析
一般我们在实际分析时，要根据实际问题得到相应的shp矢量文件，在sf中，可以使用read_sf()读取文件。这里我以海南省十八个市县的地理边界数据为例导入
hainan <- read_sf("C:\\Users\\DELL\\Desktop\\数据\\十八个市县.shp")
hainan
| 省 | 省代码 | 市 | 市代码 | 市类型 | 省类型 | geometry |
|---|---|---|---|---|---|---|
| <chr> | <dbl> | <chr> | <dbl> | <chr> | <chr> | <MULTIPOLYGON [arc_degree]> |
| 海南省 | 460000 | 白沙县 | 469025 | 省直辖县 | 省 | MULTIPOLYGON (((109.2941 18... |
| 海南省 | 460000 | 保亭县 | 469029 | 省直辖县 | 省 | MULTIPOLYGON (((109.5084 18... |
| 海南省 | 460000 | 昌江县 | 469026 | 省直辖县 | 省 | MULTIPOLYGON (((109.0407 19... |
| 海南省 | 460000 | 澄迈县 | 469023 | 省直辖县 | 省 | MULTIPOLYGON (((110.1328 19... |
| 海南省 | 460000 | 儋州市 | 460400 | 地级市 | 省 | MULTIPOLYGON (((109.4279 19... |
| 海南省 | 460000 | 定安县 | 469021 | 省直辖县 | 省 | MULTIPOLYGON (((110.2206 19... |
| 海南省 | 460000 | 东方市 | 469007 | 省直辖县 | 省 | MULTIPOLYGON (((109.121 18.... |
| 海南省 | 460000 | 海口市 | 460100 | 地级市 | 省 | MULTIPOLYGON (((110.5089 20... |
| 海南省 | 460000 | 乐东县 | 469027 | 省直辖县 | 省 | MULTIPOLYGON (((109.4009 18... |
| 海南省 | 460000 | 临高县 | 469024 | 省直辖县 | 省 | MULTIPOLYGON (((109.851 19.... |
| 海南省 | 460000 | 陵水县 | 469028 | 省直辖县 | 省 | MULTIPOLYGON (((110.0107 18... |
| 海南省 | 460000 | 琼海市 | 469002 | 省直辖县 | 省 | MULTIPOLYGON (((110.6064 19... |
| 海南省 | 460000 | 琼中县 | 469030 | 省直辖县 | 省 | MULTIPOLYGON (((109.8724 19... |
| 海南省 | 460000 | 三亚市 | 460200 | 地级市 | 省 | MULTIPOLYGON (((109.4009 18... |
| 海南省 | 460000 | 屯昌县 | 469022 | 省直辖县 | 省 | MULTIPOLYGON (((110.1109 19... |
| 海南省 | 460000 | 万宁市 | 469006 | 省直辖县 | 省 | MULTIPOLYGON (((110.1542 19... |
| 海南省 | 460000 | 文昌市 | 469005 | 省直辖县 | 省 | MULTIPOLYGON (((111.2004 19... |
| 海南省 | 460000 | 五指山市 | 469001 | 省直辖县 | 省 | MULTIPOLYGON (((109.7113 18... |
我们可以在这个数据基础进行相应的操作，后续考虑出一些使用矢量数据实战分析的文章,下面我们先简单看一下该地理数据的地图
ggplot(hainan)+
geom_sf(aes(fill=市),show.legend = FALSE)+
coord_sf()
💟文章推荐
如果想了解更多ggplot2数据可视化技巧，欢迎访问下列文章
🌝玩转数据可视化之R语言ggplot2：（三）ggplot2实现将多张图放在一起，包括并排和插图绘制（快速入门）
🌜玩转数据可视化之R语言ggplot2:：（四）单一基础几何图形绘制
🌟玩转数据可视化之R语言ggplot2：（五）分组画图
☀️玩转数据可视化之R语言ggplot2：（六）统计变换绘图：包括加权绘图、数据分布图、曲面图、图形重叠处理等
⭐玩转数据可视化之R语言ggplot2：（七）对图形添加注释和标签（包含标题、坐标轴、参考线和高亮等注释方法）
