R 语言：如何从地图图片上提取数据并重新绘图？

最近有个小伙伴遇到了一个这样的问题，他想使用一篇文献中的数据，不过这篇文献并没有提供具体数据，而仅仅是提供了一张地图图片，例如这样的（这个是我随便从一篇论文中找到的）：

那么遇到这种情况，我们该如何从图片中提取数据并重新绘制这幅图呢（毕竟论文是不允许使用截图的）。对于这个地图显然数据已经不太可能提取出来了（因为作者是对数据进行了分段），所以只能说看看能不能提取颜色，然后重新绘制。

在开始内容前，我们先设置下字体和绘图主题：

# 设置字体
library(showtext)
library(ggplot2)
library(tidyverse)
showtext_auto(enable = T)
font_add("cnfont", regular = "song.otf")
cnfont <- "cnfont"
theme_set(hrbrthemes::theme_ipsum(base_family = cnfont))

通常我们可以使用识色软件识别每个城市的颜色（这里的地图是城市地图），然后列个表。但是这很难实现，毕竟很难一一准确的识别 300 多个城市的颜色。这里我提供一种思路供大家参考。

这种方法仅限于我们能猜到这幅地图的坐标参考系的情况，例如上面的地图的坐标参考系是：“+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs”（这个很常用，自然资源部的很多地图都是用的这个）。

我的思路是先把这个图片作为栅格数据读取，然后再使用城市地理矢量数据提取每个区域的平均颜色，最后再把这些颜色绘制到地图上。

首先我们通过恰当的裁剪把中国地图的部分裁剪出来：

然后就可以把这个图片作为栅格数据读取了：

library(terra)
rast("pic1.png") -> rst
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
crs(rst) <- mycrs

这里我们虽然设定了 crs（坐标参考系），但是还没有设定范围，下面我们需要找到这个栅格数据的范围（extent）。

cnrange.geojson 文件是我使用 https://tidyfriday.cn/geojson 网站绘制的一个中国地图大致范围（不含南海）：

这里的主要目的就是剔出南海部分，避免对等下范围的提取产生影响：

library(sf)
read_sf("cnrange.geojson") -> cnrange
read_sf("2021行政区划/市.shp") %>%
  st_intersection(cnrange) -> city
city %>%
  st_transform(mycrs) -> city
# 简化下方便绘图展示
city %>%
  st_simplify(dTolerance = 2000) -> citysim

# 提取 citysim 的范围生成一个 sf 对象
st_as_sfc(st_bbox(citysim), crs = mycrs) -> rangesf

# city 的范围就是刚刚读取的栅格数据的正确范围：
st_bbox(city) -> citybbox

然后就可以用这个 citybbox 给 rst 设定范围了：

ext(rst) <- c(citybbox[1], citybbox[3], citybbox[2], citybbox[4])

然后我们绘制一幅图看看效果：

terra::plot(vect(rangesf), col = "red")
terra::plot(rst$pic1_1, add = T)
terra::plot(vect(citysim), add = T)

这样我们就可以分城市提取各个城市的平均颜色：

terra::extract(rst$pic1_1, vect(city), fun = "mean") -> citydf1
terra::extract(rst$pic1_2, vect(city), fun = "mean") -> citydf2
terra::extract(rst$pic1_3, vect(city), fun = "mean") -> citydf3

rst 实际上总共是 4 个 layer，其中前三个分别是 RGB 颜色的 R 值、G 值 和 B 值，然后我们再把提取得到的三个数据框合并起来组合成颜色：

library(tidyverse)
bind_cols(citydf1, citydf2, citydf3) %>%
  as_tibble() %>%
  select(ID = `ID...1`, contains("pic")) -> citydf

library(farver)
citydf %>%
  mutate_at(vars(contains("pic")), as.integer) -> citydf
citydf$hex <- encode_colour(citydf[,2:4])

然后我可以绘制一副简单的地图看看效果：

citydf %>%
  bind_cols(citysim) %>%
  st_sf() %>%
  ggplot() +
  geom_sf(aes(fill = I(hex)))

这样我们就把这幅图复现出来了，不过还差了图例。

图例可以通过提取颜色绘制出来。

大家可以使用这个图片识色软件：(讲义文档里面有)。

然后提取对应的颜色即可。

如果上面的应用过期了，可以试试这个：https://tidyfriday.cn/colorpicker/

# 图例
tribble(
  ~class, ~color,
  "No data", "#FFFFFF",
  "0 ~ 500", "#94D1EF",
  "500 ~ 1000", "#6FB2ED",
  "1000 ~ 2000", "#4092E7",
  "2000 ~ 3000", "#4092E7",
  "3000 ~ 4000", "#204FC1",
  "4000 ~ 5000", "#182DAB",
  "> 5000", "#0C0B93"
) -> colordf

# 随便绘制一幅图
tibble(x = 1:8, y = 1:8, z = colordf$class) %>%
  mutate(z = factor(z, levels = colordf$class)) %>%
  ggplot(aes(x, y)) +
  geom_col(aes(fill = z)) +
  scale_fill_manual(values = colordf$color) +
  labs(fill = "City level energy consumption\n(million tons)") +
  guides(fill = guide_legend(nrow = 4, byrow = F,
                             keywidth = unit(1, "cm"),
                             keyheight = unit(0.5, "cm"))) -> p

# 提取图例
library(cowplot)
get_legend(p) -> legend
ggdraw(legend)

