用R绘制世界地图及中国地图

在如今这个数据可视化无处不在的世界里，R语言不再只是统计分析的工具，它更像是一支画笔，能为你描绘出一幅幅地图。如果你曾经幻想过用代码勾勒出世界的轮廓，或者展现中国的山川河流，那么今天的探索正适合你~

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.2.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require(maps)

## Loading required package: maps

## Warning: package 'maps' was built under R version 4.2.3

require(viridis)

## Loading required package: viridis

## Warning: package 'viridis' was built under R version 4.2.3

## Loading required package: viridisLite

## Warning: package 'viridisLite' was built under R version 4.2.3

## 
## Attaching package: 'viridis'

## The following object is masked from 'package:maps':
## 
##     unemp

1.世界地图绘制

绘制简单地图

world <- map_data("world")
worldplot <- ggplot() +
  geom_polygon(data = world, aes(x=long, y = lat, group = group),fill="lightgray",color="white") + 
  coord_fixed(1.3)
worldplot

 根据指标填充颜色，如想绘制全球妊娠期糖尿病发病情况的分布比较图，可下载IDF文件。

加载数据

world_gdm <- read.csv("IDF_2021.csv", header = TRUE) #下载文件到本地文件夹

head(world_gdm)

##   X.                       Region            Country.Territory X2000 X2011
## 1  1                       Africa                       Africa   N/A   N/A
## 2  2                       Europe                       Europe   N/A   N/A
## 3  3 Middle East and North Africa Middle East and North Africa   N/A   N/A
## 4  4  North America and Caribbean  North America and Caribbean   N/A   N/A
## 5  5    South and Central America    South and Central America   N/A   N/A
## 6  6              South-East Asia              South-East Asia   N/A   N/A
##   X2021 X2030 X2045   Type                      Country
## 1    13   N/A   N/A Region                       Africa
## 2    15   N/A   N/A Region                       Europe
## 3  14.1   N/A   N/A Region Middle East and North Africa
## 4  20.7   N/A   N/A Region  North America and Caribbean
## 5  15.8   N/A   N/A Region    South and Central America
## 6  25.9   N/A   N/A Region              South-East Asia

##查看两个数据集中变量名称不一致的部分
diff <- setdiff(world$region, world_gdm$Country)
diff

##  [1] "Antarctica"                          "French Southern and Antarctic Lands"
##  [3] "Barbuda"                             "Saint Barthelemy"                   
##  [5] "Brunei"                              "Ivory Coast"                        
##  [7] "Cape Verde"                          "Czech Republic"                     
##  [9] "Canary Islands"                      "Falkland Islands"                   
## [11] "Guernsey"                            "Heard Island"                       
## [13] "Cocos Islands"                       "Christmas Island"                   
## [15] "Chagos Archipelago"                  "Jersey"                             
## [17] "Siachen Glacier"                     "Nevis"                              
## [19] "Kosovo"                              "Saint Martin"                       
## [21] "Montserrat"                          "Norfolk Island"                     
## [23] "Bonaire"                             "Sint Eustatius"                     
## [25] "Saba"                                "Pitcairn Islands"                   
## [27] "Madeira Islands"                     "Azores"                             
## [29] "Western Sahara"                      "South Sandwich Islands"             
## [31] "South Georgia"                       "Saint Helena"                       
## [33] "Ascension Island"                    "Saint Pierre and Miquelon"          
## [35] "Swaziland"                           "Sint Maarten"                       
## [37] "Turks and Caicos Islands"            "Tobago"                             
## [39] "Vatican"                             "Grenadines"                         
## [41] "Wallis and Futuna"

通常需要对不一致的变量进行重命名，以使得world map数据和GDM发病数据中的国家名称一致，这里仅用作示例，不做修改。

## 将数据设置为数值型变量
world_gdm$X2021 <- as.numeric(as.character(world_gdm$X2021))

## Warning: NAs introduced by coercion

world_gdm$prevalence<-world_gdm$X2021
world_gdm$prevalence<-gsub("-","",world_gdm$prevalence)

合并数据集

worldSubset <- inner_join(world, world_gdm, by = c("region"="Country"))
head(worldSubset)

##        long      lat group order region subregion X.
## 1 -69.89912 12.45200     1     1  Aruba      <NA> 18
## 2 -69.89571 12.42300     1     2  Aruba      <NA> 18
## 3 -69.94219 12.43853     1     3  Aruba      <NA> 18
## 4 -70.00415 12.50049     1     4  Aruba      <NA> 18
## 5 -70.06612 12.54697     1     5  Aruba      <NA> 18
## 6 -70.05088 12.59707     1     6  Aruba      <NA> 18
##                        Region Country.Territory X2000 X2011 X2021 X2030 X2045
## 1 North America and Caribbean             Aruba     -     -    NA     -     -
## 2 North America and Caribbean             Aruba     -     -    NA     -     -
## 3 North America and Caribbean             Aruba     -     -    NA     -     -
## 4 North America and Caribbean             Aruba     -     -    NA     -     -
## 5 North America and Caribbean             Aruba     -     -    NA     -     -
## 6 North America and Caribbean             Aruba     -     -    NA     -     -
##      Type prevalence
## 1 Country       <NA>
## 2 Country       <NA>
## 3 Country       <NA>
## 4 Country       <NA>
## 5 Country       <NA>
## 6 Country       <NA>

