使用 R 语言进行地理计算:入门

使用 R 语言进行地理计算:入门

这是《R 语言地理数据科学》系列课程的第一讲,通常这种课程是面向 R 语言高级用户的,但是考虑到大家的 R 语言可能还没入门,所以我打算尽可能细致的讲解,让这门课即使是刚入门的 R 用户也能轻松理解,不过尽管如此,我依旧建议你结合着我的另一个课程:《R 数据科学》一起学习。

视频讲解

本系列的课程的参考资料较多,我还会结合大量实际案例,因此此处我暂时不一一列举了。我会在每次课程的最后添加参考资料的链接。

为什么要使用 R 语言进行地理分析?

实际上我并没有学习过任何桌面 GIS 软件,例如 QGIS、ArcMap、GRASS、SAGA等,但是编程语言相对桌面软件的最大优势就是可重复性了。

所以我也只能从 R 的角度去看这个问题,R 的易用性、跨平台性以及活跃的社区都是这个问题的答案。基于我自己的使用经验,使用 R 进行地理计算可以满足我们 99% 的工作和学习需要。另外结合 R 语言强大的绘图性能,地理数据可视化也就不是问题了。

例如下图展示了夜间灯光数据和中国三个城市的位置:

1
2
3
4
5
6
7
library(leaflet)
popup = c("北京市", "拉萨市", "香港特别行政区")
leaflet(width = "100%", height = "400px") %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%
addMarkers(lng = c(116.41228, 91.089778, 114.13432),
lat = c(40.185543, 30.036993, 22.381274),
popup = popup)

可重复性:别人可以使用你的代码和数据得到相同的结果。

R 语言中有诸多用于地理计算和可视化的包,在这里我主要使用 sf、raster 和 ggplot2 包,另外我会尽可能介绍更多好用的包。我就不再赘述使用 R 语言进行地理分析的历史了,反正用 sf 包就对了。在数据量不大的时候这个包的性能表现还不错。

使用 R 语言进行地理分析能做什么?

这里展示下我之前的教程中涉及到的一些。

R 中的地理数据

准备工作

本节课主要使用下面两个 R 包:

1
2
install.packages('sf')
install.packages('raster')

地理矢量数据

地理矢量数据模型的基础是位于地理参考系统(crs)中的点。点可以表示为独立的要素(例如某个公交站台的位置),也可以组合成复杂的几何形状,例如直线和多边形。

例如 crs = 4326:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(sf)
library(ggplot2)
library(tidyverse)

worldmp <- read_sf("world.geo.json") %>%
st_transform(4326)

tibble(lon = c(0, 116),
lat = c(0, 40),
name = c("crs:4326 原点", "北京")) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) -> da4326

ggplot(da4326) +
geom_sf_label(aes(label = name), family = cnfont,
nudge_x = -12, nudge_y = -12,
color = "red") +
geom_sf(color = "red", size = 3) +
geom_sf(data = worldmp, fill = NA) +
coord_sf(crs = 4326) +
labs(title = "CRS: 4326") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

我觉得费力的去解释 4326 意义不大,因为它表示 WGS 84 坐标系,显然又变成另一个费解的词语,因为后面我都直接说 crs 的数值了。我国官方要求使用的坐标系是 CCGS 2000,对应的 crs 是:4490,这一坐标系的坐标原点和北京的位置:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
worldmp <- read_sf("world.geo.json") %>% 
st_transform(4490)

tibble(lon = c(0, 116),
lat = c(0, 40),
name = c("crs:4490 原点", "北京")) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(4490) -> da4490

da4490
da4326
ggplot(da4490) +
geom_sf_label(aes(label = name), family = cnfont,
nudge_x = -12, nudge_y = -12,
color = "red") +
geom_sf(color = "red", size = 3) +
geom_sf(data = worldmp, fill = NA) +
coord_sf(crs = 4490) +
labs(title = "CRS: 4490") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

是不是发现两个坐标系为什么差异,所以通常我们不需要特意将 4326 下的坐标转到 4490 下。两个坐标系下的坐标仅有细微的差异。

有很多 crs 可供选择:

1
2
3
# install.packages('rgdal')
rgdal::make_EPSG() %>%
DT::datatable(width = "100%")

sf 包也提供了类似的函数:

1
2
sf_proj_info(type = "proj") %>% 
DT::datatable(width = "100%")

中间的部分去哪了?

完整的讲义材料可以在加入我的知识星球(线上培训班)后从知识星球下载附件获取~
要了解如何加入,可以阅读关于界面或者添加我的微信咨询。

sf 对象的一个优势就是对单位的支持,例如我们计算世界各国的面积:

1
st_area(worldmp)

我们可以把它合并到 worldmp 数据框里面:

1
2
worldmp$area <- set_units(st_area(worldmp), km2)
worldmp

作业:

  1. 绘制一幅世界地图展示各国/地区的面积大小;
  2. 从世界地图数据里面筛选出中国的并绘图展示;
  3. 尝试探索更多的 crs 下世界地图的形状;
  4. 如果时间允许的话,可以试着从全球的夜间灯光数据里提取美国的部分并绘图展示。

参考资料

知识星球链接:https://t.zsxq.com/ZnmIyvB

#

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×