数据地图是一种经典的图示方法，在R软件中，各种程序包所提供的函数在绘制数据地图时比较方便。

基本地图的绘制

世界地图

方法一
R包中存储着常见地图的数据，比如maps包中包含世界地图、美国地图等

library(maps)
map("world", fill = TRUE, col = rainbow(200),
    ylim = c(-60, 90), mar = c(0, 0, 0, 0))
title("世界地图")

方法二

library(maps)
data("world.cities")
bigcities <- subset(world.cities, pop > 5000000)
qplot(long, lat, data = bigcities,colour=country.etc,size=pop)+ borders("world", size= 0.5)

中国地图

利用ggplot2绘制美国地图比较容易，因为美国地图是程序包maps中自带的地图，但其他国家的地图数据则需要从外部导入。

方法一
首先，安装相应的包

install.packages("mapproj")
install.packages(“sp”)
install.packages(“maptools”)
library(maptools)
library(ggplot2)
library(plyr)
library(sp)

下载中国地图的GIS数据
这是一个压缩包，完全解压后一般包含三个文件（bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx），将这三个文件解压到R的工作空间下，或者直接输入路径：

x <- readShapePoly("C:/Users/SAMSUNG/Desktop/china-province-border-data/bou2_4p.shp")
plot(x)

如果安装maptools后无法用readShapePoly（）函数读取.shp文件，可以试试rgdal程序包，然后用rgdal包中的readOGR（）函数读取.shp文件：

library(rgdal)
china_map <- readOGR("C:/Users/SAMSUNG/Desktop/china-province-border-data/bou2_4p.shp")
china_map1 <- fortify(china_map)

#绘制全国地图-不着色
ggplot(china_map1,aes(x=long,y=lat,group=group))+
   geom_polygon(fill="white",colour="black")+
   coord_map("polyconic")
   
#绘制全国地图-分省着色
xs<-data.frame(x,id=seq(0:924)-1)
china_map_data<-join(china_map1, xs, type = "full")  

ggplot(china_map_data, aes(x = long, y = lat))+
  geom_polygon(aes(group = group,fill=NAME),color="grey40" )+ #线条色
  coord_map("polyconic")+
  scale_fill_manual(values=colours(),guide=FALSE) #分省着色

方法二
从gadm.org网站上得到中国的省区地理数据，并加载到R软件内存中。GADM是世界行政区域（或行政区域界线）位置的空间数据库，可专门用于地理信息系统和类似的软件，并且提供Rdata格式的数据，非常便于在R中绘制数据地图。
但是。。。你会发现一些令人气愤的事，如下图【笔者已经无fuck说了

代码如下，可自己补全路径

install.packages(“sp”) 
library(sp)  
load(url("http://gadm.org/.../***.RData"))
#将每个省人口数据按顺序存放在数据框gadm中，生成一个变量pop
gadm$pop=c(1961,1293,7185,3571,2470,4374,2745,3831,2301,7866,5442,
           5950,3689,4456,9579,9402,5723,6570,10432,4602,867,2884,
           8041,3474,4596,300,3732,2557,562,630,2181,706)
#利用空间绘图命令进行绘图
spplot(gadm,"pop",col.regions = rev(terrain.colors(gadm$pop)),main="中国各省人口数据")

方法三

library(mapdata)
ch_cities <- subset(world.cities, country.etc=="China")
ggplot(ch_cities, aes(long, lat)) + geom_point(colour= alpha("red",0.5))+ borders("china")

方法四

map("china", col = "red4", ylim = c(18, 54), panel.first = grid())

实例一：NBA球队

NBA比赛球队30支,现想了解詹姆斯（LBJ）对阵哪个球队时表现更好，或得分更高，故通过绘制数据地图直观地得到结论。

library(ggplot2)
lbj=read.table("C:/Users/SAMSUNG/Desktop/R数据分析/data/lbj.txt",header=T,quote="'") 

attach(lbj)
head(lbj)  #查看数据集的前5行
state_map=map_data("state")  #获取美国地图的数据信息

p=ggplot(lbj,aes(map_id=state))+geom_map(aes(fill=AvgPTS),map=state_map)+
  expand_limits(x=state_map$long,y=state_map$lat)+
  scale_fill_continuous(limits=c(19,max(AvgPTS)),high='red3',low='yellow',guide="colorbar")+
  options(title='詹姆斯客场平均得分')

attach(state_map)
state.uni=unique(region)  #存放各州的名称
xx=0;yy=0  #事先建立变量xx和yy，下面用循环找到每个州对应的坐标值
for(i in 1:length(state.uni))
{
  xx[i]=mean(long[region==state.uni[i]])
  yy[i]=mean(lat[region==state.uni[i]])
}

order=0  #按变量state.uni的顺序找到数据集lbj中各州的位置，存放于变量order
for(i in 1:length(state.uni)){
  order[i]=which(state==state.uni[i])
}
labels=Opp[order]  #通过位置找到各州对应的球队名称
p+annotate("text",x=xx,y=yy,label=labels)  #最后绘图并添加注释

运行结果如下：

区域颜色越深，表示詹姆斯在该球队所得的分数越高。
用ggplot2画美国地图是相当的方便。

实例二：新冠肺炎

方法一
数据是截至2020年2月28日00：00，各省份的累计确诊病例数

library(sp)
library(maptools);
x=readShapePoly('.../bou2_4p.shp');
getColor=function(mapdata,provname,provcol,othercol)
{
f=function(x,y) ifelse(x %in% y,which(y==x),0);
colIndex=sapply(mapdata$NAME,f,provname);
fg=c(othercol,provcol)[colIndex+1];
return(fg);
#mapdata 存放地图数据的变量x
#provname 需要改变颜色的地区的名称
#provcol 对应于provname的代表颜色的向量（名称或数字）
#othercol 是其它地区的颜色
}
provname=c("北京市","天津市","河北省","山西省","内蒙古自治区",
"辽宁省","吉林省","黑龙江省","上海市","江苏省",
"浙江省","安徽省","福建省","江西省","山东省",
"河南省","湖北省","湖南省","广东省",
"广西壮族自治区","海南省","重庆市","四川省","贵州省",
"云南省","西藏自治区","陕西省","甘肃省","青海省",
"宁夏回族自治区","新疆维吾尔自治区","台湾省",
"香港特别行政区");
pop=c(410,136,318,133,75,121,93,480,337,631,1205,990,296,935,756,1272,65914,1017,1348,252,168,
576,538,146,174,1,245,91,18,72,76,32,93);
provcol=grey((1-pop/max(pop)))
plot(x,col=getColor(x,provname,provcol,"white"),xlab="",ylab="");

颜色数据部分运行如下：

方法二
将csv文件中的省份数据和程序包中所给的地图数据作矩阵合并

data<-read.csv("C:/Users/SAMSUNG/Desktop/data.csv",T) #读取人口数据
names(data)
china_data <- join(china_map_data, data, type="full") #矩阵合并
names(china_data)

ggplot(data,aes(x=long,y=lat,id=id,fill=number))+ #以人口数填充色彩
  geom_polygon(colour="grey40")+ #省界颜色
  scale_fill_gradient(low="white",high="steelblue")+#指定渐变填充色
  coord_map("polyconic")+#指定投影方式为polyconic
  theme(               #清除不需要的元素
    panel.grid = element_blank(),
    panel.background = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    legend.position = c(0.2,0.3)
  )

可能是因为GADM提供的id不同，矩阵合并时出现了问题，这个地方暂时还没解决，先留着。

空间数据可视化&地图绘制&R语言可复现

绘制地理空间数据是一项常见的可视化任务，需要专门的工具。通常，问题可以分解为两个问题：

1. 使用一个数据源绘制地图
2. 将来自另一个信息源的元数据添加到地图中。

有许多可以用于绘制地图的工具，一些流行的 GIS 软件允许点击交互进行地图开发和空间分析。这些工具具有不需要学习代码以及在地图上手动选择和放置图标和功能的便利性等优点，例如：

ArcGIS - 由 ESRI 公司开发的商业 GIS 软件，非常流行但相当昂贵

QGIS - 一个免费的开源 GIS 软件，几乎可以做 ArcGIS 可以做的任何事情

使用 R 作为 GIS 一开始似乎更令人生畏，因为它不是“点击”，而是有一个“命令行界面”（您必须编写代码才能获得所需的结果）。但是，如果您需要重复生成地图或创建可重现的分析，这是一个主要优势。

本文主要解释使用R语言绘制地图，主要用到的包有:

• rio : 导入数据
• tidyverse : 清洗、处理、绘图（包含ggplot2包）
• sf : 使用简单要素格式管理空间数据
• tmap : 生成地图，适用于交互式和静态地图
• OpenStreetMap : 在ggplot中免费添加 OSM 国内外地图作为底图

本目文录：

1. 综述

• 空间数据类型 Spatial data
  • 矢量数据 Vector data
  • 栅格数据 Raster data
• 空间数据基本格式
• 基本地图类型

2. 绘制地图准备工作

• 必备环节
• 所需R包
• 实例案例数据
• 行政边界数据
• 人口数据
• 卫生设施数据上

