使用 R 语言绘制中国乡镇地图及根据经纬度判断所处的乡镇

在之前的课程中，我给大家分享过使用 R 语言绘制省市区县地图的方法：

• 使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）：https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c

另外也给大家分享过使用 R 语言进行地理编码以及根据经纬度判断所处省市区县的方法：

• 使用 R 语言进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县：https://rstata.duanshu.com/#/brief/course/6f8633b486ed49e398c486d0d9f0f59a

乡镇地图的绘制

今天我们再来补充讲解下如何绘制乡镇地图以及如何根据经纬度判断所处的乡镇。

附件中的 town4326 文件夹是 WGS84 坐标系下的乡镇地图数据；town_mini 文件夹是经过编辑后的小地图版本的乡镇地图数据；town_long 则是长版本的。

town_mini 和 town_long 都是“+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”坐标系的，也就是和之前的地图数据一样。

首先加载相关 R 包并读取数据：

library(tidyverse)
library(sf)

# 中国地图通常使用这样的坐标系
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"

# 小地图版本
# 读取编辑后的乡镇地图数据
read_sf("town_mini/town_mini.shp") -> townmini

# 线条数据可以使用之前课程中提供的省市区县的：
read_sf("chinaprov2021mini/chinaprov2021mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> provlinemap

# 读取绘图数据
haven::read_dta("2019年各乡镇平均人口密度.dta") %>%
  rename(townname = 乡镇, towncode = 乡镇代码) -> df

df

#> # A tibble: 43,366 × 3
#>    townname     towncode       mean
#>    <chr>        <chr>         <dbl>
#>  1 一元街办事处 420102003000 9076.
#>  2 一六镇       431022108000  295.
#>  3 一六镇       440232103000   89.2
#>  4 一农场       130209401000  346.
#>  5 一品街道     500113008000  557.
#>  6 一市镇       330226104000  426.
#>  7 一平垣乡     141002203000  677.
#>  8 一平浪镇     532331105000   72.9
#>  9 一座营镇     220581109000  252.
#> 10 一心乡       230624200000   45.7
#> # ℹ 43,356 more rows

然后就可以绘图了，这部分细节较多，建议结合视频讲解学习：

library(ggspatial)
library(ggnewscale)
townmini %>%
  left_join(df) %>%
  mutate(mean = if_else(is.na(mean), -1, mean)) %>%
  mutate(group = cut(mean, breaks = c(-1, -0.1, 68, 125, 197, 312,
                                      458, 663, 1360, 50000),
                     include.lowest = T,
                     labels = c("无数据", "<=68",
                                "68~125", "125~197",
                                "197~312", "312~458",
                                "458~663", "663~1360",
                                ">1360"))) -> mapdata

mapdata %>%
  ggplot() +
  geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black",
               "省份" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3,
               "省份" = 0.3)
  ) +
  scico::scale_fill_scico_d(palette = "bamako",
                          name = "人口密度(人/m2)",
                          direction = -1) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "2019 年各乡镇平均人口密度",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：中国科学院资源环境科学与数据中心") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p1

ggsave("pic1.png", width = 10, height = 8.5, device = png)

长版地图的绘制代码类似，数据换成长版的即可：

# 长版
read_sf("town_long/town_long.shp") -> townlong
townlong %>%
  left_join(df) %>%
  mutate(mean = if_else(is.na(mean), -1, mean)) %>%
  mutate(group = cut(mean, breaks = c(-1, -0.1, 68, 125, 197, 312,
                                      458, 663, 1360, 50000),
                     include.lowest = T,
                     labels = c("无数据", "<=68",
                                "68~125", "125~197",
                                "197~312", "312~458",
                                "458~663", "663~1360",
                                ">1360"))) -> mapdata

# 线条
read_sf("chinaprov2021long/chinaprov2021long_line.shp") %>%
  filter(class %in% c("九段线", "海岸线")) %>%
  select(class) -> provlinemap

ggplot(mapdata) +
  geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3)
  ) +
  scico::scale_fill_scico_d(palette = "bamako",
                            name = "人口密度(人/m2)",
                            direction = -1) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  coord_sf(xlim = c(-2625585.8, 2206964.7) * 1.2) +
  labs(title = "2019 年各乡镇平均人口密度",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：中国科学院资源环境科学与数据中心") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p2

ggsave("pic2.png", width = 8, height = 7, device = png)

根据经纬度判断所处的乡镇

首先读取乡镇矢量数据：

read_sf("town4326/town4326.shp") -> town

工企样本.dta 是大概 1 万条工企数据的样本，含有经纬度变量，使用下面的代码就可以判断每个坐标点所处的乡镇了：

haven::read_dta("工企样本.dta") %>%
  filter(!is.na(经度)) %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326, remove = F) -> dfres

如果数据量不大，可以直接用下面的代码判断：

# dfres %>%
#   st_intersection(town) -> dfres

数据量大的时候会非常慢，建议转换成 mycrs 坐标系再运算：

# 建议转换成 mycrs 坐标系再判断：
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"
dfres %>%
  st_transform(mycrs) %>%
  st_intersection(st_transform(town, mycrs)) -> dfres2

dfres2 %>%
  st_drop_geometry()
#> # A tibble: 10,218 × 7
#>     group       gqid 企业名称                  经度  纬度 Name       code
#>  *  <dbl>      <dbl> <chr>                    <dbl> <dbl> <chr>      <chr>
#>  1 719089 1998169629 天津酿酒厂                117.  39.2 丁字沽街道 1201060040…
#>  2 721426 2008228990 大连水泥集团有限公司      122.  39.3 七顶山街道 2102130210…
#>  3 721426 2009042065 大连水泥集团有限公司      122.  39.3 七顶山街道 2102130210…
#>  4 721426 2014087208 大连水泥集团有限公司      122.  39.3 七顶山街道 2102130210…
#>  5 719068 1998169308 天津市神光新技术开发公司  117.  39.1 万兴街道   1201040050…
#>  6 718486 1999000502 北京金盾印刷厂            116.  39.9 万寿路街道 1101080010…
#>  7 718486 2004003087 北京金盾印刷厂            116.  39.9 万寿路街道 1101080010…
#>  8 718486 2007333827 北京金盾印刷厂            116.  39.9 万寿路街道 1101080010…
#>  9 718486 2003004268 北京金盾印刷厂            116.  39.9 万寿路街道 1101080010…
#> 10 718486 2013140022 北京金盾印刷厂            116.  39.9 万寿路街道 1101080010…
#> # ℹ 10,208 more rows

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制中国乡镇地图及根据经纬度判断所处的乡镇