获取经纬度
library(sf)
library(terra)
library(tidyverse)
library(readxl)
library(duckdb)
library(httr) # 必须加载，提供GET()函数
library(jsonlite) # 解析JSON返回结果
library(arrow)
# 读取接种单位数据
jzdw <- read_xlsx('./接种单位列表.xlsx') |>
select(接种单位名称, 接种单位所属区县名称) |>
distinct()
get_jzdw_location_amap <- function(
jzdw_df,
amap_key,
max_retry = 3,
base_timeout = 15
) {
# 步骤1：拼接地址
jzdw_df <- jzdw_df |>
mutate(
接种单位所属区县名称 = ifelse(
is.na(接种单位所属区县名称) | 接种单位所属区县名称 == "",
"",
接种单位所属区县名称
),
完整地址 = paste0("甘肃省", 接种单位所属区县名称, 接种单位名称) |>
str_squish(),
地址编码 = URLencode(完整地址, reserved = TRUE),
行号 = row_number() # 添加行号，方便跟踪进度
)
# 辅助函数：安全提取字段（确保返回标量字符，而非列表）
safe_extract <- function(x, field) {
# 如果x是列表/数据框，提取字段；否则返回NA
val <- tryCatch(
{
if (is.list(x) && field %in% names(x)) {
# 提取后强制转为字符，若为空则返回NA
res <- as.character(x[[field]])
if (length(res) == 0 || res == "" || is.null(res)) NA else res
} else {
NA_character_
}
},
error = function(e) {
NA_character_
}
)
# 确保最终返回单个值（非列表）
if (length(val) > 1) val[1] else val
}
# 步骤2：单个地址编码（核心修复：字段类型统一）
geocode_single <- function(address_encoded, key, row_num, full_address) {
# 1秒调用间隔（保留）
Sys.sleep(1)
# 兜底返回值（所有字段均为字符/数值标量）
default <- tibble(
经度 = as.numeric(NA),
纬度 = as.numeric(NA),
匹配地址 = NA_character_,
country = NA_character_,
province = NA_character_,
city = NA_character_,
district = NA_character_,
调用状态 = "失败-未尝试"
)
# 重试机制
for (retry in 1:max_retry) {
tryCatch(
{
# 构造URL
url <- sprintf(
"https://restapi.amap.com/v3/geocode/geo?address=%s&key=%s&output=JSON",
address_encoded,
key
)
# 发送请求（延长超时时间）
resp <- GET(
url,
timeout(base_timeout),
add_headers(
'User-Agent' = 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36'
)
)
# 检查响应状态
if (http_error(resp)) {
message(sprintf(
"[第%s行][重试%s] %s - HTTP错误：%s",
row_num,
retry,
full_address,
http_status(resp)$message
))
next
}
# 解析JSON（强制转为列表，避免特殊格式）
json_text <- content(resp, "text", encoding = "UTF-8")
json_data <- fromJSON(json_text, simplifyVector = TRUE)
if (is.null(json_data) || json_data$status != "1") {
message(sprintf(
"[第%s行][重试%s] %s - API返回失败：%s",
row_num,
retry,
full_address,
json_data$info
))
next
}
# 提取地理编码信息（确保是标量）
geocodes <- json_data$geocodes
if (length(geocodes) == 0 || !is.list(geocodes)) {
message(sprintf(
"[第%s行][重试%s] %s - 无有效地理编码",
row_num,
retry,
full_address
))
next
}
first_code <- geocodes[1, ] # 只取第一个结果，强制转为单行
# 解析经纬度（确保是数值标量）
lng_lat <- str_split(safe_extract(first_code, "location"), ",")[[1]]
lng <- if (length(lng_lat) >= 1) as.numeric(lng_lat[1]) else NA_real_
lat <- if (length(lng_lat) >= 2) as.numeric(lng_lat[2]) else NA_real_
# 成功返回结果（所有字段均为标量，类型统一）
result <- tibble(
经度 = lng,
纬度 = lat,
匹配地址 = safe_extract(first_code, "formatted_address"),
country = safe_extract(first_code, "country"),
province = safe_extract(first_code, "province"),
city = safe_extract(first_code, "city"),
district = safe_extract(first_code, "district"), # 核心修复：确保是字符标量
调用状态 = "成功"
)
message(sprintf("[第%s行] %s - 处理成功", row_num, full_address))
return(result)
},
error = function(e) {
# 捕获超时/连接错误
message(sprintf(
"[第%s行][重试%s] %s - 错误：%s",
row_num,
retry,
full_address,
e$message
))
if (retry == max_retry) {
default$调用状态 <- paste0("失败-最终超时/错误：", e$message)
}
}
)
}
# 所有重试都失败，返回兜底值
return(default)
}
# 步骤3：批量处理（确保数据类型一致）
message(sprintf(
"开始处理%s条地址，每次调用间隔1秒，最多重试%s次",
nrow(jzdw_df),
max_retry
))
jzdw_with_loc <- jzdw_df |>
rowwise() |>
do({
res <- geocode_single(.$地址编码, amap_key, .$行号, .$完整地址)
bind_cols(., res)
}) |>
ungroup() |>
select(
行号,
接种单位名称,
接种单位所属区县名称,
完整地址,
经度,
纬度,
匹配地址,
country,
province,
city,
district,
调用状态
)
# 输出失败统计
fail_count <- sum(jzdw_with_loc$调用状态 != "成功")
message(sprintf(
"处理完成！成功：%s条，失败：%s条",
nrow(jzdw_df) - fail_count,
fail_count
))
return(jzdw_with_loc)
}
# 替换为你的高德API Key
api_key <- "XXXX"
# 调用函数（调整参数：最多重试3次，超时时间15秒）
jzdw_with_location <- get_jzdw_location_amap(
jzdw,
api_key,
max_retry = 3, # 每个地址最多重试3次
base_timeout = 15 # 单次请求超时时间延长到15秒
)
jzdw_with_location |>
write_parquet('./接种单位位置.parquet')
jzdwcs <- read_parquet('./接种单位位置.parquet')
get_jzdw_location_baidu <- function(
jzdw_df,
baidu_key,
max_retry = 3,
base_timeout = 15
) {
# 步骤1：拼接地址（适配百度地图，格式不变）
jzdw_df <- jzdw_df |>
mutate(
接种单位所属区县名称 = ifelse(
is.na(接种单位所属区县名称) | 接种单位所属区县名称 == "",
"",
接种单位所属区县名称
),
完整地址 = paste0("甘肃省", 接种单位所属区县名称, 接种单位名称) |>
str_squish(),
地址编码 = URLencode(完整地址, reserved = TRUE),
行号 = row_number() # 添加行号，方便跟踪进度
)
# 辅助函数：安全提取字段（确保返回标量字符，而非列表）
safe_extract <- function(x, field) {
val <- tryCatch(
{
if (is.list(x) && field %in% names(x)) {
res <- as.character(x[[field]])
if (length(res) == 0 || res == "" || is.null(res)) NA else res
} else {
NA_character_
}
},
error = function(e) {
NA_character_
}
)
if (length(val) > 1) val[1] else val
}
# 步骤2：单个地址编码（适配百度地图API）
geocode_single <- function(address_encoded, key, row_num, full_address) {
# 1秒调用间隔（遵守百度API频率限制）
Sys.sleep(1)
# 兜底返回值（类型统一）
default <- tibble(
经度 = as.numeric(NA), # 百度返回百度坐标系(BD09)经度
纬度 = as.numeric(NA), # 百度返回百度坐标系(BD09)纬度
匹配地址 = NA_character_,
country = NA_character_,
province = NA_character_,
city = NA_character_,
district = NA_character_,
调用状态 = "失败-未尝试"
)
# 重试机制
for (retry in 1:max_retry) {
tryCatch(
{
# 百度地图地理编码API URL（核心替换）
# 文档：http://lbsyun.baidu.com/index.php?title=webapi/guide/webservice-geocoding
url <- sprintf(
"http://api.map.baidu.com/geocoding/v3/?address=%s&output=json&ak=%s&ret_coordtype=bd09ll",
address_encoded,
key
)
# 发送请求（百度API用http，非https）
resp <- GET(
url,
timeout(base_timeout),
add_headers(
'User-Agent' = 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36'
)
)
# 检查响应状态
if (http_error(resp)) {
message(sprintf(
"[第%s行][重试%s] %s - HTTP错误：%s",
row_num,
retry,
full_address,
http_status(resp)$message
))
next
}
# 解析JSON（适配百度返回格式）
json_text <- content(resp, "text", encoding = "UTF-8")
json_data <- fromJSON(json_text, simplifyVector = TRUE)
# 百度API返回状态码：0=成功，其他=失败
if (is.null(json_data) || json_data$status != 0) {
err_msg <- ifelse(
is.null(json_data$message),
"未知错误",
json_data$message
)
message(sprintf(
"[第%s行][重试%s] %s - API返回失败：%s",
row_num,
retry,
full_address,
err_msg
))
next
}
# 提取百度返回的结果（核心适配：百度的结果在result字段下）
result <- json_data$result
if (length(result) == 0 || !is.list(result)) {
message(sprintf(
"[第%s行][重试%s] %s - 无有效地理编码",
row_num,
retry,
full_address
))
next
}
# 解析百度坐标系经纬度（BD09）
lng <- as.numeric(safe_extract(result$location, "lng")) # 经度
lat <- as.numeric(safe_extract(result$location, "lat")) # 纬度
# 提取行政区域信息（百度的地址组件在addressComponent下）
addr_comp <- result$addressComponent
country <- safe_extract(addr_comp, "country") # 国家
province <- safe_extract(addr_comp, "province") # 省份
city <- safe_extract(addr_comp, "city") # 城市
district <- safe_extract(addr_comp, "district") # 区县
formatted_addr <- safe_extract(result, "formatted_address") # 匹配地址
# 成功返回结果（类型统一）
result_tbl <- tibble(
经度 = lng,
纬度 = lat,
匹配地址 = formatted_addr,
country = country,
province = province,
city = city,
district = district,
调用状态 = "成功"
)
message(sprintf("[第%s行] %s - 处理成功", row_num, full_address))
return(result_tbl)
},
error = function(e) {
# 捕获超时/连接错误
message(sprintf(
"[第%s行][重试%s] %s - 错误：%s",
row_num,
retry,
full_address,
e$message
))
if (retry == max_retry) {
default$调用状态 <- paste0("失败-最终超时/错误：", e$message)
}
}
)
}
# 所有重试都失败，返回兜底值
return(default)
}
# 步骤3：批量处理
message(sprintf(
"开始处理%s条地址，每次调用间隔1秒，最多重试%s次",
nrow(jzdw_df),
max_retry
))
jzdw_with_loc <- jzdw_df |>
rowwise() |>
do({
res <- geocode_single(.$地址编码, baidu_key, .$行号, .$完整地址)
bind_cols(., res)
}) |>
ungroup() |>
select(
行号,
接种单位名称,
接种单位所属区县名称,
完整地址,
经度,
纬度,
匹配地址,
country,
province,
city,
district,
调用状态
)
# 输出失败统计
fail_count <- sum(jzdw_with_loc$调用状态 != "成功")
message(sprintf(
"处理完成！成功：%s条，失败：%s条",
nrow(jzdw_df) - fail_count,
fail_count
))
return(jzdw_with_loc)
}
baidu_api_key <- "XXXXXXXXXXXXX" # 请替换成实际的Key
# 调用函数
jzdw_with_location_baidu <- get_jzdw_location_baidu(
jzdw,
baidu_api_key,
max_retry = 3,
base_timeout = 15
)
jzdw_with_location_baidu |>
write_parquet('./jzdw_with_location_baidu.parquet')
# ==================== 百度坐标转换函数（BD09 → WGS84） ====================
# 辅助函数：弧度转换
rad <- function(d) {
return(d * pi / 180.0)
}
# 步骤1：BD09转GCJ02（百度坐标系转火星坐标系）
bd09_to_gcj02 <- function(bd_lng, bd_lat) {
# 空值处理
if (is.na(bd_lng) || is.na(bd_lat)) {
return(list(lng = NA, lat = NA))
}
x_pi <- pi * 3000.0 / 180.0
pi_val <- pi
x = as.numeric(bd_lng) - 0.0065
y = as.numeric(bd_lat) - 0.006
z = sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi)
theta = atan2(y, x) - 0.000003 * cos(x * x_pi)
gcj_lng = z * cos(theta)
gcj_lat = z * sin(theta)
return(list(lng = gcj_lng, lat = gcj_lat))
}
# 步骤2：GCJ02转WGS84（火星坐标系转WGS84）
gcj02_to_wgs84 <- function(gcj_lng, gcj_lat) {
# 空值处理
if (is.na(gcj_lng) || is.na(gcj_lat)) {
return(list(lng = NA, lat = NA))
}
a = 6378245.0 # 长半轴
ee = 0.00669342162296594323 # 偏心率平方
pi_val = pi
gcj_lng = as.numeric(gcj_lng)
gcj_lat = as.numeric(gcj_lat)
# 超出中国范围直接返回原坐标（简化处理）
if (
gcj_lng < 72.004 ||
gcj_lng > 137.8347 ||
gcj_lat < 0.8293 ||
gcj_lat > 55.8271
) {
return(list(lng = gcj_lng, lat = gcj_lat))
}
dlat = transformlat(gcj_lng - 105.0, gcj_lat - 35.0)
dlng = transformlng(gcj_lng - 105.0, gcj_lat - 35.0)
radlat = gcj_lat / 180.0 * pi_val
magic = sin(radlat)
magic = 1 - ee * magic * magic
sqrtmagic = sqrt(magic)
dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi_val)
dlng = (dlng * 180.0) / (a / sqrtmagic * cos(radlat) * pi_val)
wgs_lat = gcj_lat - dlat
wgs_lng = gcj_lng - dlng
return(list(lng = wgs_lng, lat = wgs_lat))
}
# 辅助函数：纬度转换
transformlat <- function(lng, lat) {
pi_val = pi
ret = -100.0 +
2.0 * lng +
3.0 * lat +
0.2 * lat * lat +
0.1 * lng * lat +
0.2 * sqrt(abs(lng))
ret = ret +
(20.0 * sin(6.0 * lng * pi_val) + 20.0 * sin(2.0 * lng * pi_val)) *
2.0 /
3.0
ret = ret +
(20.0 * sin(lat * pi_val) + 40.0 * sin(lat / 3.0 * pi_val)) * 2.0 / 3.0
ret = ret +
(160.0 * sin(lat / 12.0 * pi_val) + 320 * sin(lat * pi_val / 30.0)) *
2.0 /
3.0
return(ret)
}
# 辅助函数：经度转换
transformlng <- function(lng, lat) {
pi_val = pi
ret = 300.0 +
lng +
2.0 * lat +
0.1 * lng * lng +
0.1 * lng * lat +
0.1 * sqrt(abs(lng))
ret = ret +
(20.0 * sin(6.0 * lng * pi_val) + 20.0 * sin(2.0 * lng * pi_val)) *
2.0 /
3.0
ret = ret +
(20.0 * sin(lng * pi_val) + 40.0 * sin(lng / 3.0 * pi_val)) * 2.0 / 3.0
ret = ret +
(150.0 * sin(lng / 12.0 * pi_val) + 300.0 * sin(lng / 30.0 * pi_val)) *
2.0 /
3.0
return(ret)
}
# 最终封装：BD09直接转WGS84
bd09_to_wgs84 <- function(bd_lng, bd_lat) {
# 先转GCJ02
gcj <- bd09_to_gcj02(bd_lng, bd_lat)
# 再转WGS84
wgs <- gcj02_to_wgs84(gcj$lng, gcj$lat)
return(wgs)
}
# ==================== 对数据集进行坐标转换 ====================
# 假设你的数据集是 jzdw_with_location_baidu
# 新增WGS84坐标列
jzdw_with_location_wgs84 <- jzdw_with_location_baidu %>%
rowwise() %>%
mutate(
# 调用转换函数
wgs84_coords = list(bd09_to_wgs84(经度, 纬度)),
# 提取WGS84经度和纬度
经度_WGS84 = wgs84_coords$lng,
纬度_WGS84 = wgs84_coords$lat
) %>%
ungroup() %>%
# 移除临时列，保留有用字段
select(-wgs84_coords) %>%
# 调整列顺序（把WGS84坐标放在原有坐标后）
select(
行号,
接种单位名称,
接种单位所属区县名称,
完整地址,
经度,
纬度,
经度_WGS84,
纬度_WGS84, # BD09和WGS84坐标对比
匹配地址,
country,
province,
city,
district,
调用状态
)
评论
发表评论