QGIS & R 语言：如何从地图图片上提取数据并重新绘图？（使用 QGIS 进行地图配准）

之前给大家讲解过如何使用 R 语言从地图图片中提取数据的方法，不过那套方法的使用情形过于狭窄（必须要知道地图的坐标参考系才行）。最近发现还是使用 QGIS 进行地图的配准更为方便，今天我们就结合 QGIS 和 R 语言的优势来学习下如何使用 QGIS 进行地图的配准以及如何使用 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 多个城市的颜色。这里我提供一种思路供大家参考。

我的思路是先把这个图片作为栅格数据读取（需要预先使用 QGIS 对地图图片进行配准），然后再使用城市地理矢量数据提取每个区域的颜色众数，最后再把这些颜色绘制到地图上。

首先我们需要在电脑上安装 QGIS，这是个免费的开源软件，可以从这里下载安装：https://qgis.org/zh-Hans/site/

第一次打开 QGIS，我们可以先把界面语言设置成中文的：

首先是找到偏好设置：

然后在 常规 里面设置语言：

设置完之后需要重启才能生效。

然后再次打开 QGIS，找到配准工具：

再回到前面的论文图片，我们截取第一幅保存为：pic_origin.png

在 QGIS 的配准工具界面加载这幅图片：

然后就可以添加地面控制点了（具体再在视频讲解里面演示吧），为了方便查看地图上不同点的坐标，我们可以使用下面这个工具：

# 交互式查看中国地图
library(sf)
library(tidyverse)
read_sf("2021行政区划/省.shp") -> prov

library(leaflet)
# 智图地图底图
url <- "http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}"
leaflet() %>%
  addTiles(url) -> map
mapview::mapview(prov, map = map)

地面控制点至少需要选择六个，这里我选择了 12 个：

然后在设置里面我们选择变换方法为 多项式3：

这里的变换方法可以多多尝试，知道得到最合适的。最后点击开始配准就可以得到配准后的 tif 文件了，我们可以在 R 语言中读取比较下：

library(terra)
read_sf("2021行政区划/市.shp") %>%
  st_simplify(dTolerance = 2000) %>%
  vect() -> citysim

rast("pic_origin_已更改.tif") -> rst
plot(rst$pic_origin_已更改_1)
plot(citysim, add = T)

可以看到配准的效果还是不错的。下面我们就可以使用 R 语言提取里面每个城市中颜色的众数了。exactextractr 包的 exact_extract() 函数可以提取每个城市区域颜色值的众数（使用 fun = “mode”）：

# 分城市提取每个城市的值
library(exactextractr)
read_sf("2021行政区划/市.shp") -> city
bind_cols(
  city %>%
    st_drop_geometry() %>%
    select(-contains("类型")),
  exact_extract(rst$pic_origin_已更改_1, city,
                fun = "mode", progress = F) %>%
    as_tibble() %>%
    rename(R = value),
  exact_extract(rst$pic_origin_已更改_2, city,
                fun = "mode", progress = F) %>%
    as_tibble() %>%
    rename(G = value),
  exact_extract(rst$pic_origin_已更改_3, city,
                fun = "mode", progress = F) %>%
    as_tibble() %>%
    rename(B = value)
) -> citydf
citydf

#> # A tibble: 371 × 7
#>    省     省代码 市     市代码     R     G     B
#>    <chr>   <dbl> <chr>   <dbl> <dbl> <dbl> <dbl>
#>  1 安徽省 340000 安庆市 340800   177   208   235
#>  2 安徽省 340000 蚌埠市 340300   176   206   235
#>  3 安徽省 340000 亳州市 341600   176   206   235
#>  4 安徽省 340000 池州市 341700   176   206   235
#>  5 安徽省 340000 滁州市 341100   176   206   235
#>  6 安徽省 340000 阜阳市 341200   176   206   235
#>  7 安徽省 340000 合肥市 340100   176   174   235
#>  8 安徽省 340000 淮北市 340600   177   206   235
#>  9 安徽省 340000 淮南市 340400   176   206   235
#> 10 安徽省 340000 黄山市 341000   176   206   235
#> # … with 361 more rows

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

library(farver)
library(tidyverse)
citydf %>%
  mutate_at(c("R", "G", "B"), as.integer) -> citydf

citydf$hex <- encode_colour(citydf[, 5:7])

如果颜色种类较少的话，可以手动根据颜色对数据进行分组，但是这里颜色种类较多，不好一一分辨颜色了，所以我们就直接用这个数据绘图吧。

city %>%
  st_simplify(dTolerance = 2000) -> citysim

citydf %>%
  left_join(citysim) %>%
  st_sf() -> citydf2

citydf2 %>%
  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)

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

# 再次绘制主图
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"
citydf2 %>%
  st_transform(mycrs) -> citydf

# 省份边界
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_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){
  mainchina[mainchina$id == i,]$geometry <- st_combine(
    c(mainchina[mainchina$id == i,]$geometry,
      nanhai[nanhai$id == i,]$geometry)
  )
}

mainchina

# 小地图框格
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 短书平台获取附件：QGIS & R 语言：如何从地图图片上提取数据并重新绘图？（使用 QGIS 进行地图配准）