3. 画坐标系
4. 空间数据连接类型

• 点数据与多边形区域连接
• 最近的邻域连接
• 缓冲区连接
• 其他空间连接类型

5. 等值域图 Choropleth 地图
6. 使用ggplot绘制地图
7. 底图 Basemaps

• 免费获取国内外地图OpenStreetMap
• 等高线密度热图 Contoured density heatmap
• 时间序列热图 Time series heatmap

1.综述

您的数据的空间方面可以提供很多关于爆发情况的见解，并回答以下问题：

• 当前的疾病热点在哪里？
• 随着时间的推移，热点发生了怎样的变化？
• 卫生设施的使用情况如何？是否需要任何改进？

此 GIS 页面的当前重点是解决应用流行病学家在疫情应对中的需求。我们将探索使用 tmap 和 ggplot2 包的基本空间数据可视化方法。我们还将通过 sf 包介绍一些基本的空间数据管理和查询方法。最后，我们将使用 spdep 包简要介绍空间统计的概念，例如空间关系、空间自相关和空间回归。

1.1空间数据类型 Spatial data

**地理信息系统 (GIS) **- GIS 是用于收集、管理、分析和可视化空间数据的框架或环境，GIS 中使用的空间数据的两种主要形式是矢量Vector和栅格Raster数据：

矢量数据 Vector data

GIS 中使用的最常见的空间数据格式，矢量数据由顶点vertices和路径paths的几何特征组成。矢量空间数据可以进一步分为三种广泛使用的类型：

• 点 Points - 点由表示坐标系中特定位置的坐标对 (x,y) 组成。点是最基本的空间数据形式，可用于在地图上表示一个病例（即患者家）或位置（即医院）。
• 线 Lines - 一条线由两个相连的点组成。线有长度，可用于表示道路或河流等事物。
• 多边形 Polygons - 多边形由至少三个由点连接的线段组成。多边形特征具有长度（即区域的周长）以及面积测量值。多边形可用于标注区域（即村庄）或结构（即医院的实际区域）。

栅格数据 Raster data

一种替代格式 空间数据，栅格数据是一个单元格矩阵（例如像素），每个单元格包含高度、温度、坡度、森林覆盖等信息。这些通常是航空照片、卫星图像等。栅格也可以用作“底图base maps”放在矢量数据下方。

1.2空间数据基本格式

为了在地图上直观地表示空间数据，GIS 软件要求您提供足够的信息，说明不同要素应在何处相互关联。如果您使用的是矢量数据，这对于大多数用例都是正确的，则此信息通常会存储在 shapefile 中：

Shapefiles - shapefile 是一种常见的数据格式，用于存储由线、点或多边形组成的“矢量”空间数据.单个 shapefile 实际上是至少三个文件的集合 - .shp、.shx 和 .dbf。所有这些子组件文件必须存在于给定目录（文件夹）中，以便 shapefile 可读。这些相关文件可以压缩到 ZIP 文件夹中，以便通过电子邮件发送或从网站下载。 shapefile 将包含有关要素本身的信息，以及它们在地球表面的位置。这很重要，因为虽然地球是一个球体，但地图通常是二维的。关于如何的选择 “展平flatten”空间数据会对最终地图的外观和解释产生重大影响。

坐标参考系统 (CRS) - CRS 是一种基于坐标的系统，用于定位地球表面的地理特征。它有几个关键组件：

• 坐标系 Coordinate System - 有许多不同的坐标系，因此请确保您知道您的坐标来自哪个系统。纬度/经度的度数很常见，但您也可以看到 UTM 坐标。

• 单位 Units - 了解坐标系的单位（例如十进制度 decimal degrees、米）

• 基准 Datum - 地球的特定模型版本。多年来，这些已被修订，因此请确保您的地图图层使用相同的基准。

• 投影 Projection - 对用于将真正的圆形地球投影到平坦表面（地图）的数学方程式的参考。

请记住，您可以在不使用下面显示的映射工具的情况下汇总空间数据。有时只需要一张按地理（例如地区、国家等）的简单表格！

1.3用于空间数据可视化的基本地图类型

分区统计地图 Choropleth map - 一种专题地图，其中颜色、阴影或图案用于表示与其属性值相关的地理区域。例如，较大的值可以用比较小值更深的颜色来表示。这种类型的地图在可视化变量以及它如何在定义的区域或地缘政治区域之间变化时特别有用。

案例密度热图 Case density heatmap - 一种主题图，其中颜色用于表示值的强度，但是，它不使用定义的区域或地缘政治边界对数据进行分组。这种类型的地图通常用于显示“热点”或具有高密度或点集中的区域。

点密度图 Dot density map - 一种专题地图类型，使用点来表示数据中的属性值。这种类型的地图最适合用于可视化数据的分散情况并直观地扫描集群。

比例符号地图（分级符号地图）Proportional symbols map (graduated symbols map) - 类似于等值线地图的专题地图，但不是使用颜色来指示属性的值，而是使用与值相关的符号（通常是圆形）。例如，较大的值可以用比较小值更大的符号来表示。当您想要跨地理区域可视化数据的大小或数量时，最好使用这种类型的地图。

您还可以组合几种不同类型的可视化来显示复杂的地理模式。例如，下图中的病例（点）根据其最近的医疗机构（参见图例）进行着色。大红色圆圈表示一定半径的医疗机构集水区，而鲜红色的案例点表示在任何范围之外的区域：

2.绘制地图准备工作

2.1制作地图时的几个关键项目，其中包括：

• 数据集——可以是空间数据格式（例如 shapefile，如上所述），也可以不是空间格式（例如，只是 csv）。
• 如果您的数据集不是空间格式，您还需要一个参考数据集 reference dataset。参考数据由数据的空间表示和相关属性组成，其中包括包含特定要素的位置和地址信息的材料。
  • 如果您使用预定义的地理边界（例如，行政区域），通常可以从政府机构或数据共享组织免费下载参考 shapefile。如有疑问，最好从谷歌“[regions] shapefile”开始。
  • 如果您有地址信息，但没有纬度和经度，您可能需要使用地理编码引擎来获取空间参考数据以供记录。

关于如何将数据集中的信息呈现。有许多不同类型的地图，重要的是要考虑哪种类型的地图最适合您的需求。

2.2导入所需要的包

此代码块显示分析所需的包的加载。在本文，我们强调来自 pacman 的 p_load()，它会在必要时安装包并加载它以供使用。您还可以使用基本 R 中的 library() 加载已安装的包。

if(!require(pacman))install.packages("pacman")

pacman::p_load(
  rio,           # to import data
  here,          # to locate files
  tidyverse,     # to clean, handle, and plot the data (includes ggplot2 package)
  sf,            # to manage spatial data using a Simple Feature format
  tmap,          # to produce simple maps, works for both interactive and static maps
  janitor,       # to clean column names
  OpenStreetMap, # to add OSM basemap in ggplot map
  spdep,          # spatial statistics
  showtext        # 支持ggplot显示中文
  ) 

2.3示例案例数据

出于演示目的，我们将使用来自模拟的埃博拉流行病线列表数据框中的 1000 个病例的随机样本（在计算上，使用较少的病例更容易在本手册中显示）。linelist_cleaned.rds,数据集获取方式见文末。

由于我们对案例进行随机抽样，因此当您自行运行代码时，您的结果可能与此处演示的结果略有不同。

使用 rio 包中的 import() 函数导入数据（它处理许多文件类型，如 .xlsx、.csv、.rds ）。

# 使用import函数导入数据
linelist = import(here("data", "Epidemiologist_R","linelist_cleaned.rds"))

接着使用sample()函数随机抽取1000个样本

# 随机抽取1000个行标签
sample_rows = sample(nrow(linelist), 1000)
# 子集只包含样本行和所有的列
linelist = linelist[sample_rows,]

现在我们想将linelist(dataframe格式)转换为“sf”类（空间特征）的对象，其中linelist 有两列“lon”和“lat”代表每个案例居住地的经度和纬度。

我们使用包 sf（空间特征）及其函数 st_as_sf() 来创建我们称为 linelist_sf 的新对象。这个新对象看起来与 linelist 基本相同，但 lon 和 lat 列已指定为坐标列，并且已为显示点时分配了坐标参考系 (CRS)。 4326 表示我们的坐标标识为基于 1984 年世界大地测量系统 (WGS84) - 这是 GPS 坐标的标准。

# 创建一个空间数据集
linelist_sf <- linelist %>%
     sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

这就是原始 linelist 数据框的样子。在此演示中，我们将仅使用列 date_onset 和geometry（由上面的经度和纬度字段构成，是数据框中的最后一列）。

head(linelist_sf)

