-
2022-05-27 14:48:10
最近需要对国内疫情分布情况绘制可视化地图,查找资料R中地图绘制思路,显示在R中绘制地图主要有三种方式:第一种是利用某些特定R包中自带的地图数据进行绘图;第二种从其他途径获取地理信息数据,调用相应的软件包对数据进行读取,进而绘图;第三种是基于某些供应商的tiles与Google、NASA、高德等网络在线地图相关联,调用其地图数据为自己绘图所用。下面进行举例说明:
1.【绘图前准备】爬取丁香园每日疫情数据
##加载程序包,设置路径## setwd("f://data") library(rvest) library(stringr) library(dplyr) library(ggplot2) library(leaflet) library(leafletCN) library(RColorBrewer) library(chinamap) library(sp) ###爬取丁香园网站国内34个省份疫情实时数据### url<-"https://ncov.dxy.cn/ncovh5/view/pneumonia" web_mess<-read_html(url) web_node<-web_mess%>%html_nodes("#getAreaStat") web_text<-html_text(web_node,trim=T) ###数据正则化处理### web_str<-str_extract_all(web_text,"provinceName\":\".{2,9}\",\"provinceShortName\":\".{2,9}\",\"currentConfirmedCount\":.{1,7},") pp<-str_extract_all(web_str[[1]],"provinceName\":\".{2,9}\"") province<-data.frame() for(i in 1:length(web_str[[1]])){ province[i,1]<-str_sub(pp[[i]][1],16,str_length(pp[[i]][1])-3) } cc<-str_extract_all(web_str[[1]],"\"currentConfirmedCount\":.{1,7},") num<-as.vector(34,mode = "numeric") for(i in 1:length(web_str[[1]])){ num[i]<-str_sub(cc[[i]][1],str_locate(cc[[i]][1],":")[,1]+1,str_locate(cc[[i]][1],",")[,1]-1) } mydata<-data.frame(province,as.numeric(num))%>%setNames(c("area","confirm")) case.cate<-function(x){###将数值型变量转换为分类变量函数### for(i in 1:length(x)){ if(x[i]==0){ x[i]<-"0" }else if(x[i]>0&x[i]<=10){ x[i]<-"1~10" }else if(x[i]>10&x[i]<=50){ x[i]<-"10~50" }else if(x[i]>50&x[i]<=100){ x[i]<-"50~100" }else if(x[i]>100&x[i]<=200){ x[i]<-"100~200" }else if(x[i]>200&x[i]<=500){ x[i]<-"200~500" }else if(x[i]>500&x[i]<=1000){ x[i]<-"500~1000" }else if(x[i]>1000&x[i]<=2000){ x[i]<-"1000~2000" }else if(x[i]>2000&x[i]<=300000){ x[i]<-"2000~300000" }else{ x[i]<-"1000000~" } } return(x) } lapply(mydata$confirm,case.cate)%>%as.character()%>%as.factor()->mydata$nvalue substr(mydata$area,1,2)->mydata$id_area##最终的绘图数据##
2. 结合前两种方式绘制静态地图
第一种方式是利用Y叔提供的chinamap包,里面包含了全国各省份的地理信息数据,get_map_china( )一句命令直接获取,没有安装的可以从Github上安装,但是这个包也存在一些不足,里面没有南海九段线的数据,故通过第二种方式从外部途径获取南海的地理信息数据。
remotes::install_github(“GuangchuangYu/chinamap”)
绘制静态地图
###绘制静态地图### cn<-get_map_china() factor(mydata$nvalue,levels = c("0","1~10","10~50","50~100","100~200","200~500","500~1000","1000~2000","2000~300000","1000000~"))->mydata$category##绘图前因子变量需要对其level进行重排,以便地图数据成功对应## str_sub(cn$province,1,2)->cn$id_area mapdata<-left_join(cn,mydata,by=c("id_area"="id_area")) ##自定义背景元素清除函数## theme_white<-function(){ theme(axis.line=element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background = element_blank()) } colo<-colorRampPalette(brewer.pal(9,"Reds"))(10)##调配填充颜色,由于数据变量是分类型,scale_fill_brewer(Palette="")中提供的颜色种类最多9种,故需要自行调配## source("lines_nanhai.R")##纳入南海九段线数据## read.csv("省会坐标.csv")->city##纳入省会标记## ##绘制静态地图## ggplot()+geom_polygon(data=mapdata,aes(x=long,y=lat,fill=category,group=group),colour="black",size=0.25)+scale_fill_manual(values = colo)+theme(legend.key.size=unit(0.3,"cm"),legend.position=c(0.2,0.2),legend.key.width=unit(0.3,"cm"),legend.title = element_text(size=8,hjust = 0))+theme_white()+labs(fill="现存确诊人数")+geom_sf(data=lines_nanhai,aes(geometry=geometry))+geom_text(data=city,aes(x=lon,y=lat,label=province),size=2)
结果如下图
如图所示,ggplot绘制的图形可视化是非常nice的,但是发现图例中存在NA,经检查发现原来是chinamap包中并没有提供澳门地区的地理信息数据。前两种方式绘制的图形通常不具备交互能力,无非进行拖拽,实际应用中比较受限,故考虑第三种方式绘图。3. 应用第三种方式绘制交互地图
本来想利用Remap包调用百度地图,一方面是需要事先申请注册地图认证,获取API秘钥AK;另一方面搭配的系统总是会存在崩溃状况。偶然发现李誉辉的一篇关于leaflet包的介绍,其在图形交互性能上能够得到有效实现。
绘制交互地图
###绘制疫情交互式地图### ##1.绘制确诊病例分布热图## regionNames("china")%>%substr(1,2)%>%data.frame()%>%setNames("id_area")->pr mydata<-left_join(pr,mydata,by=c("id_area"="id_area")) dat<-cbind(data.frame(regionNames("china")),mydata[,3:4])%>%setNames(c("area","confirm","category")) factor(dat$category,levels = c("0","1~10","10~50","50~100","100~200","200~500","500~1000","1000~2000","2000~300000","1000000~"))->dat$category##因子变量需要对其level进行重排,以便地图数据成功对应## map<-leafletGeo("china",dat)##将数据与地图数据关联## pal<-colorFactor(palette = c("Reds"),dat$category)##调取地图填充颜色fillcolor(),由于数据范围差异较大,不宜用连续性变量[colorNumeric]进行填充,故转化为分类变量填充## ##依据数据填充色彩,添加popup参数标签及鼠标显示,添加图例## casemap<-leaflet(map)%>%amap()%>%addPolygons(stroke = T,smoothFactor = 1,fillOpacity = 0.9,opacity = 0.8,weight = 1,color="gray",fillColor = ~pal(dat$category),popup=~htmltools::htmlEscape(value),popupOptions = popupOptions(closeButton = F,minWidth = 20),highlightOptions = highlightOptions(bringToFront = T,color = "blue",weight = 1))%>%addLegend(pal=pal,values = dat$category,position="bottomright",title = "现存本土确诊病例",opacity=1) ##2.添加风险区域标记## ##[1]纳入风险区域经纬度## read.csv("风险地区.csv",header=T)->fxarea str_sub(fxarea$name,1,3)->fxarea$area name<-as.character() for(i in 1:length(fxarea$name)){ name[i]<-str_sub(fxarea$name[i],4,str_length(fxarea$name[i])) } name->fxarea$name content<-data.frame() ##[2]将网页链接转化为popup参数提示框中能够识别的HTML语言## for(i in 1:length(fxarea$link)){ content[i,1]<-paste0("<br/>", "<b><a href='",fxarea$link[i],"'>",fxarea$area[i],"</a></b>","<br/>",fxarea$name[i]) } content$V1->fxarea$content ##自定义标记样式## mask<-makeIcon(iconUrl = "F://data//MASS//train//mask-100.png",iconWidth = 22,iconHeight = 32,iconAnchorX = 12,iconAnchorY =22)##标记可以来源于网络也可以来源于路径下存储的图片## ###绘制标记地图(添加popup参数提示框以及label参数标记标签### casemap%>%addTiles()%>%addMarkers(lng = ~fxarea$lng,lat=~fxarea$lat,popup = ~fxarea$content,label=~fxarea$cate,labelOptions = labelOptions(noHide = F,textOnly = F,style = list("color"="red","font-family"="serif","border-color"="blue")),icon=mask,popupOptions=popupOptions(closeButton = F,autoPan = F))
绘图中几个常用函数说明:
leaflet包
1.leaflet(dat)-----调用地图信息数据,可以是包含经纬度(lon,lat)的数据框,也可以是源于sp包传递的地理信息数据,最常用的是调用leafletGeo()合成的地图数据
2.addPolygons( )----主要参数:- stroke = T 显示绘制地图中多边形阴影,
- smoothFactor = 1 指定多边形边线平滑水平,
- fillOpacity = 0.9,opacity = 0.8 指定多边形透明度,
- weight = 1 指定多变形像素,
- color=“gray” 指定多边形边框颜色 ,
- fillColor =~pal(category) 指定数框中的变量填充多边形颜色(pal调色),
- popup=~htmltools::htmlEscape(value) html语言指定多边形中文本框显示内容,
- popupOptions = popupOptions(closeButton = F,minWidth = 20) 调解文本框位置及大小,
- highlightOptions = highlightOptions(bringToFront = T,color = “blue”,weight = 1))调解鼠标移动到多边形处显示情况。
3.addLegend(pal=pal 指定图例颜色,values = category 指定图例与数据框变量对应, position=“bottomright”,title = “现存本土确诊病例”,opacity=1 图例透明度)
4.addTiles()-----添加地理位置标记
5.addMarkers()----主要参数:- lng =~ lng,lat=~lat,需要标记位置的数据框中经纬度,
- popup=~content,popupOptions=popupOptions(closeButton = F,autoPan = F))指定标记位置文本框中显示内容,
- label=~fxarea$cate,labelOptions = labelOptions(noHide = F,textOnly = F,style = list(“color”=“red”,“font-family”=“serif”,“border-color”=“blue”)),指定标记显示内容, 调解鼠标移动到标记位置时显示情况
- mask<-makeIcon(iconUrl = “F://data//MASS//train//mask-100.png”,iconWidth = 22,iconHeight = 32,iconAnchorX = 12,iconAnchorY =22) 自定义标记样式
- icon=mask,改变标记样式
leafletCN包
是基于leaflet包做的一个大中华扩展,可以获取国内省市县的地理信息数据,在绘制国内地图时常用(举例):
1.regionNames(“china”)---------返回国内各省份所在名称
2.demomap(“北京市”)--------绘制北京市地图
3.geojsonMap(beijing_data,“北京市”)------依据数据框中数据为北京市地图填充颜色
4.leafletGeo(“china”,data_china)------将中国地图与数据框结合在一起,以便于leafet调用
5.amap()----调用高德地图图层结果如下图
更多相关内容 -
R语言:数据地图
2020-02-28 14:51:15R语言:绘制新冠肺炎数据地图世界地图的绘制中国地图的绘制方法一:方法二: 数据地图是一种经典的图示方法,因此 世界地图的绘制 R包中存储着常见地图的数据,比如maps包中包含世界地图、美国地图等 library(maps) ...数据地图是一种经典的图示方法,在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语言可复现
2022-04-16 22:21:39将来自另一个信息源的元数据添加到地图中。 有许多可以用于绘制地图的工具,一些流行的 GIS 软件允许点击交互进行地图开发和空间分析。这些工具具有不需要学习代码以及在地图上手动选择和放置图标和功能的便利性等...空间数据可视化&地图绘制&R语言可复现
绘制地理空间数据是一项常见的可视化任务,需要专门的工具。通常,问题可以分解为两个问题:
- 使用一个数据源绘制地图
- 将来自另一个信息源的元数据添加到地图中。
有许多可以用于绘制地图的工具,一些流行的 GIS 软件允许点击交互进行地图开发和空间分析。这些工具具有不需要学习代码以及在地图上手动选择和放置图标和功能的便利性等优点,例如:
ArcGIS - 由 ESRI 公司开发的商业 GIS 软件,非常流行但相当昂贵
QGIS - 一个免费的开源 GIS 软件,几乎可以做 ArcGIS 可以做的任何事情
使用 R 作为 GIS 一开始似乎更令人生畏,因为它不是“点击”,而是有一个“命令行界面”(您必须编写代码才能获得所需的结果)。但是,如果您需要重复生成地图或创建可重现的分析,这是一个主要优势。
本文主要解释使用
R
语言绘制地图,主要用到的包有:- rio : 导入数据
- tidyverse : 清洗、处理、绘图(包含ggplot2包)
- sf : 使用简单要素格式管理空间数据
- tmap : 生成地图,适用于交互式和静态地图
- OpenStreetMap : 在ggplot中免费添加 OSM 国内外地图作为底图
本目文录:
- 综述
- 空间数据类型 Spatial data
- 矢量数据 Vector data
- 栅格数据 Raster data
- 空间数据基本格式
- 基本地图类型
- 绘制地图准备工作
- 必备环节
- 所需R包
- 实例案例数据
- 行政边界数据
- 人口数据
- 卫生设施数据上
- 画坐标系
- 空间数据连接类型
- 点数据与多边形区域连接
- 最近的邻域连接
- 缓冲区连接
- 其他空间连接类型
- 等值域图 Choropleth 地图
- 使用ggplot绘制地图
- 底图 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:
- 导入 shapefile
- 清理列名
- 过滤行以仅保留感兴趣的区域
要导入
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 的语法类似,详情点击。- 设置 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)相交:
- 从行列表(点)开始
- 空间连接到边界,将连接类型
join
设置为“st_intersects” - 使用 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
连接方法来可视化离个别病例最近的医疗机构。- 我们从
shapefile
linelist
linelist_sf
开始 - 我们在空间上加入
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 文章中查看有关地图投影和坐标系的更多信息。该博客讨论了不同类型的地图投影,以及如何根据感兴趣的区域和地图/分析的背景选择合适的投影。
- 首先,在每个卫生设施周围创建一个半径约为 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')
- 其次,我们使用
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 开始,它是每个行政单位及其案例数量的汇总表。
- 人口数据 sle_adm3_pop 使用 dplyr 的 left_join() 连接,基于 case_adm3 数据帧中的 admin3pco 列和 sle_adm3_pop 数据帧中的 adm_pcode 列的共同值。请参阅加入数据页面)。
- select() 应用于新的数据框,只保留有用的列 - total 是总人口
- 每 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
下面我们描述了如何在底图上实现案例的等高线密度热图,从一个线列表(每个案例一行)开始。
- 从 OpenStreetMap 创建底图切片,如上所述
- 使用 latitude 和 longitude 列从 linelist 中绘制案例
- 使用 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
-
地图R语言
2022-05-26 08:19:41leaflet包 leaflet包相对其它地图包,有很多优点和缺点, 首先,绘制地图简单快捷,因为都是基于供应商的tiles,一行代码就可以render出基本widget地图。 支持管道传参,一个图层...对于在地图上添加markers图标,shapleaflet包
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
-
【R语言作图】如何在地图上任意位置画饼图直方图等
2021-01-13 13:07:52N久不分享学习心得的小Q终于有时间冒个泡了,这次分享一个小Q...2.已有的map数据中一个国家对应一个坐标,一个国家边界,利用这些已有数据+用户数据构建新的画图数据(其他新添加的图均是如此)。rworldmap包的介绍... -
如何在地图上添加柱形图(R语言)
2021-01-17 13:15:17(百度百科) #加载必要的包 library(mapproj) library(ggplot2) library(rgdal) library(plyr) #注意:由于重庆的行政编吗为500000,在R语言内会自动显示为5e5,后面利用行政编码进行数据联结会出错 options(scipen... -
R语言常用的数据可视化包总结
2022-04-08 15:24:19本文中的内容来自书籍《R语言数据可视化实战 (微视频全解版) ——大数据图表从入门到精通》。书中有更多的R语言数据可视化分析使用案例,可学习更多的关于R语言数据可视化的内容。 数据可视化,是关于数据视觉... -
绘图杂记【5】R语言绘制热力地图(重庆热力地图)
2020-12-13 11:42:28text(family="myFont",size=19,hjust = 0.5), ) # 主题设置https://blog.csdn.net/u014801157/article/details/24372531 关于颜色 参考:R语言的颜色应用 颜色数据集RColorBrewer library(RColorBrewer) brewer.pal.... -
R语言-地图绘制的思路
2020-02-26 16:49:47R中的画地图的方法不外乎两种,一种是利用包里GIS方面的数据,在R中直接画出来,另一种是从其他地方拿到数据,在R中通过某些包解析后再展现成。这两种方法只是数据来源不同,具体的画图以及美化方法无异。 第一种,... -
R语言做地理分析:数据导入,合并,绘图,导出
2021-09-28 14:02:38学习笔记,原文链接????:... 0 数据准备 0.1 伦敦边界数据 (.shp) 0.2 伦敦2008-2020垃圾乱丢事件统计数据 (.csv) https://pan.baidu.com/s/1mBri -
R 语言数据分析/数据挖掘常用包
2018-07-30 21:31:33数据分析的流程主要包括以下几个4大步骤,而熟练掌握以下每个步骤的常用包,就能提高数据科学的效率与质量。 1.数据读取 2.数据清洗 3.数据可视化 4.数据挖掘 细分为以下流程: 1.1 数据导入 1.2 数据整理... -
送书|北大出版:R语言数据分析与可视化从入门到精通
2020-10-23 15:10:00生物信息学习的正确姿势NGS系列文章包括NGS基础、高颜值在线绘图和分析、转录组分析(Nature重磅综述|关于RNA-seq你想知道的全在这)、ChIP-seq分析(ChIP-se... -
R语言在线地图神器:Leaflet for R包(一)
2016-10-14 16:14:04做Javascript相关地图开发的码农,特别关心可视化和开源的同学,都听说过Leaflet这样一个神包(神马,你没有听说过……好吧,当我没说,你自己先搜索一下……) 用官方(自吹自擂)的话来说,Leaflet包是号称最受... -
利用R语言如何画出广州房价地图
2015-05-03 01:37:39R软件的ggplot包升级以后有绘制地图的新功能,其图形元素主要是通过geom_map来实现。由于系统内maps包所自带的地图数据没有广州市的数据。其它图家地图数据则要从外部导入,本文则尝试从外部导入广州数据,然后用... -
利用R语言画出两地路线图
2017-07-25 10:46:50发一篇教大家如何用R语言画路线图的攻略~ ~\(≧▽≦)/~ 题目是这样的: 利用R的地图包,提取武汉地图,并在图上标注一条从华中师范大学北门到黄鹤楼的路线. -
【R 数据科学】R语言进行数据科学整理最有用的包大全
2017-06-09 11:05:12一、数据科学工作流程1.1 数据导入 1.2 数据整理 1.3 反复理解数据 1.4 数据可视化 1.5 数据转换 1.6 统计建模 1.7 作出推断(比如预测) 1.8 沟通交流 ...在R和python上都可使用 readr:实现 -
R语言数据清理:视频游戏数据案例研究
2019-06-12 14:08:03在本教程中,您将学习如何将中等大小的数据集(如游戏元数据)转换为有用的格式,以便使用R进行进一步分析。 您将了解整洁数据集遵循的关键原则,为什么跟踪它们有用,以及如何清理您给出的数据。整理也是了解新... -
《R语言与数据挖掘》⑤高级绘图工具【lattice包】【ggplot2】【交互式】
2021-12-20 16:36:20书籍:《R语言与数据挖掘》 作者:张良均 出版社:机械工业出版社 ISBN:9787111540526 本书由北京华章图文信息有限公司授权杭州云悦读网络有限公司电子版制作与发行 版权所有·侵权必究 lattice包 lattice包的图形... -
技巧 | 在R语言中使用高德地图的API进行地理/逆地理编码(地址与经纬度的相互转换)...
2022-03-07 00:54:06高德地图和百度地图都提供了坐标拾取系统,通过坐标查询或坐标反查操作可以查询一个地址对应的经纬度或经纬度对应的地址名称。但是,手动查询的方式效率很低,也不能进行批量查询。本篇就来介绍在R语言... -
R数据分析:方法与案例详解--自学笔记
2021-11-18 20:06:36第二章 数据结构与基本运算 2.1 数据类型 数值型(numeric) 整数 小数 科学数 字符型(character) == 夹杂单引号或者双引号之间==“MR” 逻辑型 ==只能读取T (TRUE)或 F (FALSE)值 复数型 a+bi 原始型(raw) 以... -
R语言与地图(一)
2020-12-18 21:34:19原标题:R语言与地图(一)Hello,大家好! 很高兴能借助这个平台每周5分钟与大家分享 一点儿数据分析 的那些事。数据分析的工具很多,而R语言因其免费开源、易于解释和统计分析很受高校和业界的青睐。今天介绍 R软件在... -
这个报表工具绝了!能做GIS数据地图,还能集成R语言!
2019-11-04 15:47:28最近更是新出了几款插件,能在线预览文档,能做GIS数据地图,还能集成R语言! 这还是我认识的报表工具? 没错,就是FineReport ! 一起来看看这几个锦上添花的功能吧! 1、集成R语言 在数据分析、数据统计和... -
利用R语言绘制世界航班路线图
2018-04-17 12:53:29taoyan:伪码农,R语言爱好者,爱开源。 个人博客: https://ytlogos.github.io/ 、 **一、 简介** 本文基于NASA的夜间地图... -
空间数据分析与R语言实践
2019-11-21 12:49:37T贝尔实验室开发的一种用来进行数据探索、统计分析、作图的解释型语言。它的丰富的数据类型(向量、数组、列表、对象等)特别有利于实现新的统计算法,其交互式运行方式及强大的图形及交互图形功能使得我们可以方便... -
R语言入门(1)-初识R语言
2021-05-21 18:08:57首先在用户根目录下cat查看一下, 发现没有.Renviron文件, 这个是R语言的环境配置文件.2. 那么就用echo语句追加一句"LANGUAGE=en" 到.Renviron文件, 如果没有这个文件, echo语句会自动创建.3. 然后再cat查看一下, ... -
揭秘!文字识别在高德地图数据生产中的演进
2020-07-30 14:45:57本文分享文字识别技术在高德地图数据生产中的演进与实践,介绍了文字识别自研算法的主要发展历程和框架,以及未来的发展和挑战。一 背景作为一个DAU过亿的国民级软件,高德地图每天为用户提供海量的 -
R语言入门
2021-05-22 06:09:45相关学习网站R语言常用领域探索性数据分析统计推断回归分析1.线性模型拟合数据。(预测变量、结果变量),通过预测变量来预测结果变量。2.预测:利用建立好的线性模型预测未来情况机器学习-分类问题训练模型+预测开发... -
带南海九段线分位数地图可视化(R语言版)
2019-11-13 01:27:41内容是使用R语言进行带南海九段线分位数地图可视化。虾神的原博文地址如下(Python版)。 Python实现带南海九段线分位数地图完整可视化版本(附代码及数据) 1999-2017年中国各省旅游外汇收入分析及可视化(附代码及... -
技巧 | 如何使用R语言的基础绘图系统的拼图功能
2021-05-30 01:09:04我们知道ggplot2工具包有很多方便的拼图拓展包,如cowplot、patchwork等,而本篇就来介绍在使用R语言的基础绘图系统如何进行拼图。需要明确的是,基础绘图系统的拼图功能不需要...