使用 R 语言爬取宗教活动场所基本信息

最近有个小伙伴想爬取国家宗教事务局官网上的宗教场所基本信息数据，今天就让我们一起来学习下这个数据的爬取，然后我们再根据地址变量使用高德地图和百度地图地理编码接口解析经纬度并绘制地图展示。

国家宗教事务局: https://www.sara.gov.cn/gjzjswj/zjjcxxcxxt/zjhdcsjbxx/index.shtml

这个网站上的数据是这样的：

是一个表格，进行网页分析可以发现这个表格数据可以使用如下 curl 语句请求到：

curl $'https://www.sara.gov.cn/eportal/ui?portal.url=/portlet/place-info-front\u0021queryList.portlet&moduleId=b68f71decbfb4ff4b274cac2bf67204b&pageId=614c812b4af64b4e91e31a7a7bc9f12a' \
  -H 'Accept: application/json, text/javascript, */*; q=0.01' \
  -H 'Accept-Language: zh-CN,zh;q=0.9,en;q=0.8' \
  -H 'Connection: keep-alive' \
  -H 'Content-Type: application/x-www-form-urlencoded; charset=UTF-8' \
  -H 'Cookie: JSESSIONID=11BB3892315623DF0B6C70A0E234A9A1; __jsluid_s=6e96424b5d1ddbcd09accce1cf4632ef; https_waf_cookie=0d1569c8-c3c1-4a9ebbbb6a7104f3119d0974b4ebb790fce1; HWWAFSESID=9f1a2739aebe5a90c6; HWWAFSESTIME=1675322361777' \
  -H 'Origin: https://www.sara.gov.cn' \
  -H 'Referer: https://www.sara.gov.cn/gjzjswj/zjjcxxcxxt/zjhdcsjbxx/index.shtml' \
  -H 'Sec-Fetch-Dest: empty' \
  -H 'Sec-Fetch-Mode: cors' \
  -H 'Sec-Fetch-Site: same-origin' \
  -H 'User-Agent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/109.0.0.0 Safari/537.36' \
  -H 'X-Requested-With: XMLHttpRequest' \
  -H 'request-by: ajax-request-tag' \
  -H 'sec-ch-ua: "Not_A Brand";v="99", "Google Chrome";v="109", "Chromium";v="109"' \
  -H 'sec-ch-ua-mobile: ?0' \
  -H 'sec-ch-ua-platform: "macOS"' \
  --data-raw 'province=&city=&town=&faction=&religionName=&keyWord=&currentPage=1&pageSize=15' \
  --compressed

我们可以把这个 curl 语句改写成 R 语言的，不过直接改写比较麻烦，可以使用这个网页工具：https://curlconverter.com/r/

# 这些代码大家需要根据从自己电脑上复制 curl 代码并转换
library(tidyverse)
library(jsonlite)

# curl 语句转 R 语言：https://curlconverter.com/r/
library(httr)

cookies = c(
  `JSESSIONID` = "11BB3892315623DF0B6C70A0E234A9A1",
  `__jsluid_s` = "6e96424b5d1ddbcd09accce1cf4632ef",
  `https_waf_cookie` = "0d1569c8-c3c1-4a9ebbbb6a7104f3119d0974b4ebb790fce1",
  `HWWAFSESID` = "9f1a2739aebe5a90c6",
  `HWWAFSESTIME` = "1675322361777"
)