## Simple feature collection with 6 features and 28 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.26978 ymin: 8.453176 xmax: -13.21946 ymax: 8.484195
## Geodetic CRS:  WGS 84
##      case_id generation date_infection date_onset date_hospitalisation date_outcome outcome gender age age_unit age_years age_cat age_cat5                             hospital infector source wt_kg ht_cm ct_blood fever chills cough aches vomit temp time_admission      bmi days_onset_hosp                   geometry
## 5680  78703f         23     2014-12-10 2014-12-13           2014-12-16   2014-12-20   Death      m  16    years        16   15-19    15-19                              Missing   5bf94a  other    57   167       17   yes     no    no    no   yes 39.3          15:01 20.43817               3 POINT (-13.25821 8.453866)
## 3020  525a03         14     2014-10-21 2014-11-05           2014-11-07   2014-11-16   Death      f   5    years         5     5-9      5-9 St. Mark's Maternity Hospital (SMMH)   5bc9be  other    46    71       23   yes     no   yes    no   yes 38.9          11:25 91.25174               2 POINT (-13.22085 8.453176)
## 5154  19d49c         15     2014-09-06 2014-09-09           2014-09-16         <NA>   Death      f  10    years        10   10-14    10-14                        Port Hospital   bf6471  other    39   101       17   yes     no   yes    no    no 38.5          06:39 38.23155               7 POINT (-13.26978 8.480407)
## 3963  76a224         21     2015-01-05 2015-01-20           2015-01-21   2015-01-27 Recover   <NA>  15    years        15   15-19    15-19                    Military Hospital   674201  other    59   152       22   yes     no   yes    no    no 40.0          09:16 25.53670               1 POINT (-13.22583 8.484195)
## 2268  5355c3         16           <NA> 2014-10-04           2014-10-05   2014-10-17    <NA>      f  24    years        24   20-29    20-24                        Port Hospital     <NA>   <NA>    63   153       23   yes     no   yes   yes    no 39.7          12:52 26.91273               1 POINT (-13.21946 8.481999)
## 3539  185ac2         20     2014-11-19 2014-12-10           2014-12-11   2014-12-20   Death      m   5    years         5     5-9      5-9                    Military Hospital   4d1c91  other    48    93       22   yes     no   yes    no    no 38.8          13:25 55.49775               1 POINT (-13.25619 8.483108)

# 运行以下代码，得到交互式表格
# DT::datatable(head(linelist_sf, 10), rownames = FALSE, 
#               options = list(pageLength = 8, scrollX=T), class = 'white-space: nowrap' )

2.4塞拉利昂行政边界 Admin boundary shapefile

预先，我们已从或者，您可以通过我们的 The epidemiologist R下载所有示例数据或,者获取方式见文末。

现在我们将执行以下操作以在 R 中保存 Admin Level 3 shapefile：

1. 导入 shapefile
2. 清理列名
3. 过滤行以仅保留感兴趣的区域

要导入 shapefile，我们使用 sf 中的 read_sf() 函数。它通过 here() 提供文件路径。 - 在我们的例子中，该文件位于我们的 R 项目中的“data”、“gis”和“shp”子文件夹中，文件名为“sle_adm3.shp”（有关更多信息，请参见导入和导出和 R 项目页面）。您需要提供自己的文件路径。接下来我们使用 janitor 包中的 clean_names() 来标准化 shapefile 的列名。我们还使用 filter() 只保留 admin2name 为“Western Area Urban”或“Western Area Rural”的行。

# ADM3 level clean
sle_adm3_raw = st_read(here("data","Epidemiologist_R", "sle_adm3.shp"))

## Reading layer `sle_adm3' from data source `/Users/cpf/Documents/paper/writting_blog/data/Epidemiologist_R/sle_adm3.shp' using driver `ESRI Shapefile'
## Simple feature collection with 167 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.30901 ymin: 6.923379 xmax: -10.27056 ymax: 9.999253
## Geodetic CRS:  WGS 84

sle_adm3 <- sle_adm3_raw %>%
  clean_names() %>% # standardize column names
  filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # filter to keep certain areas

您可以在下面看到导入和清理后 shapefile 的外观。

head(sle_adm3)

## Simple feature collection with 6 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.29416 ymin: 8.094272 xmax: -12.91333 ymax: 8.494073
## Geodetic CRS:  WGS 84
##   objectid     admin3name admin3pcod    admin3ref_n         admin2name admin2pcod admin1name admin1pcod   admin0name admin0pcod       date   valid_on valid_to shape_leng   shape_area rowcacode0 rowcacode1 rowcacode2   rowcacode3                       geometry
## 1      155     Koya Rural   SL040101     Koya Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.63822264 0.0136601326        SLE     SLE004  SLE004001 SLE004001001 MULTIPOLYGON (((-13.02082 8...
## 2      156 Mountain Rural   SL040102 Mountain Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.29268365 0.0031823435        SLE     SLE004  SLE004001 SLE004001002 MULTIPOLYGON (((-13.21496 8...
## 3      157 Waterloo Rural   SL040103 Waterloo Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.72265524 0.0136413976        SLE     SLE004  SLE004001 SLE004001003 MULTIPOLYGON (((-13.12215 8...
## 4      158     York Rural   SL040104     York Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 1.23916989 0.0198344753        SLE     SLE004  SLE004001 SLE004001004 MULTIPOLYGON (((-13.24441 8...
## 5      159      Central I   SL040201      Central I Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.06882179 0.0001882636        SLE     SLE004  SLE004002 SLE004002001 MULTIPOLYGON (((-13.22646 8...
## 6      160         East I   SL040203         East I Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.05749169 0.0001429506        SLE     SLE004  SLE004002 SLE004002003 MULTIPOLYGON (((-13.2129 8....

# 
# DT::datatable(head(sle_adm3, 10), rownames = FALSE, 
#               options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

2.5人口数据 Popluation Ddata

塞拉利昂：ADM3 的人口

这些数据可以再次通过 Epirhandbook R 包下载，使用 import() 来加载 .csv 文件。我们还将导入的文件传递给 clean_names() 以标准化列名语法。

# Population by ADM3
sle_adm3_pop <- import(here("data", "Epidemiologist_R", "sle_admpop_adm3_2020.csv")) %>%
  clean_names()

这是人口文件的样子,可以查看每个辖区的男性人口、女性人口、总人口以及按年龄组分列的人口列。

head(sle_adm3_pop)

##        adm0_en adm0_pcode  adm1_en adm1_pcode adm2_en adm2_pcode       adm3_en adm3_pcode female   male  total t_00_04 t_05_09 t_10_14 t_15_19 t_20_24 t_25_29 t_30_34 t_35_39 t_40_44 t_45_49 t_50_54 t_55_59 t_60_64 t_65_69 t_70_74 t_75_79 t_80plus
## 1 Sierra Leone         SL Southern       SL03      Bo     SL0301        Badjia   SL030101   4440   4630   9070    1201    1418    1084    1117     847     777     555     539     382     310     240     142     144      95      84      51       86
## 2 Sierra Leone         SL Southern       SL03      Bo     SL0301         Bagbo   SL030102  14275  14585  28859    3820    4513    3450    3556    2695    2471    1765    1713    1219     988     762     450     458     300     266     162      274
## 3 Sierra Leone         SL Southern       SL03      Bo     SL0301 Bagbwe(Bagbe)   SL030103  11645  11687  23331    3088    3648    2788    2874    2179    1998    1427    1385     984     798     615     363     370     243     215     130      222
## 4 Sierra Leone         SL Southern       SL03      Bo     SL0301       Bo Town   SL030191  93653 100759 194412   25732   30401   23239   23949   18154   16648   11891   11541    8208    6652    5127    3032    3087    2021    1796    1089     1847
## 5 Sierra Leone         SL Southern       SL03      Bo     SL0301         Boama   SL030104  25066  26037  51104    6763    7991    6109    6295    4772    4376    3125    3034    2157    1748    1348     797     812     531     472     287      486
## 6 Sierra Leone         SL Southern       SL03      Bo     SL0301    Bumpe Ngao   SL030105  24214  25154  49369    6535    7720    5901    6082    4610    4228    3019    2930    2084    1689    1302     770     784     513     456     277      469

# 
# DT::datatable(head(sle_adm3_pop, 10), rownames = FALSE, 
#               options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

2.6卫生设施 health Facilities

塞拉利昂：来自 OpenStreetMap 的卫生设施数据

我们使用 read_sf() 导入设施点 shapefile，再次清理列名，然后过滤以仅保留标记为“医院hospital”、“诊所clinic”或“医生doctors”的点。

sle_hf <- sf::read_sf(here("data", "Epidemiologist_R", "sle_hf.shp")) %>% 
  clean_names() %>%
  filter(amenity %in% c("hospital", "clinic", "doctors"))

head(sle_hf)

