由于高德地图接口的更新，之前的代码不再适用了，主要是两个更新：
- 高德地图的接口会根据一个地址返回多个结果；
- 高德地图不再建议一次解析多个地址了，会返回混乱的结果。
- 如果需要进行大量的解析，可以联系李老师帮忙，费用参考这个：https://rstata.duanshu.com/#/course/6b44ce2712af47b6ba7b2a2fd74a4e68
之前给大家讲解过根据地址解析经纬度并根据经纬度判断所处的省市区县的方法，不过那个方法效率较低，面对数百万地址的工作就显得力不从心了。今天再给大家介绍一套无比高效的解析流程。
这套工作流程包含下面几个过程：
- 使用高德地图地理编码接口解析（多线程）；
- 反复重复上述过程，直到筛选出那些无法使用高德地图地理编码接口解析成功的地址；
- 然后使用百度地图地理编码接口解析上面的操作没有成功的部分；
- 合并所有的解析结果；
- GCJ02 坐标转换成 WGS84 坐标；
- 使用地理计算根据经纬度判断所处的省市区县。
下面就让我们以金融机构网点数据的地理编码为例进行讲解。
加载所需的 R 包
library(tidyverse)
library(jsonlite)
library(parallel)
|
读取金融机构网点数据
readxl::read_xlsx("金融机构网点数据.xlsx") -> df
df %>%
slice(1:1000) -> df
df %>%
mutate(id = row.names(.)) %>%
select(id, everything()) -> df
|
我们下面就以这一千个地址为例进行讲解。
选择感兴趣的变量：
df %>%
select(id, 机构名称, 机构所在地, 机构地址, 地址代码) -> dfsim
|
申请高德地图和百度地图的密钥
- 高德地图文档：https://lbs.amap.com/api/webservice/guide/api/georegeo
- 百度地图文档：https://lbsyun.baidu.com/index.php?title=webapi/guide/webservice-geocoding-base
高德地图的 city 参数既可以是城市名称也可以是城市区划代码，百度地图的 city 参数只能是城市名称，所以我们需要准备两个 city 变量。
准备数据
dfsim %>%
mutate(address = 机构地址,
citygd = 地址代码) %>%
mutate(citybd = str_match(机构所在地, "-+(.*)")[,2],
citybd = if_else(is.na(citybd), 机构所在地, citybd)) %>%
select(id, 机构名称, address, citybd, citygd) -> dfsim
dfsim
|
使用高德地图解析
在这套工作流程中我们会优先使用高德地图地理编码接口进行解析，因为高德地图支持更高的 QPS（每秒请求的次数）。不过不幸运的是现在的高德和百度都只提供每日 5000 次的免费额度，难以进行大量的解析。
这里插播一条广告。如果大家需要大量的解析，可以考虑联系李老师帮忙解析。费用可以参考这个：https://rstata.duanshu.com/#/course/6b44ce2712af47b6ba7b2a2fd74a4e68
例如解析 1 条地址：
dfsim %>%
mutate(address = paste0(address, 机构名称),
address = str_replace_all(address, "#", "号")) -> dfsim2
paste0("https://restapi.amap.com/v3/geocode/geo?address=", dfsim2$address[1],
"&key=50828d53749c02431a65192f7c092a77&city=", dfsim2$citygd[1]) -> url
fromJSON(url) -> list
list$geocodes %>%
as_tibble() %>%
select(formatted_address, province, city,
district, location, level) -> geodf
geodf$id <- dfsim2$id[1]
geodf
|
这样就完成了 1 个地址的解析，得到了三个结果。我们先把所有的结果都保留，最后再筛选。
由此我们就可以设计一个循环处理所有的了：
dir.create("rds1")
lapply(1:nrow(dfsim2), function(x){
if(!file.exists(paste0("rds1/", x, ".rds"))) {
suppressMessages({
suppressWarnings({
try({
paste0("https://restapi.amap.com/v3/geocode/geo?address=", dfsim2$address[x],
"&key=50828d53749c02431a65192f7c092a77&city=", dfsim2$citygd[x]) -> url
fromJSON(url) -> list
list$geocodes %>%
as_tibble() %>%
select(formatted_address, province, city,
district, location, level) -> geodf
geodf$id <- dfsim2$id[x]
geodf %>%
mutate_all(as.character) %>%
write_rds(paste0("rds1/", x, ".rds"))
})
})
})
}
}) -> res
|
多线程解析
搭配上多线程，这个解析效率就会被极大的提高，当然这个线程数量得根据你自己电脑的实际情况来设定。
# 查看自己电脑的最大线程数
detectCores()
# 创建线程，也别太高：
makeCluster(36) -> cl
# 分发对象
clusterExport(cl, "dfsim2")
|
clusterEvalQ(cl, expr = ({
library(tidyverse)
library(jsonlite)
}))
parLapply(cl, 1:nrow(dfsim2), function(x){
if(!file.exists(paste0("rds1/", x, ".rds"))) {
suppressMessages({
suppressWarnings({
try({
paste0("https://restapi.amap.com/v3/geocode/geo?address=", dfsim2$address[x],
"&key=50828d53749c02431a65192f7c092a77&city=", dfsim2$citygd[x]) -> url
fromJSON(url) -> list
list$geocodes %>%
as_tibble() %>%
select(formatted_address, province, city,
district, location, level) -> geodf
geodf$id <- dfsim2$id[x]
geodf %>%
mutate_all(as.character) %>%
write_rds(paste0("rds1/", x, ".rds"))
})
})
})
}
}) -> res
|
合并结果：
parLapply(cl, fs::dir_ls("rds1"), readr::read_rds) %>%
bind_rows() -> df1
df1
|
这一轮中解析成功的：
df1 %>%
inner_join(dfsim2) %>%
dplyr::filter(str_detect(location, ",")) %>%
dplyr::filter(str_sub(city, 1, 2) == str_sub(citybd, 1, 2)) %>%
dplyr::filter(!level %in% c("省", "市", "县")) -> df1success
|
筛选出没成功的：
dfsim2 %>%
anti_join(select(df1success, id)) -> dfsim2
dfsim2
|
也可以对解析没成功的地址进行调整，反复使用上面的代码尝试，最后实在无法成功的可以用百度地图再试试。
使用百度地图解析
以“新疆石河子市东环路17号”为例：
paste0("https://api.map.baidu.com/geocoding/v3/?address=", "新疆石河子市东环路17号",
"&output=json&ak=h4jN43d3ajnnXGzBHmIowERWwM2o6PIe&ret_coordtype=gcj02ll&city=", "石河子市") -> url
url %>%
fromJSON() -> lst
tibble(
location = paste0(lst$result$location, collapse = ","),
level = lst$result$level,
precise = lst$result$precise,
confidence = lst$result$confidence,
comprehension = lst$result$comprehension
)
|
于是就可以编写程序使用百度地图地理编码接口解析最后的了，百度地图的 QPS 较低，可以设置 5 进程：
dfsim2 %>%
mutate(address = str_remove_all(address, " ")) -> dfsim2
makeCluster(5) -> cl
clusterExport(cl, "dfsim2")
dir.create("rds2")
parLapply(cl, 1:nrow(dfsim2), function(x){
if(!file.exists(paste0("rds2/", dfsim2$id[x], ".json"))) {
try({
download.file(paste0("https://api.map.baidu.com/geocoding/v3/?address=", dfsim2$address[x],
"&output=json&ak=h4jN43d3ajnnXGzBHmIowERWwM2o6PIe&ret_coordtype=gcj02ll"),
paste0("rds2/", dfsim2$id[x], ".json"))
})
}
}) -> res
|
合并结果：
fs::dir_ls("rds2") %>%
parLapply(cl, ., function(x){
jsonlite::fromJSON(x) -> lst
tibble::tibble(
file = x,
location = paste0(lst$result$location, collapse = ","),
level = lst$result$level,
precise = lst$result$precise,
confidence = lst$result$confidence,
comprehension = lst$result$comprehension
)
}) %>%
bind_rows() -> dfres2
dfres2 %>%
mutate(file = str_match(file, "/(\\d+)")[,2]) %>%
rename(id = file) -> dfres2
dfsim2 %>%
anti_join(dfres2 %>% select(id))
dfres2 %>%
select(id, location) -> df2success
|
合并两轮结果
然后就可以一股脑的合并了：
bind_rows(df1success, df2success) %>%
select(id, location) %>%
group_by(id) %>%
slice(1) %>%
ungroup() -> dfall
dfall
|
GCJ02 转 WGS84
解析的结果是 GCJ02 坐标系的，不能直接使用，需要转换成 WGS84。
坐标转换.R 文件中有转换函数，具体函数的编写方法不重要，会使用即可：
source("坐标转换.R")
dfall %>%
separate(location, into = c("lng", "lat"), sep = ",", remove = F) %>%
mutate(经度 = as.numeric(lng),
纬度 = as.numeric(lat),
value2 = map2_chr(经度, 纬度, GCJ02_WGS84)) %>%
select(-contains("度")) %>%
separate(value2, into = c("经度", "纬度"), sep = ",") %>%
type_convert() %>%
select(-lat, -lng) %>%
select(id, lat = 纬度, lng = 经度, everything()) -> dfall2
dfall2
|
这样我们就完成了这一千个地址解析经纬度的工作。
根据经纬度判断所处的省市区县
虽然上面的代码返回的有省市区县的结果，但是不够“整洁”，所以还不如用地理计算判断下。
这里以 2021 年的数据为例，这里分别把地图数据和坐标数据都转换成了 mycrs 坐标系，这样进行判断的会更快：
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("2021行政区划/县.shp") %>%
st_transform(mycrs) -> county
dfall2 %>%
select(-location) %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326, remove = F) %>%
st_transform(mycrs) %>%
st_intersection(county) %>%
st_drop_geometry() %>%
select(-contains("类型")) %>%
rename(经度 = lng, 纬度 = lat) -> dfall3
dfall3
|
最后再合并原数据即可：
df %>%
select(id, 机构名称) %>%
left_join(dfall3 %>% mutate(id = as.character(id))) -> dfres
dfres
dfres %>%
haven::write_dta("经纬度解析结果.dta")
|
这样我们就完成了经纬度解析。当然时间操作中经常会遇到更糟糕的情况，例如地址很乱、只有企业名称等。这个时候要非常谨慎，特别是对于只有企业名称的地址，建议先根据企业名称匹配工商注册信息，再根据工商注册信息中的地址解析经纬度。
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县
评论