本文引用的所有信息均为公开信息，仅代表作者本人观点，与就职单位无关。
本文首先以实际示例介绍三类典型的空间数据，希望读者能对不同类型的空间数据的特点和使用场景有个简单印象。紧接着，介绍各类地理空间数据的可视化方法及 R 包实现，笔者从个人经验出发详细介绍了 4 个，相信足够让读者应付大部分使用场景。之后，再介绍 R 语言中的空间数据类型及其基本操作，空间数据操作是基本功，也是比较难的，没放在空间数据可视化前面是怕吓跑了一些读者，笔者所知也十分有限，也怕写不好。再基于一个真实数据用不同的 R 包实现类似的可视化效果，是几代空间数据可视化方法的并列介绍，最后是一些心得、总结和展望。
本文有几个写作意图：
空间数据基本类型有点、线、面（或称多边形）等，还分矢量数据和栅格数据，地理可视化的展示形式有气泡图、热力图、瓦片图、地形图等。下面通过几个实际的空间数据，以图形的展示方式以期帮助读者建立数据直觉，陆续使用 lattice [1]、ggplot2 [2] 和 leaflet [3] 来绘制静态和交互图形。
下面以 RgoogleMaps 包内置的数据集 incidents 为例介绍空间点数据，并作一些简单的分析。
图 2.1: 2012 年旧金山警方出警地点的空间分布
除了 RgoogleMaps 提供的 plotmap() 函数，当然，还可以用 leaflet 来绘制交互式的散点图，因为交互，可以放入更多的信息。如图 2.2，鼠标悬浮在散点上可以看到事件发生的警区 PdDistrict，点击可以看到警方的处置方式 Resolution，这是设置 addCircles() 的参数 label 和 popup 实现的，实际上，它们还可以放入更加丰富的信息。
图 2.2: 2012 年旧金山警方出警地点的空间分布
接下来再细致一点地探索 incidents 数据集，各个警区出警次数的分布情况，
图 2.3: 事件发生次数随警区的分布
进一步，可将图 2.3 所示分布绘制在地图上，如图 2.4 展示事件发生次数随警区的空间分布。
图 2.4: 事件发生次数随警区的空间分布：左图非暴力事件，右图暴力事件
图 2.5: 事件发生次数随时段的分布
每次接报出警，警方会根据具体情况快速对事件定级分类，然后派出相应的警力以及处置方式。可以见事件性质 Category、是否发生暴力 violent 和警方处置方式 Resolution 应当有很强的关联。既然 Category， violent 和 Resolution 有很强的关联，展示关联关系非常重要，有的分类 Category 就应该给予足够的处置力度，特别对于暴力事件，关注警方的措施是否有不当或处置不力的？以此，还可以推测旧金山的治安状况。
图 2.6: 事件性质和警方处置情况
从图 2.6 看出，对于不同的事件 Category 有不同的执行措施 Resolution，执行措施和是否属于暴力有关系。一般来讲，暴力事件都应该有处置，发生暴力事件只有 5 类，按发生次数依次是打架斗殴、抢劫、性骚扰、绑架、盗窃，从图中也可看出盗窃是最多的犯罪事件。下表 2.1 是对图中分类的中文翻译，笔者英文水平一般，翻译仅供参考。
ggplot2 提供的 map_data() 函数可以从 maps 包中提取和转化地图数据，而 coord_map() 函数结合 mapproj 包可以实现坐标转化。下面以 R 内置的数据集 state.x77 为例，它记录了美国 1977 年 50 个州的人口普查情况，Population 各州的人口数量 （1975 年 7 月 1 日估计），Income 人均收入（1974 年），Illiteracy 文盲比例（人口百分比），平均寿命（1969-1971 年统计），Murder 谋杀和非过失杀人率（单位：十万分之一，1976 年统计），HS Grad 高中毕业生百分比（1970 年统计），Frost 首都或大城市最低气温低于冰点的平均天数（1931-1960 年统计），Area 土地面积（单位：平方英里）。如图 2.7， maps 包内置的美国各州的地图数据缺少夏威夷和阿拉斯加州。
图 2.7: 1975 年美国各个州的人口
图 2.8: 1975 年美国各个州的人口
图 2.9: terra 包绘制锡安国家公园地形
地理可视化的展示形式有多种，散点图、气泡图、峰窝图和热度图等，还可以有更加复合的专题地图，如在地图上省、市、区、县一级别，显示各个年龄段的人口比例，类似将饼图贴在地图上。为简单起见，主要以散点图为例介绍入门级别的内容。
图 3.1: echarts4r 使用内置矢量地图数据
图 3.2: echarts4r 使用内置矢量地图数据
下面介绍在瓦片地图数据的基础上展示 quakes 数据，如图 3.3 所示，这里采用地理图层 e_leaflet() 调用开放的 OpenStreetMap 瓦片服务。
图 3.3: echarts4r 使用 OpenStreetMap 地图
只要提供瓦片地图服务的模版，就可以切换地图底图数据，图 3.4 使用了谷歌的地图服务。
图 3.4: echarts4r 使用谷歌地图
类似谷歌的要求，使用 Mapbox 地图服务也需要先申请个人访问令牌，申请下来后，也把它设置到 R 语言环境变量里 MAPBOX_TOKEN，以便后续使用。
图 3.5: echarts4r 使用 MAPBOX 地图
在国内，最好使用百度、高德这样的国家认可的地图服务提供商，这也很容易，只需用高德的 tile 服务数据替换 OSM 的。
图 3.6: echarts4r 使用高德地图
图 3.7: QGIS 软件内的 OpenStreetMap 地图
图 3.8: 地震震级映射连续的颜色值
如图 3.9 所示，用不同颜色表示不同种类的鸢尾花。
图 3.9: 三种鸢尾花用不同颜色标记
实际上，还可以将地震震级分段，用梯度变化的颜色表示震级的大小，如图 3.10 所示。
图 3.10: 散点图
绘制热力图就这么简单，如图 3.11 所示，在图 3.10 的基础上添加了热力图，实际上就是根据点的密度添加色彩，密度值根据核密度估计算法统计出来。
图 3.11: 热力图
图 3.12: plotly 包使用 Mapbox 地图
图 3.13: plotly 包使用 Mapbox 卫星地图
图 3.14: plotly 包使用 Mapbox 绘制热力图
plotly 包内置了两个地图数据集，一个是美国州级多边形数据"USA-states"，另一个是国家级多边形数据，来自 Natural Earth，需要 3 个字母的国家代码传给参数 locations。
图 3.15: plotly 包使用 Mapbox 绘制面量图
多边形矢量地图数据，多边形的边界常常是行政区域，choropleth 地图是填充颜色的多边形，颜色深浅表示数据指标，地图的作用在于展示指标的空间变化，在美国，郡县一级的行政区大致相当于咱们市一级行政单位，比如图 3.16 展示美国各个县的失业率，大城市失业率相对较高。值得注意，数据量较大的时 plotly 渲染速度较慢。
图 3.16: plotly 包使用 Mapbox 绘制面量图
2001 年 lattice 包就发布到 CRAN 了，已经过去 20 多年了，功能相当全面成熟稳定，内置于 R 软件。下面介绍 lattice 包的数据可视化能力，继续以 quakes 数据集为例展开介绍。
图 3.17: 散点图
图 3.18: 热力图
图 3.19: 气泡图
图 3.20: 蜂窝图
图 3.21: 等高图
图 3.22: lattice 包绘制三维气泡图
图 3.23: rgl 包绘制三维气泡图
结合三维气泡图，也许读者已经发现震级和震深有一些关系，不妨再看看仔细它们之间关系，如图 3.24。
图 3.24: 抖动图
无论震级如何，地震似乎集中发生在表层或深层的位置，中间带反而比较稀疏，可见震级和地震的位置关系不大，地不地震、震级如何主要是由地质条件决定的。
本节开始介绍一些空间数据相关的基础知识，R 语言有一些扩展包是专门处理空间数据的，提供基础数据结构和类型，了解这些是为空间数据操作、分析、建模打基础。空间数据操作（Spatial Data Manipulation），如计算空间区域的面积、空间点的路径、区域叠加、基于空间关系的数据查询等。空间数据分析（Spatial Data Analysis），空间数据的描述性分析和探索性分析，帮助发现数据中的模式规律。空间统计分析（Spatial Statistical Analysis），传统的统计模型认为观测值之间常常是独立的，而空间数据自带空间相关性，需要专门的统计分析方法。空间数据建模（Spatial Data Modeling），面对不同的空间数据类型需要不同的建模方法，比如空间点模式分析（Spatial Point Pattern Analysis）、格网 / 面状数据分析（Area Data/Grid Data/lattice Data Analysis）、地统计分析（Geostatistics）等。
sp [7] 是一个 R 包，提供一种存储和操作空间数据的对象类型，具体地，提供点 SpatialPoints / SpatialMultiPoints、线 SpatialLines、多边形 SpatialPolygons、栅格 SpatialGrid / SpatialPixels 等基础空间数据类型以及相应构造方法，还有聚合 aggregate()，合并 merge()，转化 spTransform() 和绘图 spplot() / bubble() / image() 等空间数据操作和可视化函数。
可以看出 sp 打造的一套空间数据基础类型及其衍生类型，知道数据类型才知道后续该怎么对其进行操作，调用合适的方法。比如最常使用的绘图函数 plot() 已经汇集很多方法了，它也是一个绘图对象，对于不同类型的数据对象，它会调用相应的绘图方法。在 R 语言内，一切都是对象，一切操作都是函数调用。
sp 提供最基础的空间数据类是 Spatial 类，下面是一个实例。
空间对象 S 有两个重要的部分，地理边界和坐标参考系，这两个信息可以分别用函数 bbox() 和 slot(, "proj4string") 单独提取。
实际上，除了函数 bbox() 和 slot(, "proj4string")，Spatial 类还有很多方法。
查看转化后的数据情况：
直接调用 plot() 方法，获取图 4.1。实际上，会调用 SpatialPoints 类的 plot() 方法，plot() 本身也是一个泛型函数，R 语言中的泛型有点类似面向对象编程中的多态。
图 4.1: 空间数据绘图
调 sp 数据对象特有的方法 spplot()，获取图 4.2。
图 4.2: 空间数据绘图
地震主要是由于地壳运动产生的，两条地震带由内圈到外圈，地震深度有明显的阶梯变化，暗示着断层的方向和位置。
数据集 alaska 只有 1 行 6 列，使用的坐标参考系是 Alaska Albers (EPSG:3467)。各个变量分别是 GEOID 地理标识符（geographic identifiers），NAME 州的名称，REGION 区域的名称，AREA 面积（平方千米），total_pop_10 2010 年人口，total_pop_15 2015 年人口，geometry 表示 sf 多边形类 MULTIPOLYGON 的集合属性，阿拉斯加地区是由多个多边形组成的，通过图 4.3 可以直观地看出来。
图 4.3: sf 包绘制多边形数据
图 4.4: 2014 年世界各国人均 GDP
spData 包内置的 world 数据集将台湾和中国大陆的地区用不同的颜色标记，严格来讲，这是不符合国家地图规范的，不能出现在期刊、书籍等正式的出版物里。
先简要介绍 Base R 内置的一般数据操作，也是最常见的，然后介绍 sf 包和空间数据操作。sf 包对 Base R 内置的一些函数做了扩展以支持空间数据的合并、聚合等。
通过排列组合 R 软件内置的一些基础数据操作函数能够满足大部分使用场景，这些基础函数都来自 base 包和 stats 包，由 R 语言 核心团队维护，约 10 来个函数：
函数 subset() 提取数据框中的某些列和行，也可以用 $ 或 [ 实现。
函数 transform() 可以在原来的数据框上添加新的一列，这个新的列通常是根据已有的列加减乘除等计算而来。
函数 order() 按照数据框的某个（些）列排序，如先按照第一列从小到大排序，再按照第二列从小到大排序。
函数 aggregate() 分组聚合函数，也是很常用的。
函数 reshape() 将长格式的数据框变形为宽格式的数据框，反之亦可。
函数 merge() 合并两个数据框，有内关联、左关联、右关联等合并方式，合并操作可以和 SQL 中的关联操作对应。
函数 duplicated() 查出是否有重复的记录，后续可以决定是否去掉。
下面以 R 软件内置的数据集 mtcars 为例，介绍数据操作的一个典型过程，拆分、计算、合并，早年 Hadley Wickham 创建 plyr 包 [9] 就是干这个事情，后来，陆续出现很多 R 包，如 reshape 包（已退休）、reshape2（已退休）、dplyr（活跃）、purrr（活跃）、tidyr（活跃），一直处于不太稳定的状态，此处不再细说。
下面举个简单的例子，计算两点之间的距离，首先构造几何对象「点」，然后计算两点之间的距离。
默认采用的是笛卡尔坐标系，二维、三维、四维都可以，两个点的情况下，相当于计算平面上两个点的距离，平面上两个点之间直线最短，即欧式距离：
就这个例子而言，st_distance() 函数的作用相当于：
sf 包提供的函数 st_crs() 可以非常方便地检索坐标系的信息，以 WGS 84 为例。
可知 WGS 84 坐标系采用的基准面是椭球面，椭球的赤道半径为 6378137 米，扁率为 1/298.257223563。
空间数据对象常常围绕点、线、多边形或多维数组（用来表示卫星图像、时空数据）等展开，可以把这些几何对象看作是主键一样的东西，在同一坐标参考系下，一个地理位置或一块区域是不会有两个相同的地理标记的，其它观测变量是维度或指标，也可以看作这块区域的属性，比如说邵东市地区（一个多边形区域），地区拥有的人口数量、土地面积、人均 GDP 等。
平面几何，两点之间直线最短，距离就是欧氏距离。
球面几何，两点之间最短距离是大圆的弧长，点的坐标视为经纬度参与计算。
可见，R 包 sf 和 Presto 计算的结果完全一致，Presto 采用的空间坐标参考系也是 WGS 84。
nc2 是一种宽格式数据，有两列 SID74 和 SID79，现在希望在地图上以 ggplot2 的分面绘图方式同时展示两个变量的空间分布，因此，需要先将数据由「宽格式」转为「长格式」。转化前，先来看下 nc2 的数据类型：
图 6.1: sf 数据可视化
图 6.2: leaflet 数据可视化
实现类似功能的还有 mapdeck 包。
在 sf 包出来前，maptools 曾经是空间数据操作的主力，是读取 *.shp 文件的主要方式。配合 sp 包，是实现空间数据分析和可视化的主要帮手。
sp 包提供的函数 spplot() 方法绘制类似的分面图，这其中借助了 lattice 包的绘图能力，效果如 6.3 所示。
图 6.3: lattice 数据可视化
ggmap 包获取的 Google 瓦片地图和 leaflet 包 使用的 OpenStreetMap 地图在坐标参考系上是一样的，都是 EPSG:3857，想把数据绘制到地图上，需要先将数据的坐标参考系转换到 EPSG:3857，坐标单位从经纬度随之变为米。
图 6.4: sp 数据可视化
latticeExtra [11] 在 lattice 的基础上，封装了一些高级函数，围绕本篇主题，下面介绍一下 mapplot() 函数，读者不妨类比 spplot() 函数，可知它是专门为 map 数据类型提供绘图服务。 map 数据类型在 R 语言中存在很多年了，算得上是古老了，很久以前绘制地图的 R 包主要有：提供基础地图数据的 maps 包， 提供更多常用地图数据的 mapdata 包， 实现坐标投影转化的 mapproj 包 和数据读写、操作的 maptools 包。
参数 x 提供符合 R 语言语法的公式，比如 y ~ x1 + x2。参数 data 提供数据集，包含观测变量。参数 map 提供地理边界信息，地理单元的名称必须和公式中变量 y 的值匹配。
图 6.5: latticeExtra 数据可视化
本文意在提供一个基于 R 语言的空间数据结构、可视化的简单概览，帮助读者了解空间数据可视化的轮廓，提供 R 语言在此领域的探索和实践。
从 sp 到 sf，从 raster 到 terra，是 R 语言在空间数据领域前进的一大步，在使用的直观感受上，数据读写、操作的性能得到极大的提升，整个空间数据处理领域的 R 包依赖也将大为减少。
ggplot2 及其生态是很好的，可能大家不知道的是，lattice 在绘图能力上丝毫不逊色于 ggplot2，甚至还有超越。本篇主要介绍空间数据，下面以 R 软件内置的地形数据 volcano 为例。
本文只涉及空间数据领域的皮毛，说白了就是导入几种格式的空间数据画点小地图，也未涉及百万乃至百亿级的大规模空间数据的可视化，以及支持可视分析需要的空间数据存储、检索、操作、分析、建模和应用的知识和技术。此外，对于不同类型的空间数据有不同的分析方法和统计模型，比如点数据的模式分析，面数据的层次分析，栅格数据的图像分析等，本文亦未涉及。再有，面向具体业务问题，从问题发现到数据模型的转化、界定、拆解、优化和应用所需要的领域知识，本文也未涉及。笔者所知所学有限，空间数据领域是如此广阔，需要大家一起发力。
各种花里胡哨引入一堆 R 包和函数，直接告诉你答案，又不告诉你为什么，三天不学习就跟不少变化。↩︎
基本动作就是内功心法，无论面对的数据操作如何复杂，最终都能一一化解，以不变应万变。↩︎