## Simple feature collection with 6 features and 11 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.264 ymin: 8.38842 xmax: -13.14276 ymax: 8.490478
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 12
##       osm_id source                                                                 addrfull building healthcare operatorty addrcity name                                  amenity  healthca_1 capacitype             geometry
##        <dbl> <chr>                                                                  <chr>    <chr>    <chr>      <chr>      <chr>    <chr>                                 <chr>    <chr>      <chr>               <POINT [°]>
## 1 3197529881 UNMEER ; https://data.hdx.rwlabs.org/dataset/ebola-treatment-centers   <NA>     <NA>     <NA>       <NA>       <NA>     China-SL Friendship Hospital, Jui     hospital <NA>       <NA>        (-13.14276 8.38842)
## 2 3197529888 UNMEER ; https://data.hdx.rwlabs.org/dataset/ebola-treatment-centers   <NA>     <NA>     <NA>       <NA>       <NA>     Lakka Hospital ETU                    hospital <NA>       <NA>            (-13.264 8.397)
## 3 3212073781 https://data.hdx.rwlabs.org/dataset/sierra-leone-ebola-care-facilities <NA>     <NA>     <NA>       <NA>       <NA>     Princess Christian Maternity Hospital hospital <NA>       <NA>       (-13.21935 8.489861)
## 4 3339906977 MSF-CH                                                                 <NA>     <NA>     <NA>       <NA>       <NA>     COMMUNITY CLINIC                      clinic   <NA>       <NA>       (-13.16954 8.441156)
## 5 3341831443 MSF-CH                                                                 <NA>     <NA>     <NA>       <NA>       <NA>     Den Clinic                            clinic   <NA>       <NA>       (-13.24653 8.483513)
## 6 3341855623 <NA>                                                                   <NA>     <NA>     <NA>       <NA>       <NA>     MABELL HEALTH CENTER                  clinic   <NA>       <NA>       (-13.22632 8.490479)

# DT::datatable(head(sle_hf, 10), rownames = FALSE, 
#               options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

3.画坐标系 Plotting coordinates

这种情况下，绘制 X-Y 坐标（经度/纬度、点）的最简单方法是直接从我们在准备部分创建的 linelist_sf 对象中将它们绘制为点。

tmap 包为静态static（“绘图”模式）和交互式interactive（“视图”模式）提供了简单的映射功能，只需几行代码。 tmap 的语法与 ggplot2 的语法类似，详情点击。

1. 设置 tmap 模式。在这种情况下，我们将使用“plot”模式，它会产生静态static输出。

# view动态 or plot静态
tmap_mode("plot")

下面，这些点是单独绘制的。tm_shape() 与 linelist_sf 对象一起提供。然后我们通过 tm_dots() 添加点，指定大小和颜色。因为 linelist_sf 是一个 sf 对象，我们已经指定了包含纬度/经度坐标和坐标参考系（CRS）的两列：

# 绘制点
tm_shape(linelist_sf) + tm_dots(size = 0.08, col = 'blue')

单独来看，这些点并不能告诉我们太多。所以我们还应该绘制行政边界：

我们再次使用 tm_shape()，但我们不提供案例点 shapefile，而是提供管理边界 shapefile（多边形）。

使用 bbox = argument（bbox 代表“边界框”），我们可以指定坐标边界。首先我们显示没有bbox的地图显示，然后再使用它。

# Just the administrative boundaries (polygons)
# # 只是行政边界（多边形）
tm_shape(sle_adm3) +               # admin boundaries shapefile
  tm_polygons(col = "#F7F7F7")+    # show polygons in light grey
  tm_borders(col = "#000000",      # show borders with color and line weight
             lwd = 2) +
  tm_text("admin3name")            # column text to display for each polygon

# Same as above, but with zoom from bounding box
# # 同上，但从边界框缩放
tm_shape(sle_adm3,
         bbox = c(-13.3, 8.43,    # corner 左上角
                  -13.2, 8.5)) +  # corner
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")

现在把点和多边形边界放在一起

# All together
tm_shape(sle_adm3, bbox = c(-13.3, 8.43, -13.2, 8.5)) +     # 边界
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")+
tm_shape(linelist_sf) + 
  tm_dots(size=0.08, col='blue', alpha = 0.5) + # 点
  tm_layout(title = "Distribution of Ebola cases")   # give title to map

3.空间连接 Spatial joins

您可能熟悉将数据从一个数据集连接到另一个数据集。空间连接具有类似的目的，但利用了空间关系。您可以利用它们的空间关系，而不是依赖列中的公共值来正确匹配观测值，例如一个要素在另一个要素内，或与另一个要素最近的邻居，或在距另一个要素一定半径的缓冲区内，等等。

sf 包提供了多种空间连接方法。请参阅此参考中有关 st_join 方法和空间连接类型的更多文档。

3.1多边形中的点 Points in polygon

空间分配行政单位到案例 Spatial assign administrative units to cases

这是一个有趣的难题：案例行列表不包含有关案例行政单位的任何信息。尽管在初始数据收集阶段收集此类信息是理想的，但我们也可以根据其空间关系（即点与多边形相交）将行政单元分配给各个案例。

下面，我们将在空间上将案例位置（points）与 ADM3 边界（多边形polygons）相交：

1. 从行列表（点）开始
2. 空间连接到边界，将连接类型join设置为“st_intersects”
3. 使用 select() 仅保留某些新的管理边界列

linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects)

sle_adms 中的所有列都已添加到行列表中！现在，每个案例都有列详细说明其所属的管理级别。在此示例中，我们只想保留两个新列（admin level 3），因此我们select(）旧列名和另外两个感兴趣的列：

linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects) %>% 
  
  # Keep the old column names and two new admin ones of interest
  select(names(linelist_sf), admin3name, admin3pcod)

下面，仅出于显示目的，您可以看到前十个案例以及已附加的管理级别 3 (ADM3) 管辖区，基于点在空间上与多边形形状相交的位置。

# Now you will see the ADM3 names attached to each case
linelist_adm %>% select(case_id, admin3name, admin3pcod)

## Simple feature collection with 1000 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.27095 ymin: 8.447961 xmax: -13.20589 ymax: 8.490442
## Geodetic CRS:  WGS 84
## First 10 features:
##      case_id     admin3name admin3pcod                   geometry
## 5680  78703f       West III   SL040208 POINT (-13.25821 8.453866)
## 3020  525a03 Mountain Rural   SL040102 POINT (-13.22085 8.453176)
## 5154  19d49c       West III   SL040208 POINT (-13.26978 8.480407)
## 3963  76a224        East II   SL040204 POINT (-13.22583 8.484195)
## 2268  5355c3        East II   SL040204 POINT (-13.21946 8.481999)
## 3539  185ac2        West II   SL040207 POINT (-13.25619 8.483108)
## 5418  798126       West III   SL040208  POINT (-13.26513 8.47565)
## 3874  e3ba2c        West II   SL040207 POINT (-13.23237 8.460837)
## 2481  b345c5        West II   SL040207 POINT (-13.22861 8.469415)
## 2482  3ca785       West III   SL040208 POINT (-13.25258 8.457118)

现在我们可以按行政单位来描述我们的案例——这是在空间连接之前我们无法做到的！

# Make new dataframe containing counts of cases by administrative unit
# 创建包含行政单位案件计数的新数据框

case_adm3 <- linelist_adm %>%          # begin with linelist with new admin cols
  as_tibble() %>%                      # convert to tibble for better display
  group_by(admin3pcod, admin3name) %>% # group by admin unit, both by name and pcode 
  summarise(cases = n()) %>%           # summarize and count rows
  arrange(desc(cases))                     # arrange in descending order

head(case_adm3)

## # A tibble: 6 × 3
## # Groups:   admin3pcod [6]
##   admin3pcod admin3name     cases
##   <chr>      <chr>          <int>
## 1 SL040102   Mountain Rural   281
## 2 SL040208   West III         214
## 3 SL040207   West II          204
## 4 SL040204   East II          108
## 5 SL040201   Central I         54
## 6 SL040203   East I            54

我们还可以按行政单位创建病例计数条形图 bar plot。

在此示例中，我们以 linelist_adm 开始 ggplot()，以便我们可以应用 fct_infreq() 之类的因子函数，它按频率对条形图进行排序（有关提示，请参阅因子页面）。

ggplot(
    data = linelist_adm,                       # begin with linelist containing admin unit info
    mapping = aes(
      x = fct_rev(fct_infreq(admin3name))))+ # x-axis is admin units, ordered by frequency (reversed)
  geom_bar()+                                # create bars, height is number of rows
  coord_flip()+                              # flip X and Y axes for easier reading of adm units
  theme_classic()+                           # simplify background
  labs(                                      # titles and labels
    x = "Admin level 3",
    y = "Number of cases",
    title = "Number of cases, by adminstative unit",
    caption = "As determined by a spatial join, from 1000 randomly sampled cases from linelist"
  )

3.2最近的邻居 Nearest neighboor

寻找最近的医疗机构/集水区

了解与疾病热点相关的卫生设施的位置可能很有用。

我们可以使用 st_join() 函数（sf 包）中的 st_nearest_feature 连接方法来可视化离个别病例最近的医疗机构。

1. 我们从 shapefile linelist linelist_sf 开始
2. 我们在空间上加入 sle_hf，这是卫生设施和诊所的位置（点）

# Closest health facility to each case
linelist_sf_hf <- linelist_sf %>%                  # begin with linelist shapefile  
  st_join(sle_hf, join = st_nearest_feature) %>%   # data from nearest clinic joined to case data 
  select(case_id, osm_id, name, amenity) %>%       # keep columns of interest, including id, name, type, and geometry of healthcare facility
  rename("nearest_clinic" = "name")                # re-name for clarity

head(linelist_sf_hf)

## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.26978 ymin: 8.453176 xmax: -13.21946 ymax: 8.484195
## Geodetic CRS:  WGS 84
##      case_id     osm_id                  nearest_clinic  amenity                   geometry
## 5680  78703f 3341831443                      Den Clinic   clinic POINT (-13.25821 8.453866)
## 3020  525a03 6979628967 Shriners Hospitals for Children hospital POINT (-13.22085 8.453176)
## 5154  19d49c 3341831443                      Den Clinic   clinic POINT (-13.26978 8.480407)
## 3963  76a224 4540108136                       panasonic   clinic POINT (-13.22583 8.484195)
## 2268  5355c3 3342039261   GINER HALL COMMUNITY HOSPITAL   clinic POINT (-13.21946 8.481999)
## 3539  185ac2 3341831443                      Den Clinic   clinic POINT (-13.25619 8.483108)

# DT::datatable(head(linelist_sf_hf, 10), rownames = FALSE, 
#               options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

我们可以看到，约 30% 的病例中，“Den Clinic”是最近的医疗机构。

# Count cases by health facility
hf_catchment <- linelist_sf_hf %>%   # begin with linelist including nearest clinic data
  as.data.frame() %>%                # convert from shapefile to dataframe
  count(nearest_clinic,              # count rows by "name" (of clinic)
        name = "case_n") %>%         # assign new counts column as "case_n"
  arrange(desc(case_n))              # arrange in descending order

hf_catchment                         # print to console

##                          nearest_clinic case_n
## 1                            Den Clinic    369
## 2       Shriners Hospitals for Children    321
## 3         GINER HALL COMMUNITY HOSPITAL    176
## 4                             panasonic     43
## 5 Princess Christian Maternity Hospital     33
## 6                                  <NA>     23
## 7                  MABELL HEALTH CENTER     22
## 8                     ARAB EGYPT CLINIC     13

# tmap_mode("view")   # set tmap mode to interactive  

# plot the cases and clinic points 
tm_shape(linelist_sf_hf) +            # plot cases
  tm_dots(size=0.08,                  # cases colored by nearest clinic
          col='nearest_clinic') +    
tm_shape(sle_hf) +                    # plot clinic facilities in large black dots
  tm_dots(size=0.3, col='black', alpha = 0.4) +      
  tm_text("name") +                   # overlay with name of facility
tm_view(set.view = c(-13.2284, 8.4699, 13), # adjust zoom (center coords, zoom)
        set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")

3.3缓冲区 Buffers

我们还可以探索有多少病例位于距离最近的医疗机构 2.5 公里（约 30 分钟）的步行距离内。

注意：为了更准确地计算距离，最好将您的 sf 对象重新投影到相应的本地地图投影系统，例如 UTM（地球投影到平面上）。在此示例中，为简单起见，我们将坚持使用世界大地测量系统 (WGS84) 地理坐标系（地球以球形/圆形表面表示，因此单位为十进制度）。我们将使用一般转换：1 十进制度 = ~111km。

在这篇 esri 文章中查看有关地图投影和坐标系的更多信息。该博客讨论了不同类型的地图投影，以及如何根据感兴趣的区域和地图/分析的背景选择合适的投影。

1. 首先，在每个卫生设施周围创建一个半径约为 2.5 公里的圆形缓冲区。这是通过 tmap 中的函数 st_buffer() 完成的。因为 地图的单位是纬度/经度十进制度，这就是“0.02”的解释方式。如果您的地图坐标系以米为单位，则必须以米为单位提供数字。

sle_hf_2k <- sle_hf %>%
  st_buffer(dist=0.02)       # decimal degrees translating to approximately 2.5km 

tmap_mode("plot")
# Create circular buffers
# 创建循环缓冲区
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) +                    # plot clinic facilities in large red dots
  tm_dots(size=0.3, col='black')      

2. 其次，我们使用 st_join() 和 st_intersects* 的连接类型将这些缓冲区与案例（点）相交。也就是说，来自缓冲区的数据连接到它们相交的点。

# Intersect the cases with the buffers
linelist_sf_hf_2k <- linelist_sf_hf %>%
  st_join(sle_hf_2k, join = st_intersects, left = TRUE) %>%
  filter(osm_id.x== osm_id.y | is.na(osm_id.y)) %>%
  select(case_id, osm_id.x, nearest_clinic, amenity.x, osm_id.y)

现在我们可以计算结果：nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y),0])在 1000 个案例中没有与任何缓冲区相交（缺少该值），因此步行超过 30 分钟最近的医疗机构。

# Cases which did not get intersected with any of the health facility buffers
linelist_sf_hf_2k %>% 
  filter(is.na(osm_id.y)) %>%
  nrow()

## [1] 1000

我们可以将结果可视化，使得不与任何缓冲区相交的案例显示为红色。

# tmap_mode("view")

# First display the cases in points 首先以点显示案例
tm_shape(linelist_sf_hf) +
  tm_dots(size=0.08, col='nearest_clinic') +

# plot clinic facilities in large black dots 用大黑点绘制诊所设施
tm_shape(sle_hf) +                    
  tm_dots(size=0.3, col='black')+   

# Then overlay the health facility buffers in polylines 然后在折线中覆盖卫生设施缓冲区
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2) +

# Highlight cases that are not part of any health facility buffers 用红点突出显示不属于任何卫生设施缓冲区的病例
# in red dots  
tm_shape(linelist_sf_hf_2k %>%  filter(is.na(osm_id.y))) +
  tm_dots(size=0.1, col='red') +
tm_view(set.view = c(-13.2284,8.4699, 13), set.zoom.limits = c(13,14))+

# add title  
tm_layout(title = "Cases by clinic catchment area")

3.4其他空间连接类型

参数join的替代值包括文档:

• st_contains_properly
• st_contains
• st_covered_by
• st_covers
• st_crosses
• st_disjoint
• st_equals_exact
• st_equals
• st_is_within_distance
• st_nearest_feature
• st_overlaps
• st_touches
• st_within

4.等值域图 Choropleth 地图

Choropleth 地图可用于按预定义区域（通常是行政单位或卫生区域）可视化您的数据。例如，在疫情应对中，这有助于针对高发病率的特定地区分配资源。

现在我们已经为所有案例分配了行政单位名称（请参阅上面的空间连接部分），我们可以开始按区域映射案例计数（等值域图）。由于我们还有 ADM3 的人口数据，我们可以将此信息添加到之前创建的 case_adm3 表中。

我们从上一步中创建的数据框 case_adm3 开始，它是每个行政单位及其案例数量的汇总表。

1. 人口数据 sle_adm3_pop 使用 dplyr 的 left_join() 连接，基于 case_adm3 数据帧中的 admin3pco 列和 sle_adm3_pop 数据帧中的 adm_pcode 列的共同值。请参阅加入数据页面）。
2. select() 应用于新的数据框，只保留有用的列 - total 是总人口
3. 每 10,000 人口中的病例数计算为带有 mutate() 的新列

# Add population data and calculate cases per 10K population
case_adm3 <- case_adm3 %>% 
     left_join(sle_adm3_pop,                             # add columns from pop dataset
               by = c("admin3pcod" = "adm3_pcode")) %>%  # join based on common values across these two columns 基于这两列的共同值连接
     select(names(case_adm3), total) %>%                 # keep only important columns, including total population
     mutate(case_10kpop = round(cases/total * 10000, 3)) # make new column with case rate per 10000, rounded to 3 decimals

case_adm3                                                # print to console for viewing

## # A tibble: 10 × 5
## # Groups:   admin3pcod [10]
##    admin3pcod admin3name     cases  total case_10kpop
##    <chr>      <chr>          <int>  <int>       <dbl>
##  1 SL040102   Mountain Rural   281  33993       82.7 
##  2 SL040208   West III         214 210252       10.2 
##  3 SL040207   West II          204 145109       14.1 
##  4 SL040204   East II          108  99821       10.8 
##  5 SL040201   Central I         54  69683        7.75
##  6 SL040203   East I            54  68284        7.91
##  7 SL040206   West I            42  60186        6.98
##  8 SL040202   Central II        25  23874       10.5 
##  9 SL040205   East III          11 500134        0.22
## 10 <NA>       <NA>               7     NA       NA

将此表与 ADM3 多边形 shapefile 连接以进行映射

case_adm3_sf <- case_adm3 %>%                 # begin with cases & rate by admin unit
  left_join(sle_adm3, by="admin3pcod") %>%    # join to shapefile data by common column
  select(objectid, admin3pcod,                # keep only certain columns of interest
         admin3name = admin3name.x,           # clean name of one column
         admin2name, admin1name,
         cases, total, case_10kpop,
         geometry) %>%                        # keep geometry so polygons can be plotted
  drop_na(objectid) %>%                       # drop any empty rows
  st_as_sf()                                  # convert to shapefile

# tmap mode
tmap_mode("plot")               # view static map

# plot polygons
tm_shape(case_adm3_sf) + 
        tm_polygons("cases") +  # color by number of cases column
        tm_text("admin3name")   # name display

我们还可以绘制发病率图

# Cases per 10K population
tmap_mode("plot")             # static viewing mode