headers = c(
  `Accept` = "application/json, text/javascript, */*; q=0.01",
  `Accept-Language` = "zh-CN,zh;q=0.9,en;q=0.8",
  `Connection` = "keep-alive",
  `Content-Type` = "application/x-www-form-urlencoded; charset=UTF-8",
  `Origin` = "https://www.sara.gov.cn",
  `Referer` = "https://www.sara.gov.cn/gjzjswj/zjjcxxcxxt/zjhdcsjbxx/index.shtml",
  `Sec-Fetch-Dest` = "empty",
  `Sec-Fetch-Mode` = "cors",
  `Sec-Fetch-Site` = "same-origin",
  `User-Agent` = "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/109.0.0.0 Safari/537.36",
  `X-Requested-With` = "XMLHttpRequest",
  `request-by` = "ajax-request-tag",
  `sec-ch-ua` = '"Not_A Brand";v="99", "Google Chrome";v="109", "Chromium";v="109"',
  `sec-ch-ua-mobile` = "?0",
  `sec-ch-ua-platform` = '"macOS"'
)

data = list(
  `province` = "",
  `city` = "",
  `town` = "",
  `faction` = "",
  `religionName` = "",
  `keyWord` = "",
  `currentPage` = "1",
  `pageSize` = 5000
)

res <- httr::POST(url = "https://www.sara.gov.cn/eportal/ui?portal.url=/portlet/place-info-front!queryList.portlet&moduleId=b68f71decbfb4ff4b274cac2bf67204b&pageId=614c812b4af64b4e91e31a7a7bc9f12a", httr::add_headers(.headers=headers), httr::set_cookies(.cookies = cookies), body = data, encode = "form")

通过尝试我们可以发现 pageSize 参数直接改成 50000 就可以一次性请求到所有页面的数据了，这样就不用循环了。但是需要注意并非任何网站爬取都这样这样操作，有些网站会限制。

然后我们就可以使用 content() 函数提取 res 里面的数据了：

content(res, type = "text") %>%
  fromJSON() -> df
df$result %>%
  as_tibble() -> df

df
#> # A tibble: 5,000 × 10
#>    placeIn…¹ relig…² faction place…³ address perso…⁴ provi…⁵ city  town  chann…⁶
#>    <chr>     <chr>   <chr>   <chr>   <chr>   <chr>   <chr>   <chr> <chr> <chr>
#>  1 1         佛教    汉语系  广济寺  北京西… 演觉    北京市  市辖… 西城… 1935
#>  2 10        佛教    汉语系  双泉寺  北京市… 常藏    北京市  市辖… 石景… 1935
#>  3 100       佛教    汉语系  天光寺  上海市… 能照    上海市  市辖… 青浦… 1935
#>  4 1000      佛教    汉语系  宁海福… 浙江省… 释克凤  浙江省  宁波… 宁海… 1935
#>  5 10000     佛教    汉语系  会昌县… 江西省… 钟贵明  江西省  赣州… 会昌… 1935
#>  6 10001     佛教    汉语系  会昌县… 江西省… 释祖正  江西省  赣州… 会昌… 1935
#>  7 10002     佛教    汉语系  会昌县… 江西省… 蔡永增  江西省  赣州… 会昌… 1935
#>  8 10003     佛教    汉语系  会昌县… 江西省… 郭宜发  江西省  赣州… 会昌… 1935
#>  9 10004     佛教    汉语系  会昌县… 江西省… 黄绪桂  江西省  赣州… 会昌… 1935
#> 10 10005     佛教    汉语系  会昌县… 江西省… 吴才成  江西省  赣州… 会昌… 1935
#> # … with 4,990 more rows, and abbreviated variable names ¹placeInfoId,
#> #   ²religionName, ³placeName, ⁴personCharge, ⁵province, ⁶channelId

再整理下：

df %>%
  select(1:9) %>%
  set_names(c("id", "宗教", "派别", "场所名称", "地址", "负责人姓名",
              "省", "市", "县")) %>%
  mutate(id = as.numeric(id)) -> df

# 临时保存下
df %>%
  write_rds("df.rds")

这样实际上我们就已经完成了这个数据爬取，不过下面我们还可以根据地址变量解析经纬度。

关于地址解析经纬度的方法建议大家学习下平台上的地理编码教程（平台上搜索“地理编码”即可找到）。

首先第一轮我们多线程使用高德地图进行地理编码：