然后我们就可以再次绘图了。这里我们仿照样图绘制带小地图版本的中国市级地图：

# 再次绘制主图
# 省份边界
read_sf("2021行政区划/省.shp") %>%
  st_cast("MULTILINESTRING") %>%
  mutate(name = "省级边界") %>%
  st_transform(mycrs) -> provline

# 海岸线 & 九段线
read_sf("海岸线/海岸线.shp") %>%
  st_transform(mycrs) %>%
  mutate(name = "海岸线") -> hax
read_sf("九段线.geojson") %>%
  st_transform(4326) %>%
  st_transform(mycrs) %>%
  mutate(name = "九段线") -> jdx

# 绘图数据
citydf %>%
  bind_cols(city) %>%
  select(-ID, -contains("_")) -> citydf

# 合并上面的文件
citydf %>%
  bind_rows(provline) %>%
  bind_rows(jdx) %>%
  bind_rows(hax) %>%
  mutate(id = row.names(.)) %>%
  st_sf() -> tempdf

# 拆分成主图和小地图并合并
# 主图
main_bbox <- st_bbox(c(xmin = -2625586,
                       xmax = 2206965,
                       ymax = 5921583,
                       ymin = 1836814.1),
                     crs = st_crs(mycrs)) %>%
  st_as_sfc()

tempdf %>%
  st_make_valid() %>%
  st_intersection(main_bbox) -> mainchina

# 南海九段线小地图
small_bbox <- st_bbox(c(xmin = 120000,
                        xmax = 1766004.1,
                        ymax = 2557786.0,
                        ymin = 320000),
                      crs = st_crs(mycrs)) %>%
  st_as_sfc()

nanhai <- tempdf %>%
  st_make_valid() %>%
  st_intersection(small_bbox) %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs)

# 合并
for(i in nanhai$id){
  print(i)
  mainchina[mainchina$id == i,]$geometry <- st_combine(
    c(mainchina[mainchina$id == i,]$geometry,
      nanhai[nanhai$id == i,]$geometry)
  )
}

mainchina

#> Simple feature collection with 406 features and 15 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -2625586 ymin: 1875139 xmax: 2206965 ymax: 5921583
#> CRS:           +proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> # A tibble: 406 × 16
#>    hex     省     省代码 省类型 市    市代码 市类型 name  CNAME ASCRIPTION GB
#>    <chr>   <chr>   <dbl> <chr>  <chr>  <dbl> <chr>  <chr> <chr> <chr>      <chr>
#>  1 #ADCAE5 安徽省 340000 省     安庆… 340800 地级市 <NA>  <NA>  <NA>       <NA>
#>  2 #AFCDE9 安徽省 340000 省     蚌埠… 340300 地级市 <NA>  <NA>  <NA>       <NA>
#>  3 #AFCCE9 安徽省 340000 省     亳州… 341600 地级市 <NA>  <NA>  <NA>       <NA>
#>  4 #ABC8E2 安徽省 340000 省     池州… 341700 地级市 <NA>  <NA>  <NA>       <NA>
#>  5 #A7C3E1 安徽省 340000 省     滁州… 341100 地级市 <NA>  <NA>  <NA>       <NA>
#>  6 #ABC7E2 安徽省 340000 省     阜阳… 341200 地级市 <NA>  <NA>  <NA>       <NA>
#>  7 #A0BEE8 安徽省 340000 省     合肥… 340100 地级市 <NA>  <NA>  <NA>       <NA>
#>  8 #B0CDEA 安徽省 340000 省     淮北… 340600 地级市 <NA>  <NA>  <NA>       <NA>
#>  9 #ACCAEB 安徽省 340000 省     淮南… 340400 地级市 <NA>  <NA>  <NA>       <NA>
#> 10 #A2BDDA 安徽省 340000 省     黄山… 341000 地级市 <NA>  <NA>  <NA>       <NA>
#> # … with 396 more rows, and 5 more variables: lng <dbl>, lat <dbl>, FID <dbl>,
#> #   id <chr>, geometry <GEOMETRY [m]>

# 小地图框格
small_bbox %>%
  as_tibble() %>%
  st_sf() %>%
  mutate(geometry = st_cast(geometry, "MULTILINESTRING")) %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  st_set_crs(mycrs) -> small_bboxsf

# 绘图
library(ggspatial)
mainchina %>%
  dplyr::filter(!is.na(hex)) %>%
  ggplot() +
  geom_sf(aes(fill = I(hex)), color = "gray80", size = 0.01) +
  geom_sf(data = subset(mainchina, name == "省级边界"),
          color = "black", size = 0.1) +
  geom_sf(data = mainchina %>%
            dplyr::filter(name == "九段线"),
          color = "black", size = 0.2) +
  geom_sf(data = mainchina %>%
            dplyr::filter(name == "海岸线"),
          color = "#0055AA", size = 0.2) +
  geom_sf(data = small_bboxsf, fill = NA,
          color = "black", size = 0.2) +
  annotation_scale(
    width_hint = 0.2,
    text_family = cnfont
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "false",
    width = unit(1.6, "cm"),
    height = unit(2, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) +
  theme(panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) -> map

ggdraw() +
  draw_plot(map) +
  draw_plot(legend,
            0.12, 0.2, 0.15, 0.15) -> res

ggsave("dt.png", width = 9, height = 9, device = png)

这样我们就解决了这个问题。

点击这里跳转到 RStata 短书平台获取附件：R 语言：如何从地图图片上提取数据并重新绘图？