# plot
tm_shape(case_adm3_sf) +                # plot polygons
  tm_polygons("case_10kpop",            # color by column containing case rate
              breaks=c(0, 10, 50, 100), # define break points for colors
              palette = "Purples"       # use a purple color palette
              ) +
  tm_text("admin3name")                 # display text

5.使用 ggplot2

进行映射如果您已经熟悉使用 ggplot2，则可以使用该包来创建数据的静态映射。 geom_sf() 函数将根据数据中的特征（点、线或多边形）绘制不同的对象。例如，您可以在 ggplot() 中使用 geom_sf()，使用 sf 数据和多边形几何来创建等值线图。

为了说明这是如何工作的，我们可以从我们之前使用的 ADM3 多边形 shapefile 开始。回想一下，这些是塞拉利昂的 3 级管理员区域：

sle_adm3

## Simple feature collection with 12 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.29894 ymin: 8.094272 xmax: -12.91333 ymax: 8.499809
## Geodetic CRS:  WGS 84
## First 10 features:
##    objectid     admin3name admin3pcod    admin3ref_n         admin2name admin2pcod admin1name admin1pcod   admin0name admin0pcod       date   valid_on valid_to shape_leng   shape_area rowcacode0 rowcacode1 rowcacode2   rowcacode3                       geometry
## 1       155     Koya Rural   SL040101     Koya Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.63822264 1.366013e-02        SLE     SLE004  SLE004001 SLE004001001 MULTIPOLYGON (((-13.02082 8...
## 2       156 Mountain Rural   SL040102 Mountain Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.29268365 3.182343e-03        SLE     SLE004  SLE004001 SLE004001002 MULTIPOLYGON (((-13.21496 8...
## 3       157 Waterloo Rural   SL040103 Waterloo Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.72265524 1.364140e-02        SLE     SLE004  SLE004001 SLE004001003 MULTIPOLYGON (((-13.12215 8...
## 4       158     York Rural   SL040104     York Rural Western Area Rural     SL0401    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 1.23916989 1.983448e-02        SLE     SLE004  SLE004001 SLE004001004 MULTIPOLYGON (((-13.24441 8...
## 5       159      Central I   SL040201      Central I Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.06882179 1.882636e-04        SLE     SLE004  SLE004002 SLE004002001 MULTIPOLYGON (((-13.22646 8...
## 6       160         East I   SL040203         East I Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.05749169 1.429506e-04        SLE     SLE004  SLE004002 SLE004002003 MULTIPOLYGON (((-13.2129 8....
## 7       161        East II   SL040204        East II Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.08397463 1.494827e-04        SLE     SLE004  SLE004002 SLE004002004 MULTIPOLYGON (((-13.22653 8...
## 8       162     Central II   SL040202     Central II Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.04878431 6.513035e-05        SLE     SLE004  SLE004002 SLE004002002 MULTIPOLYGON (((-13.23154 8...
## 9       163       West III   SL040208       West III Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.30219417 1.698296e-03        SLE     SLE004  SLE004002 SLE004002008 MULTIPOLYGON (((-13.28529 8...
## 10      164         West I   SL040206         West I Western Area Urban     SL0402    Western       SL04 Sierra Leone         SL 2016-08-01 2016-10-17     <NA> 0.06952030 1.822183e-04        SLE     SLE004  SLE004002 SLE004002006 MULTIPOLYGON (((-13.24677 8...

我们可以使用 dplyr 的 left_join() 函数来添加我们想要映射到 shapefile 对象的数据。在这种情况下，我们将使用之前创建的 case_adm3 数据框来汇总行政区域的病例数；但是，我们可以使用相同的方法来映射存储在数据框中的任何数据。

sle_adm3_dat <- sle_adm3 %>% 
  inner_join(case_adm3, by = "admin3pcod") # inner join = retain only if in both data objects

select(sle_adm3_dat, admin3name.x, cases) # print selected variables to console

## Simple feature collection with 9 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.29894 ymin: 8.384533 xmax: -13.12612 ymax: 8.499809
## Geodetic CRS:  WGS 84
##     admin3name.x cases                       geometry
## 1 Mountain Rural   281 MULTIPOLYGON (((-13.21496 8...
## 2      Central I    54 MULTIPOLYGON (((-13.22646 8...
## 3         East I    54 MULTIPOLYGON (((-13.2129 8....
## 4        East II   108 MULTIPOLYGON (((-13.22653 8...
## 5     Central II    25 MULTIPOLYGON (((-13.23154 8...
## 6       West III   214 MULTIPOLYGON (((-13.28529 8...
## 7         West I    42 MULTIPOLYGON (((-13.24677 8...
## 8        West II   204 MULTIPOLYGON (((-13.25698 8...
## 9       East III    11 MULTIPOLYGON (((-13.20465 8...

要按区域制作病例计数柱形图，使用 ggplot2，我们可以调用 geom_col()，如下所示：

ggplot(data=sle_adm3_dat) +
  geom_col(aes(x=fct_reorder(admin3name.x, cases, .desc=T),   # reorder x axis by descending 'cases' 
                                                              #通过降序“案例”对 x 轴重新排序

               y=cases)) +                                  # y axis is number of cases by region
  theme_bw() +
  labs(                                                     # set figure text
    title="Number of cases, by administrative unit",
    x="Admin level 3",
    y="Number of cases"
  ) + 
  guides(x=guide_axis(angle=45))                            # angle x-axis labels 45 degrees to fit better

如果我们想使用 ggplot2 来制作案例计数的等值线图，我们可以使用类似的语法来调用 geom_sf() 函数：

ggplot(data=sle_adm3_dat) + 
  geom_sf(aes(fill=cases))    # set fill to vary by case count variable

然后，我们可以使用与 ggplot2 一致的语法来自定义地图的外观，例如：

ggplot(data=sle_adm3_dat) +                           
  geom_sf(aes(fill=cases)) +                        
  scale_fill_continuous(high="#54278f", low="#f2f0f7") +    # change color gradient
  theme_bw() +
  labs(title = "Number of cases, by administrative unit",   # set figure text
       subtitle = "Admin level 3"
  )

对于习惯使用 ggplot2 的 R 用户，geom_sf() 提供了一个简单直接的实现，适用于基本的地图可视化。要了解更多信息，请阅读 geom_sf() 小插图或 ggplot2 书。

6.底图 Basemaps

6.1 OpenStreetMap

下面我们将描述如何使用 OpenStreetMap 功能为 ggplot2 地图实现底图。替代方法包括使用 ggmap，它需要在 Google 上免费注册（详情）。

OpenStreetMap 是一个协作项目，旨在创建免费的可编辑世界地图。基础地理位置数据（例如城市位置、道路、自然特征、机场、学校、医院、道路等）被认为是项目的主要输出。

首先我们加载 OpenStreetMap 包，我们将从中获取我们的底图。

然后，我们创建对象映射，我们使用 OpenStreetMap 包（文档）中的函数 openmap() 定义它。我们提供以下内容：