# 解析经纬度
library(tidyverse)
library(sf)
library(parallel)
makeCluster(36) -> cl
clusterExport(cl, "%>%")
dir.create("rds")
read_rds("df.rds") -> df
# 第一轮(多运行几遍)
df %>%
  mutate(市 = if_else(市 == "市辖区", 省, 市)) %>%
  mutate(市 = if_else(市 == "直管县", 县, 市)) %>%
  mutate(地址 = if_else(!str_detect(地址, 市), paste0(市, 地址), 地址)) %>%
  select(id, add = 地址) %>%
  mutate(add = str_remove_all(add, "[\\r\\n\\t\\s]"),
         add = str_replace_all(add, "#", "号")) -> df2

df2
clusterExport(cl, "df2")
parLapply(cl, seq(1, nrow(df2), by = 10), function(x){
  if(!file.exists(paste0("rds/", x, ".rds"))) {
    try({
      df2 %>%
        dplyr::slice(x:(x + 9)) -> temp

      temp %>%
        dplyr::pull(add) %>%
        paste0(collapse = "|") %>%
        paste0("https://restapi.amap.com/v3/geocode/geo?address=", .,
               "&key=ac724d548ac22be9abd641e86eff95f8&batch=true") -> url

      # 读取解析 json 文件
      jsonlite::fromJSON(url) -> list
      list$geocodes %>%
        tibble::as_tibble() %>%
        dplyr::select(province, citycode, city, district, location, level) -> geodf

      # 与之前的 temp 合并
      geodf %>%
        dplyr::bind_cols(temp, .) %>%
        dplyr::mutate_all(as.character) %>%
        readr::write_rds(paste0("rds/", x, ".rds"))
    })
  }
}) -> res

# 合并成功的
fs::dir_ls("rds") %>%
  parLapply(cl, ., read_rds) %>%
  bind_rows() -> dfa1

dfa1

这里的 makeCluster(36) -> cl 是创建了 36 进程，这个进程数需要根据大家自己的电脑进行调整，可以使用 parallel::detectCores() 函数查看自己电脑的进程数。

然后我们筛选第一轮解析没有成功的进行第二轮解析，由于没成功的数量不多了，所以我们就直接用百度地图地理编码接口进行解析：

# 第二轮
df2 %>%
  anti_join(
    dfa1 %>%
      dplyr::filter(str_detect(location, ",")) %>%
      select(id) %>%
      mutate(id = as.numeric(id))
  ) -> df2

df2 %>%
  inner_join(df) %>%
  select(-add) %>%
  select(id, add = 地址) -> df2

# 剩下的用百度地图
df2 %>%
  mutate(location = map_chr(add, function(x) {
    try({
      paste0("https://api.map.baidu.com/geocoding/v3/?address=", x,
             "&output=json&ak=KNfhEZFNPTNL5Gyjm3O81nlE5mURbOhn&ret_coordtype=gcj02ll") %>%
        fromJSON() -> lst
      paste0(lst$result$location, collapse = ",")
    })
  })) -> dfa2

dfa2

最后还有一些没有成功的，我们单独处理下就行：

# 最后仍然失败的，单独处理下（百度下其他地址）
dfa2 %>%
  dplyr::filter(!str_detect(location, ",")) %>%
  mutate(add = "仙桃市邓李湾村") %>%
  mutate(location = map_chr(add, function(x) {
    try({
      paste0("https://api.map.baidu.com/geocoding/v3/?address=", x,
             "&output=json&ak=KNfhEZFNPTNL5Gyjm3O81nlE5mURbOhn&ret_coordtype=gcj02ll") %>%
        fromJSON() -> lst
      paste0(lst$result$location, collapse = ",")
    })
  })) -> dfa3

合并所有成功的结果：

