前言
![]()
这张图来自1000 Genomes Project Consortium. A global reference for human genetic variation[J]. Nature, 2015, 526(7571): 68-74.
这张地理分布图去除了南北极，以及国家之间的分界线，再在地图图层上加了饼图展示地理分布和分类信息。我决定我也画一个。。。
我使用过maps包等，画出来的图都觉得不够好看，这里用ggplot(主要是geom_polygon函数)用罗宾逊投影(Robinson projection) 坐标系画世界地图，当然，这里参考了网上很多文章，链接我都放在文章末尾了。
地图数据
地球上的大陆在图上是一个一个多边形，我们这里主要用geom_polygon多边形绘图函数完成工作，首先得下载地图数据。
下载后解压到工作目录下的文件夹中，如
- /workdir/ne_110m_graticules_all
- /workdir/ne_110m_land
R 包准备
我们需要rgdal包来读取地图数据文件。
1
2
3
4
| install.packages("rgdal")
install.packages("ggplot2")
install.packages("reshape")
install.packages("scatterpie")
|
画图
世界地图（罗宾逊投影）
文章的地图中，大海的颜色是white，陆地的颜色是 #B4D9A2,陆地的border是 #0F8538。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
library(rgdal)
library(ggplot2)
library(reshape)
library(scatterpie)
setwd("/workdir/")
wmap <- readOGR(dsn="ne_110m_land", layer="ne_110m_land")
wmap_df <- fortify(wmap)
wmap_robin <- spTransform(wmap, CRS("+proj=robin"))
wmap_df_robin <- fortify(wmap_robin)
theme_opts <- list(theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_rect(fill="white"),
panel.border = element_blank(),
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=22)))
ggplot(wmap_df_robin, aes(long,lat, group=group,fill=hole)) +
geom_polygon(colour="#0F8538",lwd=0.05) +
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#B4D9A2", "white"), guide="none")
ggsave("test.pdf", width=16, height=9)
|
效果图如下
![]()
世界地图（罗宾逊投影）+ 边框
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
bbox <- readOGR("ne_110m_graticules_all", layer="ne_110m_wgs84_bounding_box")
bbox_df<- fortify(bbox)
bbox_robin <- spTransform(bbox, CRS("+proj=robin"))
bbox_robin_df <- fortify(bbox_robin)
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat, group=group, fill=hole),colour="#0F8538",lwd=0.05) +
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#B4D9A2", "white"), guide="none")
ggsave("test1.pdf", width=16, height=9)
|
效果图如下
![]()
这个图中还可以添加经纬线，以及国界线。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
countries <- readOGR("ne_110m_admin_0_countries", layer="ne_110m_admin_0_countries")
countries_robin <- spTransform(countries, CRS("+init=ESRI:54030"))
countries_robin_df <- fortify(countries_robin)
bbox <- readOGR("ne_110m_graticules_all", layer="ne_110m_wgs84_bounding_box")
bbox_df<- fortify(bbox)
bbox_robin <- spTransform(bbox, CRS("+proj=robin"))
bbox_robin_df <- fortify(bbox_robin)
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat, group=group, fill=hole),colour="#0F8538",lwd=0.05) +
geom_polygon(data=countries_robin_df, aes(long,lat, group=group, fill=hole)) +
geom_path(data=countries_robin_df, aes(long,lat, group=group, fill=hole), color="black", size=0.3) +
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#B4D9A2", "white"), guide="none")
|
世界地图（罗宾逊投影）+ 边框 + 点图
背景画好了，我们得开始展现我们的样本的信息，比如样本的地理分布，该地区取样的多少用圆点的大小来衡量。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
mydata = read.csv("jingweidu_34.txt",head=T)
head(mydata)
places_robin_df <- project(cbind(mydata$jingdu, mydata$weidu),
proj="+init=ESRI:54030")
places_robin_df <- as.data.frame(places_robin_df)
names(places_robin_df) <- c("LONGITUDE", "LATITUDE")
places_robin_df$num <- mydata$geshu
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat,
group=group,
fill=hole),
colour="#0F8538",
lwd=0.05) +
geom_point(data=places_robin_df, aes(LONGITUDE,
LATITUDE,
group=NULL,
fill=NULL,
size=num),
color="#808080",
alpha=I(8/10))+
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#B4D9A2", "white"), guide="none") +
scale_size_continuous(range=c(5,20), guide="none")
ggsave("test2.pdf", width=16, height=9)
|
效果图如下
![]()
世界地图（罗宾逊投影）+ 边框 + 分组点图
点图都画出来，分组点图也就简单了，在输入数据中加入分组信息(group)。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
| mydata= read.csv("jingweidu_34_g.txt",head=T)
head(mydata)
places_robin_df <- project(cbind(mydata$jingdu, mydata$weidu),
proj="+init=ESRI:54030")
places_robin_df <- as.data.frame(places_robin_df)
names(places_robin_df) <- c("LONGITUDE", "LATITUDE")
places_robin_df$num <- mydata$geshu
places_robin_df$group <- mydata$group
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat,
group=group,
fill=hole),
colour="#0F8538",
lwd=0.05) +
geom_point(data=places_robin_df, aes(LONGITUDE,
LATITUDE,
group=NULL,
fill=NULL,
size=num,
color=group),
alpha=I(8/10))+
scale_color_manual(values=c("#808080","#D2B48C")) +
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#B4D9A2", "white"), guide="none") +
scale_size_continuous(range=c(5,20), guide="none")
ggsave("test3.pdf", width=16, height=9)
|
效果图如下
![]()
世界地图（罗宾逊投影）+ 边框 + 饼图
首先整理你的样本信息为下图样式（当然，这不是必须的，也可以直接整理成画图的格式）。
![]()
每一行代表一个样本，第三列代表分组
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
| mydata <- read.table("jingweidu_34_scatterpie.txt", head=T)
head(mydata)
mydata_reshape<- cast(mydata,Latitude+Longitude~Sub.population)
mydata_reshape$radius=(((mydata_reshape$Cluster1+
mydata_reshape$Cluster2+
mydata_reshape$Cluster3)%/%25)+700000)
head(mydata_reshape)
places_robin_df <- project(cbind(mydata_reshape$Longitude,
mydata_reshape$Latitude),
proj="+init=ESRI:54030")
places_robin_df <- as.data.frame(places_robin_df)
names(places_robin_df) <- c("LONGITUDE", "LATITUDE")
places_robin_df$Cluster1 <- mydata_reshape$Cluster1
places_robin_df$Cluster2 <- mydata_reshape$Cluster2
places_robin_df$Cluster3 <- mydata_reshape$Cluster3
places_robin_df$radius <- mydata_reshape$radius
head(places_robin_df)
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat,
group=group,
fill=hole),
colour="#0F8538",
lwd=0.05) +
geom_scatterpie(aes(x=LONGITUDE, y=LATITUDE,r=radius),
data=places_robin_df,
cols=c("Cluster1","Cluster2","Cluster3"),
color="black",
alpha=1,lwd=0)+
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#DCDCDC","#808080","#FFFF00","#B4D9A2","white"),
guide="none")
ggsave("test4.pdf", width=16, height=9)
|
效果图如下
![]()
只显示部分地图
直接对画图用的bbox_robin_df wmap_df_robin 两个数据框进行过滤就行，但是我用的这个地图的数据，两个点的距离有点远，所以边界显示效果不是很好。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|
bbox_robin_df<-bbox_robin_df[which(bbox_robin_df$lat >= -2.929533e+06),]
wmap_df_robin <- wmap_df_robin[which(wmap_df_robin$lat >= -2887685),]
places_robin_df <- places_robin_df[which(places_robin_df$LATITUDE >= -2887684),]
ggplot(bbox_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="white",colour="black",lwd=0.05) +
geom_polygon(data=wmap_df_robin, aes(long,lat,
group=group,
fill=hole),
colour="#0F8538",
lwd=0.05) +
geom_scatterpie(aes(x=LONGITUDE, y=LATITUDE,r=radius),
data=places_robin_df,
cols=c("Cluster1","Cluster2","Cluster3"),
color="black",
alpha=1,lwd=0)+
coord_equal() +
theme_opts +
scale_fill_manual(values=c("#DCDCDC","#808080","#FFFF00","#B4D9A2","white"),
guide="none")
ggsave("test5.pdf", width=16, height=9)
|
效果图如下
![]()
追加更新：改变地图中心
1
2
3
4
|
wmap_robin <- spTransform(wmap, CRS("+proj=robin"))
wmap_robin <- spTransform(wmap, CRS("+proj=robin +lon_0=150 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
|
总结
到这里我们模仿也就完成基本完成了，我觉得还有用是添加国界线，但是这里加上就少了简洁美了。
还可以优化的是，找个更好的地图数据，让边界更好看。
剩下的就是按照自己需求改改图了
参考