绘制最终数据

## 设置绘图主题
plain <- theme(
  axis.text = element_blank(),
  axis.line = element_blank(),
  axis.ticks = element_blank(),
  panel.border = element_blank(),
  panel.grid = element_blank(),
  axis.title = element_blank(),
  panel.background = element_rect(fill = "white"),
  plot.title = element_text(hjust = 0.5)
)

worldSubset$prevalence<-as.numeric(worldSubset$prevalence)

#绘图
worldHDI <- ggplot(data = worldSubset, mapping = aes(x = long, y = lat, group = group)) + 
  coord_fixed(1.3) +
  geom_polygon(aes(fill = prevalence)) +
 # scale_fill_distiller(palette ="RdBu", direction = -1) + # or direction=1
  scale_fill_gradient2(low = "#ADC6AD", mid = "#e7b800", high = "red", midpoint = 10)+
  ggtitle("Prevalence of Gestational Diabetes in 2021 by WHO region") +
  plain

worldHDI

2.绘制特定区域的地图

获取特定国家的地图数据

# Some EU Contries
some.eu.countries <- c(
  "Portugal", "Spain", "France", "Switzerland", "Germany",
  "Austria", "Belgium", "UK", "Netherlands",
  "Denmark", "Poland", "Italy", 
  "Croatia", "Slovenia", "Hungary", "Slovakia",
  "Czech republic"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)
# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  dplyr::group_by(region) %>%
  dplyr::summarise(long = mean(long), lat = mean(lat))

region.lab.data

## # A tibble: 17 × 3
##    region          long   lat
##    <chr>          <dbl> <dbl>
##  1 Austria        13.5   47.6
##  2 Belgium         4.73  50.6
##  3 Croatia        16.3   44.6
##  4 Czech Republic 15.4   49.9
##  5 Denmark        10.7   55.7
##  6 France          3.23  46.2
##  7 Germany        10.4   51.2
##  8 Hungary        19.4   47.2
##  9 Italy          11.8   42.2
## 10 Netherlands     5.32  52.1
## 11 Poland         19.1   51.4
## 12 Portugal       -7.89  39.9
## 13 Slovakia       19.6   48.8
## 14 Slovenia       14.8   46.1
## 15 Spain          -2.91  40.7
## 16 Switzerland     8.31  46.7
## 17 UK             -4.10  55.6

可视化

ggplot(some.eu.maps, aes(x = long, y = lat)) +
  geom_polygon(aes( group = group, fill = region))+
  geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  scale_fill_viridis_d()+
  theme_void()+
  theme(legend.position = "none")

3.中国地图绘制

下载中国地图json文件

DataV.GeoAtlas地理小工具系列

读取地图文件，绘制基础图形

library(sf)

## Warning: package 'sf' was built under R version 4.2.3

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(ggplot2)
library(cowplot)

## Warning: package 'cowplot' was built under R version 4.2.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'purrr' was built under R version 4.2.3

## Warning: package 'stringr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ purrr::map()       masks maps::map()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggspatial)

## Warning: package 'ggspatial' was built under R version 4.2.3

library(dplyr)
China_map=read_sf("China.json")
# 地图散点图2
p2=ggplot(China_map)+
  geom_sf(color='white',fill="lightgray",size=0.8)+#地图线条粗细
  annotation_scale(location = "bl", width_hint = 0.3) +#添加比例尺并调整位置及长度
  annotation_north_arrow(location = "tl", which_north = F, 
                         pad_x = unit(0.05, "in"), pad_y = unit(0.05, "in"),
                         style = north_arrow_nautical)+#添加指北针，指北针类型style有north_arrow_orienteering；north_arrow_fancy_orienteering；north_arrow_minimal与north_arrow_nautical四种类型
  theme_map()
p2

## Scale on map varies by more than 10%, scale bar may be inaccurate

提取各个省份经纬度坐标

city_name<-China_map$name
location_center<-China_map$center

data_province<-matrix(nrow = 34,ncol = 2)
colnames(data_province)<-c("LON","LAT")
for (i in 1:34) {
  data_province[i,1]<-location_center[[i]][1]
  data_province[i,2]<-location_center[[i]][2]
 # rownames(data_province[i,])<-city_name[i]
}
rownames(data_province)<-city_name[1:34]
#write.csv(data_province,file = "Data/Map/data_province_center.csv")
head(data_province)

##                   LON      LAT
## 北京市       116.4053 39.90499
## 天津市       117.1902 39.12560
## 河北省       114.5025 38.04547
## 山西省       112.5492 37.85701
## 内蒙古自治区 111.6708 40.81831
## 辽宁省       123.4291 41.79677

更简单点的，使用mapdata包中的map()函数，即可绘制中国地图。

获取经纬度坐标之后，就可以根据自己的数据填充感兴趣的数值了，具体方法同前，感兴趣的大家可以尝试一下。

希望可以通过今天的博文，大家可以稍微获得一点代码与地理的浪漫邂逅。

欲探索更多R包使用方法，请关注微信公众号《单细胞学会》。