# 合并所有的结果
bind_rows(
  dfa1 %>%
    dplyr::filter(str_detect(location, ",")) %>%
    mutate(id = as.numeric(id)),
  dfa2 %>%
    dplyr::filter(str_detect(location, ",")) %>%
    mutate(id = as.numeric(id)),
  dfa3
) -> dfall

不过这样得到的数据是 GCJ02 坐标的，我们还需要转换成 WGS84 坐标才能使用：

source("GCJ02坐标转WGS84.R")

dfall %>%
  select(id, location) %>%
  separate(location, into = c("经度", "纬度"), sep = ",") %>%
  mutate(经度 = as.numeric(经度),
           纬度 = as.numeric(纬度),
           value2 = map2_chr(经度, 纬度, GCJ02_WGS84)) %>%
  select(-contains("度")) %>%
  separate(value2, into = c("经度", "纬度"), sep = ",") %>%
  type_convert() %>%
  inner_join(df) -> dfall2

dfall2

有了经纬度就可以根据经纬度判断每个宗教场所所处的省市区县了：

# 根据经纬度判断所处的省市区县
read_sf("2021行政区划/县.shp") %>%
  st_transform(3055) -> county

dfall2 %>%
  rename(省_orign = 省, 市_origin = 市, 县_origin = 县) -> dfall2

# 经纬度转 sf 对象
dfall2 %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326, remove = F) %>%
  st_transform(3055) -> dfsf

dfsf %>%
  st_intersection(county) -> dfsf2

dfsf2 %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  rename(ID = id) -> df

df

df %>%
  writexl::write_xlsx("宗教活动场所基本信息（含经纬度及其所处的省市区县）.xlsx")

这个时候实际上还可以继续检查下解析地址的准确性问题，不过这里因为时间的关系我们就不再继续检查了。

还可以保存成 dta 格式的：

df %>%
  haven::write_dta("宗教活动场所基本信息（含经纬度及其所处的省市区县）.dta",
  label = "爬取&整理：微信公众号 RStata")

最后我们再把这些经纬度坐标绘制到地图上：

# 绘图展示
library(tidyverse)
library(sf)
readxl::read_xlsx("宗教活动场所基本信息（含经纬度及其所处的省市区县）.xlsx") %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326, remove = F) %>%
  mutate(type = paste0(宗教, " - ", 派别)) -> df
read_sf("2021行政区划/省.shp") -> prov
read_sf("九段线/九段线.shp") -> jdx
read_sf("海岸线/海岸线.shp") -> hax
source("theme_map.R")
library(ggspatial)

ggplot(prov) +
  geom_sf(fill = NA, linewidth = 0.2,
          color = "gray") +
  geom_sf(data = jdx, linewidth = 0.5,
          color = "black", fill = NA) +
  geom_sf(data = hax, linewidth = 0.3,
          color = "#0055AA", fill = NA) +
  geom_sf(data = df, aes(color = type), alpha = 0.5,
          shape = 15, size = 1) +
  scale_color_manual(values = c("#D7BA9F", "#800080", "#1C3181",  "#1BB6AF",  "#FFAD0A")) +
  coord_sf(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",
           xlim = c(-3500000, 3090000)) +
  scale_x_continuous(expand = c(0.001, 0.001)) +
  scale_y_continuous(expand = c(0.001, 0.001)) +
  guides(fill = guide_legend(nrow = 1,
                             label.position = "top")) +
  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(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        plot.background = element_rect(color = "gray")) +
  labs(title = "中国宗教活动场所地理分布",
       subtitle = "绘制：微信公众号 RStata",
       caption = "数据来源：国家宗教事务局\n<http://www.sara.gov.cn/zjhdcsjbxx/index.jhtml>",
       color = "宗教") +
  guides(color = guide_legend(nrow = 2)) -> p

ggsave("中国宗教活动场所地理分布.png",
       width = 9, height = 9, device = png)

knitr::plot_crop("中国宗教活动场所地理分布.png")

这样我们就完成了这个任务。

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言爬取宗教活动场所基本信息