• upperLeft 和 lowerRight 两个坐标对指定底图切片的范围 在这种情况下，我们从 linelist 行中放入了 max 和 min，因此地图将动态响应数据
• zoom = (如果为 null 它自动确定）
• type = 底图的类型 - 我们列出了几个 这里的可能性和代码当前使用第一个（[1]）“osm”
• mergeTiles = 我们选择了TRUE，所以底图都合并为一个

# load package
pacman::p_load(OpenStreetMap)

# Fit basemap by range of lat/long coordinates. Choose tile type
map <- OpenStreetMap::openmap(
  upperLeft = c(max(linelist$lat, na.rm=T), max(linelist$lon, na.rm=T)),   # limits of basemap tile
  lowerRight = c(min(linelist$lat, na.rm=T), min(linelist$lon, na.rm=T)),
  zoom = NULL,
  type = c("osm", "stamen-toner", "stamen-terrain", "stamen-watercolor", "esri","esri-topo")[1])

如果我们现在使用 OpenStreetMap 包中的 autoplot.OpenStreetMap() 绘制此底图，您会看到轴上的单位不是纬度/经度坐标。它使用不同的坐标系。要正确显示案例住所（以纬度/经度存储），必须更改此设置。

autoplot.OpenStreetMap(map)

因此，我们想使用 OpenStreetMap 包中的 openproj() 函数将地图转换为纬度/经度。我们提供底图地图，还提供我们想要的坐标参考系 (CRS)。我们通过为 WGS 1984 投影提供“proj.4”字符串来做到这一点，但您也可以通过其他方式提供 CRS。 （请参阅此页面以更好地了解 proj.4 字符串是什么）

# Projection WGS84
map_latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

现在，当我们创建绘图时，我们看到沿轴是纬度和经度坐标。坐标系已转换。现在，如果重叠，我们的案例将正确绘制！

# Plot map. Must use "autoplot" in order to work with ggplot
autoplot.OpenStreetMap(map_latlon)

7.等高线密度热图 Contoured density heatmaps

下面我们描述了如何在底图上实现案例的等高线密度热图，从一个线列表（每个案例一行）开始。

1. 从 OpenStreetMap 创建底图切片，如上所述
2. 使用 latitude 和 longitude 列从 linelist 中绘制案例
3. 使用 ggplot2 中的 stat_density_2d()将点转换为密度热图

当我们有一个具有纬度/经度坐标的底图时，我们可以使用其居住地的经度/纬度坐标将我们的案例绘制在顶部。基于函数 autoplot。 OpenStreetMap() 创建底图，ggplot2 函数将很容易添加到顶部，如下所示的 geom_point()：

# Plot map. Must be autoplotted to work with ggplot
autoplot.OpenStreetMap(map_latlon)+                 # begin with the basemap
  geom_point(                                       # add xy points from linelist lon and lat columns 
    data = linelist,                                
    aes(x = lon, y = lat),
    size = 1, 
    alpha = 0.5,
    show.legend = FALSE) +                          # drop legend entirely
  labs(x = "Longitude",                             # titles & labels
       y = "Latitude",
       title = "Cumulative cases")

上面的地图可能难以解释，尤其是在点重叠的情况下。因此，您可以改为使用 ggplot2 函数 stat_density_2d() 绘制 2d 密度图。您仍在使用 linelist 纬度/经度坐标，但会执行 2D 内核密度估计，并使用等高线显示结果 - 就像地形图一样。在此处阅读完整文档。

# begin with the basemap
autoplot.OpenStreetMap(map_latlon)+
  
  # add the density plot
  ggplot2::stat_density_2d(
        data = linelist,
        aes(
          x = lon,
          y = lat,
          fill = ..level..,
          alpha = ..level..),
        bins = 10,
        geom = "polygon",
        contour_var = "count",
        show.legend = F) +                          
  
  # specify color scale
  scale_fill_gradient(low = "black", high = "red")+
  
  # labels 
  labs(x = "Longitude",
       y = "Latitude",
       title = "Distribution of cumulative cases")

7.2 时间序列热图 Time series heatmap

上面的密度热图显示了累积案例cumulative cases。我们可以通过基于症状发作月份的热图来检查随着时间和空间的爆发，从线条列表中派生出来。

我们从linelist开始，创建一个带有发病年份和月份的新列。 base R 中的 format() 函数改变了日期的显示方式。在这种情况下，我们想要“YYYY-MM”。

linelist = linelist %>% 
  mutate(date_onset_ym = format(date_onset, "%Y-%m"))
# Examine the values 
table(linelist$date_onset_ym, useNA = "always")

## 
## 2014-04 2014-05 2014-06 2014-07 2014-08 2014-09 2014-10 2014-11 2014-12 2015-01 2015-02 2015-03 2015-04    <NA> 
##       3       7      13      38      77     181     204     134     100      73      54      42      32      42

现在，我们简单地通过 ggplot2 将 facetting 引入密度热图。应用 facet_wrap()，将新列用作行。为了清楚起见，我们将分面列的数量设置为 3。

# packages
pacman::p_load(OpenStreetMap, tidyverse)

# begin with the basemap
autoplot.OpenStreetMap(map_latlon)+
  
  # add the density plot
  ggplot2::stat_density_2d(
        data = linelist,
        aes(
          x = lon,
          y = lat,
          fill = ..level..,
          alpha = ..level..),
        bins = 10,
        geom = "polygon",
        contour_var = "count",
        show.legend = F) +                          
  
  # specify color scale
  scale_fill_gradient(low = "black", high = "red")+
  
  # labels 
  labs(x = "Longitude",
       y = "Latitude",
       title = "Distribution of cumulative cases over time")+
  
  # facet the plot by month-year of onset
  facet_wrap(~ date_onset_ym, ncol = 4)               

Reference

the handbook team. (2021). R for applied epidemiology and public health | The Epidemiologist R Handbook. Epirhandbook.com. https://epirhandbook.com/en/index.html
tmap: get started! (2020). R-Project.org. https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html
Geometric binary predicates on pairs of simple feature geometry sets — geos_binary_pred. (2017). Github.io. https://r-spatial.github.io/sf/reference/geos_binary_pred.html
Visualise sf objects — CoordSf. (2022). Tidyverse.org. https://ggplot2.tidyverse.org/reference/ggsf.html

‌

leaflet包

LeafletR主要是用R语言的语法封装了JS版的Leaflet，可以在R语言的plot窗口，利用html5技术显示各种地图，还可以绘制自己的要素图形。

它有如下功能：
交互地图浏览（缩放、平移）
使用多种底图进行任意组合
加载地图瓦片（WMTS）
点要素定位标记
多边形要素标记
线要素标记
弹出窗口
解析加载GeoJson
从R或者RSutido创建地图窗口
可以把地图嵌入 knitr/R包所生成的Markdown文档中，或者是Shiny制作的APP中。
可以直接获取通过SP包生成就（加载）的空间对象以及包含经纬度的数据框进行展示。
可以设定地图范围以及封装自定义的鼠标事件。

一般来说，它的基本使用步骤如下：

1、加载leaflet包
2、通过leaflet包创建地图控件。
3、通过图层操作的方法（如addTiles、 addMarkers、 addPolygons）来处理图层数据，并且修改地图插件的各种参数，来把图层显示在地图控件上。
4、可以重复第三步，可以增加更多的图层数据。
5、把地图部件显示出来。完成绘图。

下面来看一个例子：

m <- leaflet()
at <- addTiles(m)
addMarkers(at,lng=116.391, lat=39.912, popup="这里是北京")

leaflet()%>%addTiles()%>%addMarkers(lng=116.391, lat=39.912, popup="这里是北京")

官方的文档里面，对于这种需要重复进行容器嵌套的写法，提供了管道操作符“%>%”来实现，它的主要作用就是把前面的语句（变量）传递给下一个语句，并且作为第一个参数使用，上面那三个语句，利用管道操作符来写，如下,输出结果是完全一样的，但是语句变得简单了。

二、地图控件

在leaflet包初始化的时候，一般调用leaflet()这个方法，这个方法就是对地图控件进行初始化，会生成一个地图容器，以后所有的图层操作，都在这个容器内处理。一般来说，这个方法都被作为其他方法的第一个参数来使用。我们可以通过显示参数设定或者通过管道操作符%>%来把这个容器传递给其他的方法。

地图控件的基本方法有下面这三个：

setView() :设定地图的显示级别缩放比例、和地图的中心点。
fitBounds()：设定地图的范围，一般是一个矩形，结构是：[lng1, lat1] – [lng2, lat2]。
clearBounds()：清除地图的范围设定。

下面来看个例子：

m<- leaflet()
m<- setView(m,lng=116.38,lat=39.9,zoom=9)
addTiles(m)

或者用管道操作符

leaflet()%>%setView(lng=116.38,lat=39.9,zoom=9)%>%addTiles()

初始化改为3级的话，结果如下：

leaflet()%>%setView(lng=116.38,lat=39.9,zoom=3)%>%addTiles()

leaflet包支持各种与空间信息有关的对象，包括使用sp包定义的空间对象，和R语言中带有空间信息的数据框等，如下所示：
与R相关的：

由经纬度信息组成矩阵
带有经纬度字段的数据框。

与sp包相关的：

SpatialPoints[DataFrame]
Line/Lines
SpatialLines[DataFrame]
Polygon/Polygons
SpatialPolygons[DataFrame]

还有就是maps包里面的各种空间图形信息

这些对象都可以直接用于leaflet里面的方法，作为图层添加到地图上。下面我们通过leaflet来绘制一些地图

当然，最简单的，还是绘制点，通过经纬度，把点画上去，我们先生成一批随机点：

df = data.frame(Lat = rnorm(100), Lon = rnorm(100))

然后画上去：

m <-leaflet(df)
addCircles(m)

#写法2：

df %>%leaflet()%>%addCircles()

leaflet包相对其它地图包，有很多优点和缺点，

首先，绘制地图简单快捷，因为都是基于供应商的tiles，一行代码就可以render出基本widget地图。

支持管道传参，一个图层一个图层进行添加，代码结构更加清晰。

其次，有很多tiles供应商可以选择,包括高德、google、Stamen, Esri, OpenWeatherMap，NASA，

等好几十个tiles供应商。当然其中一些需要注册。其中的google可以绕过注册，已经很难得了。

对于在地图上添加markers图标，shapes形状，线条等，异常方便快捷，这在ggplot2中很难做到的。

支持栅格数据，rasters栅格数据是基于像素点的地图。可以看出，leaflet具有很强的包容性。

支持多种投影坐标系，甚至可以自定义坐标系，这在某些特殊场景非常重要。

当然还有更重要的是，其具有一定交互能力，可以缩放拖拽，

简单的图层切换也不需要使用Shiny。使得更容易上手。

其它特点，首先tiles是基于供应商的，必须联网，

其次对颜色支持不一样，只支持HEX颜色空间和colors()中的颜色名称。 当然内置的几个palette函数，非常特别。

1.Widget设置

Widget地图框的设置，就是确定Widget的基本参数，

包括CRS坐标系，widget的中心坐标，zoom level(缩放)的范围, widget边界坐标，data数据等。

1.1

leafletOptions()

leaflet()中有个options参数，用leafletOptions()函数来指定，可以控制widget缩放范围。

语法：
1leafletOptions(minZoom = NULL, maxZoom = NULL, crs = leafletCRS(),
2 worldCopyJump = NULL, preferCanvas = NULL, …)

关键参数：
minZoom,表示最低缩小倍数，作用于所有地图层。
maxZoom,表示最高放大倍数，作用于所有地图层。
crs, 表示指定坐标系统，
preferCanvas, 表示是否将leaflet.js路径呈现在地图上。
1library(leaflet)
2
3leaflet(options= leafletOptions(minZoom=0, maxZoom =18))

1.2

中心、缩放、边界

关键函数：
setView() ,设定地图的view(包括center位置和zoom level)
flyTo() ,切换到一个指定的location或zoom-level，使用光滑的pan-zoom
fitBounds() ,设定地图矩形区域边界。view将限制在[lng1, lat1] - [lng2, lat2]
flyToBounds() ,切换到一个指定的地图矩形区域边界，使用光滑的pan/zoom
setMaxBounds() ,限定地图矩形区域最大边界
clearBounds() ,清除地图矩形区域边界, 然后view将只会受地图层的经度和纬度数据限制。
语法：
setView(map, lng, lat, zoom, options = list())
flyTo(map, lng, lat, zoom, options = list())
fitBounds(map, lng1, lat1, lng2, lat2, options = list())
flyToBounds(map, lng1, lat1, lng2, lat2, options = list())
setMaxBounds(map, lng1, lat1, lng2, lat2)
ClearBounds (map)

参数解释：
map，表示leaflet()创建的map widget
lng, 表示map center的经度,东经为正
lat, 表示map center的纬度，北纬为正
zoom, 表示zoom level
options, 列表传参，传递zoom或pan参数。
lngl, latl, lng2, lat2, 表示widget边界的坐标。

1library(leaflet)
2
3# 设定中心坐标和zoom level
4m <- leaflet() %>% addTiles() %>% setView(-71.0382679,42.3489054, zoom =18)56
# 显示第一个view
7m %>% fitBounds(-72,40,-70,43) 
# 设定view边界89# 显示第二个view
10m %>% clearBounds() # 清除边界限制, leaflet()默认为世界地图

Data数据

这里的Data不仅仅是画地图上行政区域的数据，而且包括要在地图上呈现的数据。

大多数图层添加函数都有data参数，通常使用%>%管道符逐渐传递data参数。

leaflet()通常支持下列几种形式的数据。

矩阵数据(由经度和纬度构成)。

数据框(由经度和纬度构成)。

从sp包传递的数据，包括：

SpatialPoints(数据框类型)

Line()/Lines()

SpatialLines()(数据框类型)

Polygon()/Polygons()

SpatialPolygons()(数据框类型)

从maps包传递的数据，主要是map()函数传递的数据框。

对于经度和纬度组成矩阵或数据框类型数据，在调动data添加图层时，会根据变量名进行猜测匹配：

若变量名称为lat,或latitude等，则猜测为纬度，猜测时，不区分大小写。

若变量名称为lng, long或 longitude等，则猜测为经度，猜测时，不区分大小写。

也可以手动指定经度和纬度变量，使用~语法。

在参数传递过程中，默认后面的data参数覆盖前面data参数。

1.3.1 data中经纬度的指定/猜测/覆盖

1library(leaflet)23# 自动猜测匹配

4set.seed(123)5df <- data.frame(Lat = 1:10, Long = rnorm(10))6leaflet(df)%>%addCircles()78# 手动指定经度和纬度变量,结果一样 leaflet(df) %>% addCircles(lng = ~Long,

9# lat = ~Lat)

1011# 在add_xxx()函数中重新指定参数进行覆盖，结果一样 leaflet() %>%

12# addCircles(data = df) leaflet() %>% addCircles(data = df, lat = ~Lat, lng

13# = ~Long)

1.3.2 sp对象的data

1library(leaflet)
2library(sp)
3library(RColorBrewer)
4
5Sr1 <- Polygon(cbind(c(2,4,4,1,2), c(2,3,5,4,2))) 
#4对非重复点坐标，首尾相连
6Sr2 <- Polygon(cbind(c(5,4,2,5), c(2,3,2,2))) 
#3对非重复点坐标，画三角形
7Sr3 <- Polygon(cbind(c(4,4,5,10,4), c(5,3,2,5,5))) 
#4对非重复点坐标，画四边形
8Sr4 <- Polygon(cbind(c(5,6,6,5,5), c(4,4,3,3,4)), hole = TRUE) # hole = TRUE,表示中空
9
10
11Srs1 <- Polygons(list(Sr1),"s1") #'s1'指定ID参数，多个多边形才有ID参数，Polygon()没有ID参数
12Srs2 <- Polygons(list(Sr2),"s2")
13Srs3 <- Polygons(list(Sr4, Sr3),"s3/4") # 合并Sr3和Sr4多边形
14
15# 列表传参，传递多个多边形参数
16SpP <- SpatialPolygons(list(Srs1, Srs2, Srs3),1:3)
17
18leaflet(height ="300px") %>% addPolygons(data = SpP, fillColor = brewer.pal(3,19name ="Set1"))

1.3.3 从maps包中获取data

1library(leaflet)
2library(maps)
3library(RColorBrewer)
4
5mapStates <- map("state", fill = TRUE, plot = FALSE)
6leaflet(data= mapStates) %>% addTiles() %>% addPolygons(fillColor = brewer.pal(10,7 name="Paired"), stroke = FALSE)

1.3.4 其它参数

其它绘图参数支持R自带的数据类型，如：向量，颜色向量，数据框，同样支持用~指定。

1library(leaflet)
2
3# 随便编一个数据
4m <- leaflet() %>% addTiles()
5df <- data.frame(lat = rnorm(100), lng = rnorm(100),size= runif(100,5,20),
6color= sample(colors(),100))
7
8m <- leaflet(df) %>% addTiles()
9#circle图标的半径不随zoom变化:
10m %>% addCircleMarkers(radius = ~size,color= ~color, fill = FALSE)
11#circle图标的半径不随zoom变化:
12m %>% addCircleMarkers(radius = runif(100,4,10),color= c("red"))

2.Basemaps底图

leaflet支持Tilemap类型的底图，

leaflet支持多种免费第三方providers的tiles, 包括Stamen, Esri, OpenWeatherMap等，

用names(providers)可以查看所有的providers。

2.1

默认Tiles(OpenStreetMap)

使用addTiles()函数添加tiles并使用默认参数,默认即是OpenStreetMap，即街区Tiles。

1library(leaflet)
2
3m <- leaflet()%>%setView(lng = -71.0589, lat = 42.3601, zoom = 12)4m%>%addTiles()# 显示街区图

第三方 Tiles

调用addProviderTiles()函数，在参数providers$后面加tiles供应商的名字就行了。 需要注意：部分第三方tiles需要注册。

通过options参数调用providerTileOptions()函数可以规避部分tiles的注册。

如果有定制的tiles模板的URL链接,可以在addTiles()函数中调用。

1library(leaflet)
2
3m <- leaflet()%>%setView(lng = -71.0589, lat = 42.3601, zoom = 12)
4
5# 使用Stamen.Toner 的 tiles
6m%>%addProviderTiles(providers$Stamen.Toner)

巧妙的数据转换

通过调用函数addWMSTiles()可以添加WMS(Web Map Service)的tiles。

Web地图服务（WMS）是一种标准协议，描述如何通过Internet提供任何地理配准的地图图像，

这通常由使用来自地理信息系统数据库的数据的地图服务器生成。

协议标准由Open Geospatial Consortium（OGC）开发，

并于1999年首次发布.WMS提供了一种使用HTTP接口请求地理注册地图图像的简单方法。

WMS供应商（https://en.wikipedia.org/wiki/Web_Map_Service）

1library(leaflet)
2
3leaflet() %>% addTiles() %>% setView(-93.65,42.0285, zoom =4) %>% # 叠加一个WMS的tiles图层4addWMSTiles("http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi", layers ="nexrad-n0r-900913",
5 options = WMSTileOptions(format ="image/png", transparent = TRUE), attribution ="Weather data <U+00A9> 2012 IEM Nexrad")

Tiles图层叠加

多个Tiles图层也可以叠加， 但是这种情况通常仅用于表层tiles是半透明的情况下，

或在options参数中手动指定不透明度opacity。

opacity从0(完全透明)到1(完全不透明)。

1library(leaflet)
2
3m <- leaflet() %>% setView(lng = -71.0589, lat =42.3601, zoom =12)
4
5m %>% addProviderTiles(providers$MtbMap) %>% # 底层tiles
6 addProviderTiles(providers$Stamen.TonerLines, # 叠加一层tiles,显示公路和街道
7 options = providerTileOptions(opacity =0.35)) %>% # opacity设定非透明度
8 addProviderTiles(providers$Stamen.TonerLabels) # 叠加tiles,显示公路名，街道名，机场图标。
链接：https://blog.csdn.net/R3eE9y2OeFcU40/article/details/88324895