R可视乎｜空间地理数据可视化（1）
1. 前言
研究生讨论班第一次用 slides 作报告，主要讲了《Geospatial Health Data》[1]一书中关于空间地理数据可视化的内容。文末给出对应的 pdf 网页版本。
本篇主要介绍：用 R 包制作地图的基础内容，之后会再详细介绍数据可视化主要的 R 包和函数，敬请期待。由于本文内容较多，所以做了下思维导图：
2. 空间数据类型
一个 d 维的空间过程如下所示：
其中，s 为观测的位置，Z 指的是我们所观测的属性，例如，婴儿猝死数量或降雨量。通过域 D 的特征可将空间数据分为：区域数据（areal data）、地理统计数据（geostatistical data）、点模式数据（point patterns）。接下来分别对这三类数据进行介绍。
2.1 区域数据
区域数据中，域 D 是固定的并且被划分为具有明确边界的有限数量单元，人们常通过邮区编号、人口普查、像素报告的遥感数据等来收集获取区域数据。
例子：1974 年美国北卡罗来纳州各县婴儿猝死数量（Sudden Infant Death），如下图所示：
1974年美国北卡罗来纳州各县婴儿猝死数量
基于此，通过利用人口和其他协变量数据，可获得每个县的死亡风险估计。
2.2 地理统计数据
对于这种类型的数据，域 D 是一个连续的固定集合。连续是指 s 可以在 D 中连续地变化，Z(s)可以在 D 的任何地方被观测到，Z(s) 可以是连续的也可以是离散的；固定是指域 D 中的点是非随机的（non-stochastic）。
例子：下图是巴西巴拉那州 143 个记录站的旱季平均降雨量：
巴西巴拉那州143个记录站测量的平均降雨量
上图中的数据代表了特定站点获得的降雨量的测量值，可使用地理统计学的相关模型预测来未取样站点的降雨量。
2.3 点模式数据
与前两种数据不同，点模式数据中域 D 是随机的，s 给出了随机事件的位置。对于
，Z(s)表示事件的发生，其值可以为 1，也可以是随机地给出一些额外的信息。
例子：居住在城市中患有特定疾病的个人的地理坐标，如下图所示：
约翰·斯诺绘制的1854年伦敦霍乱爆发的地图
我们可以使用点过程方法分析这些数据来了解死亡的空间分布，并评估某位置附近是否存在过度的疾病风险。
3. 坐标参考系统
坐标参考系统（Coordinate Reference Systems，CRS），是用来表示空间数据的重要工具，通过使用坐标参考系统我们可以知道坐标的原点和测量单位。CRS主要有地理坐标参考系统（又称非投影坐标参考系统）和投影坐标参考系统两类。
地球的三维表面（左）和地球的二维表面（右）
3.1 地理坐标参考系统
- 使用经度和纬度来确定地球三维椭圆体表面上的位置。
- 纬度和经度是以十进制度（DD）或度、分、秒（DMS）为单位的角度。
- 地球表面一个点的纬度是赤道平面与通过该点和地球中心的直线之间的角度。
- 地球表面某一点的经度是指本初子午线以西或以东到一条经过该点的经线的角度。
地球的纬线（左）和经线（右）
3.2 投影坐标参考系统
投影是将地球的三维表面转化为某一个二维平面，所有的地图投影都会以某种方式扭曲地球表面，并不能同时保留所有的面积、方向、形状和距离属性。
最常用的投影方式是墨卡托投影（Universal Transverse Mercator，UTM），这种投影方式将地球划分为60个经度为6度的区域，每个区域都使用横向墨卡托投影，绘制出一个南北方向的范围。
例子：下图是CMG Lee 绘制的等距矩形世界地图的通用横轴墨卡托区域，其中不规则区域和纽约市突出显示：
CMG Lee 绘制的等距矩形世界地图上的通用横轴墨卡托区域
地球上的某一位置可由UTM区号、半球(南或北)以及区的东经和北纬的坐标(以米为单位)给出。关于这种投影进一步的细节，可查看维基百科[2]。
3.3 在 R 中设置坐标参考系统
地球的形状可以用一个扁椭球形的模型来近似，它在赤道上隆起，在两极扁平，目前世界上有很多不同的参考椭球体来使用，最常用的是全球定位系统（GPS）所使用的世界大地测量系统（WGS84）。
除此之外，还有欧洲石油调查组（EPSG）所制定的地图，由于坐标系的不同，各地的地图也会不同，例如中国：以地球几何球心为中心时，EPSG 代码为 4479；以地球椭球焦点为中心时，EPSG 代码为 4480。WGS84 的 EPSG 代码为 4326。
在 R 语言中，CRS 是用 proj4 字符串指定的，这些字符串指定了投影、椭球体和基准点的属性。例如，WGS84 经度/纬度投影被指定为
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
UTM 29 区的 proj4 字符串由以下公式给出
"+proj=utm +zone=29 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
而南方的 UTM29 区为
"+proj=utm +zone=29 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +south"
此外，如果我们希望将数据d转换为具有不同投影的数据，则可以使用 rgdal 包中的 spTransform() 函数或 sf 包中的 st_transform() 函数。
例子：创建一个由经度和纬度给出坐标的空间数据集，并使用 rgdal 将其转换为南方 UTM 35 区的坐标数据集：
library(rgdal)
# create data with coordinates given by longitude and latitude
d <- data.frame(long = rnorm(100, 0, 1), lat = rnorm(100, 0, 1))
coordinates(d) <- c("long", "lat")
# assign CRS WGS84 longitude/latitude
proj4string(d) <- CRS("+proj=longlat +ellps=WGS84
+datum=WGS84 +no_defs")
# reproject data from longitude/latitude to UTM zone 35 south
d_new <- spTransform(d, CRS("+proj=utm +zone=35 +ellps=WGS84
+datum=WGS84 +units=m +no_defs +south"))
# add columns UTMx and UTMy
d_new$UTMx <- coordinates(d_new)[, 1]
d_new$UTMy <- coordinates(d_new)[, 2]
4. 图形文件
- 空间地理数据是用一个叫做 shapefile 的数据存储格式来表示的，它可以储存位置、形状及其他地理属性。
- 一个 shapefile 是由一系列相关的文件组成，这些文件有不同的拓展名，并存储在同一个目录中。一个shapefile必须包括的三个文件为：
.shp、.shx和.dbf，可以构成 shapefile 的其他文件另有.prj、.sbn、.sbx和.shp.xml。 - 我们可以使用 rgdal 包中的
readOGR()函数，或者 sf 包中的st_read()函数来读取 shapefile 文件。
例子：用 readOGR() 读取存储在 sf 包中的北卡罗来纳州的 shapefile，如下所示：
# name of the shapefile of North Carolina of the sf package
nameshp <- system.file("shape/nc.shp", package = "sf")
library(rgdal)
map <- readOGR(nameshp, verbose = FALSE)
class(map)
head(map@data)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS
## 0 0.114 1.442 1825 1825 Ashe 37009
## 1 0.061 1.231 1827 1827 Alleghany 37005
## 2 0.143 1.630 1828 1828 Surry 37171
## 3 0.070 2.968 1831 1831 Currituck 37053
## 4 0.153 2.206 1832 1832 Northampton 37131
## 5 0.097 1.670 1833 1833 Hertford 37091
## FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79
## 0 37009 5 1091 1 10 1364 0
## 1 37005 3 487 0 10 542 3
## 2 37171 86 3188 5 208 3616 6
## 3 37053 27 508 1 123 830 2
## 4 37131 66 1421 9 1066 1606 3
## 5 37091 46 1452 7 954 1838 5
## NWBIR79
## 0 19
## 1 12
## 2 260
## 3 145
## 4 1197
## 5 1237
用 rgdal 包导入的北卡罗来纳州的地图如下图所示：
plot(map)
由 rgdal 包得到的美国北卡罗来纳州地图
- 用
st_read()读取地图：
# read shapefile with st_read()
library(sf)
map <- st_read(nameshp, quiet = TRUE)
class(map)
head(map)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS
## 0 0.114 1.442 1825 1825 Ashe 37009
## 1 0.061 1.231 1827 1827 Alleghany 37005
## 2 0.143 1.630 1828 1828 Surry 37171
## 3 0.070 2.968 1831 1831 Currituck 37053
## 4 0.153 2.206 1832 1832 Northampton 37131
## 5 0.097 1.670 1833 1833 Hertford 37091
## FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79
## 0 37009 5 1091 1 10 1364 0
## 1 37005 3 487 0 10 542 3
## 2 37171 86 3188 5 208 3616 6
## 3 37053 27 508 1 123 830 2
## 4 37131 66 1421 9 1066 1606 3
## 5 37091 46 1452 7 954 1838 5
## NWBIR79
## 0 19
## 1 12
## 2 260
## 3 145
## 4 1197
## 5 1237
用 sf 包导入的北卡罗来纳州的地图可以产生如下结果：
plot(map)
由 sf 包得到的美国北卡罗来纳州地图
小编有话说
- 本篇主要介绍：用 R 包制作地图的基础内容，包括：几种空间数据类型、不同的坐标参考系统介绍以及如何使用 R 包导入图形文件以及绘图。
- 这是空间地理数据可视化系列的第一期，主要由 林华师 制作。本系列的宗旨是带你系统学习如何使用 R 对空间地理数据进行可视化。
- 未来几期会具体介绍各类绘制空间地理数据的 R 包，敬请期待。
参考资料
[1]
《Geospatial Health Data》: https://www.paulamoraga.com/book-geospatial/sec-spatialdataandCRS.html
[2]
维基百科: www.wikipedia.org