使用 Stata 进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县

由于高德地图接口的更新，之前的代码不再适用了，主要是两个更新：

1. 高德地图的接口会根据一个地址返回多个结果；
2. 高德地图不再建议一次解析多个地址了，会返回混乱的结果；
3. Stata 中可以使用 parallel 进行多线程解析，大大提高解析速度；
4. 如果需要进行大量的解析，可以联系李老师帮忙，费用参考这个：https://rstata.duanshu.com/#/course/6b44ce2712af47b6ba7b2a2fd74a4e68

之前给大家讲解过根据地址解析经纬度并根据经纬度判断所处的省市区县的方法，不过那个方法效率较低，面对数百万地址的工作就显得力不从心了。今天再给大家介绍一套无比高效的解析流程。

这套工作流程包含下面几个过程：

1. 使用高德地图地理编码接口解析（含多线程内容）；
2. 反复重复上述过程，直到筛选出那些无法使用高德地图地理编码接口解析成功的地址；
3. 然后使用百度地图地理编码接口解析上面的操作没有成功的部分；
4. 合并所有的解析结果；
5. GCJ02 坐标转换成 WGS84 坐标；
6. 使用地理计算根据经纬度判断所处的省市区县。

下面就让我们以金融机构网点数据的地理编码为例进行讲解。

读取并处理数据

首先我们读取 xlsx 文件并选择前 10000 个观测值：

*- 设定工作目录
cd "~/Desktop/使用Stata进行地理编码：高德与百度接口/"
import excel using "金融机构网点数据.xlsx", clear first
foreach i of varlist _all {
 cap format `i' %10s
}

*- 选择前 1000 个
keep in 1/1000
save "金融机构网点数据", replace

*- 生成 id 变量来标识观测值
use 金融机构网点数据, clear
gen id = _n
order id

*- 选择感兴趣的变量
keep id 机构名称 机构所在地 机构地址 地址代码

申请高德地图和百度地图的密钥

1. 高德地图文档：https://lbs.amap.com/api/webservice/guide/api/georegeo
2. 百度地图文档：https://lbsyun.baidu.com/index.php?title=webapi/guide/webservice-geocoding-base

高德地图的 city 参数既可以是城市名称也可以是城市区划代码，百度地图的 city 参数只能是城市名称，所以我们需要准备两个 city 变量。

gen address = 机构地址
gen citygd = 地址代码
gen citybd = 机构所在地
format address city* %20s
replace citybd = ustrregexs(1) if ustrregexm(citybd, "-(.*)")
keep id 机构名称 address citybd citygd

*- 地址里面不能出现 # 符号，替换掉
replace address = subinstr(address, "#", "号", .)
replace address = subinstr(address, " ", "", .)
save dfsim, replace

使用高德地图解析

在这套工作流程中我们会优先使用高德地图地理编码接口进行解析，因为高德地图支持更高的 QPS（每秒请求的次数）。不过不幸运的是现在的高德和百度都只提供每日 5000 次的免费额度，难以进行大量的解析。

这里插播一条广告。如果大家需要大量的解析，可以考虑联系李老师帮忙解析。费用可以参考这个：https://rstata.duanshu.com/#/course/6b44ce2712af47b6ba7b2a2fd74a4e68

例如解析北京的十条地址：

use dfsim, clear
keep if citygd == "1100"
local address = ""
forval i = 1/10{
 local address = "`address'|`=address[`i']'"
}
di "`address'"
local address = substr("`address'", 2, .)
di "`address'"

*- URL 转码: percentencode.ado
percentencode `address'
*> %E5%8C%97%E4%BA%AC%E5%B8%82%E4%B8%B0%E5%8F%B0%E5%8C%BA%E5%8D%97%E5%9B%9B%E7%8E%AF%E8%A5%BF%E8%B7%AF
*> > 186%E5%8F%B7%E4%B8%80%E5%8C%BA1%E5%8F%B7%E6%A5%BC5%E5%B1%82%E4%B8%AD%E5%9B%BD%E5%86%9C%E4%B8%9A%E
*> > 5%8F%91%E5%B1%95%E9%93%B6%E8%A1%8C%E5%8C%97%E4%BA%AC%E5%B8%82%E5%88%86%E8%A1%8C

ret list
*> macros:
*>       r(percentencode) : "%E5%8C%97%E4%BA%AC%E5%B8%82%E4%B8%B0%E5%8F%B0%E5%8C%BA%E5%8D%97%E5%9B.."

*- 构造 URL
local url = "https://restapi.amap.com/v3/geocode/geo?address=" + ///
 "`r(percentencode)'" + ///
 "&key=091b12221c04bf3049e0657d9d5dea0a&batch=true&city=" + ///
 "1100"

di "`url'"
*> https://restapi.amap.com/v3/geocode/geo?address=%E5%8C%97%E4%BA%AC%E5%B8%82%E4%B8%B0%E5%8F%B0%E5%8
*> > C%BA%E5%8D%97%E5%9B%9B%E7%8E%AF%E8%A5%BF%E8%B7%AF186%E5%8F%B7%E4%B8%80%E5%8C%BA1%E5%8F%B7%E6%A5%
*> > BC5%E5%B1%82%E4%B8%AD%E5%9B%BD%E5%86%9C%E4%B8%9A%E5%8F%91%E5%B1%95%E9%93%B6%E8%A1%8C%E5%8C%97%E4
*> > %BA%AC%E5%B8%82%E5%88%86%E8%A1%8C&key=50828d53749c02431a65192f7c092a77&city=1100

*- 下载 json 文件
copy "`url'" temp.json, replace

然后我们就可以使用 insheetjson 命令处理得到的 json 文件了：

*- 解析 json 文件
*- 查看 json 文件的结构
*- 安装
*- ssc install insheetjson
*- ssc install libjson
insheetjson using "temp.json", showr flatten

*> Response from server:
*>         status = 1
*>         info = OK
*>         infocode = 10000
*>         count = 3
*>         geocodes:1:formatted_address = 北京市丰台区中国农业发展银行(北京市分行)
*>         geocodes:1:country = 中国
*>         geocodes:1:province = 北京市
*>         geocodes:1:citycode = 010
*>         geocodes:1:city = 北京市
*>         geocodes:1:district = 丰台区
*>         geocodes:1:adcode = 110106
*>         geocodes:1:location = 116.329972,39.865275
*>         geocodes:1:level = 兴趣点
*>         geocodes:2:formatted_address = 北京市丰台区南四环西路186号院
*>         geocodes:2:country = 中国
*>         geocodes:2:province = 北京市
*>         geocodes:2:citycode = 010
*>         geocodes:2:city = 北京市
*>         geocodes:2:district = 丰台区
*>         geocodes:2:adcode = 110106
*>         geocodes:2:street = 南四环西路
*>         geocodes:2:number = 186号院
*>         geocodes:2:location = 116.295491,39.825012
*>         geocodes:2:level = 门址
*>         geocodes:3:formatted_address = 北京市丰台区南四环西路186号
*>         geocodes:3:country = 中国
*>         geocodes:3:province = 北京市
*>         geocodes:3:citycode = 010
*>         geocodes:3:city = 北京市
*>         geocodes:3:district = 丰台区
*>         geocodes:3:adcode = 110106
*>         geocodes:3:street = 南四环西路
*>         geocodes:3:location = 116.283150,39.823781
*>         geocodes:3:level = 门牌号

由此我们可以编写下面的代码处理 json 文件：

*- 准备空变量
clear
gen str100 formatted_address = ""
gen str100 province = ""
gen str100 citycode = ""
gen str100 city = ""
gen str100 district = ""
gen str100 adcode = ""
gen str100 location = ""
gen str100 level = ""
*- 处理 json 数据
insheetjson formatted_address province citycode city district adcode location level using "temp.json", columns("formatted_address" "province" "citycode" "city" "district" "adcode" "location" "level") tableselector("geocodes")
compress

*>     +-----------------------------------------------------------------+
*>     |                        formatted_address               location |
*>     |-----------------------------------------------------------------|
*>  1. | 北京市丰台区中国农业发展银行(北京市分行)   116.329972,39.865275 |
*>  2. |            北京市丰台区南四环西路186号院   116.295491,39.825012 |
*>  3. |              北京市丰台区南四环西路186号   116.283150,39.823781 |
*>     +-----------------------------------------------------------------+

这里得到了 3 个结果，通常我们会选择第一条作为最优结果，另外也可以根据 level 筛选有较大可能性的结果。

处理策略可以是先把所有的 json 文件下载下来，然后再处理：

use dfsim, clear

*- 创建文件夹保存 json 文件
*- 如果数据里面没有城市变量的话，这里去掉 "&city=..." 参数就可以了
cap mkdir "json1"
forval i = 1/`=_N' {
 local address = "`=address[`i']'"
 percentencode `address'
 local url = "https://restapi.amap.com/v3/geocode/geo?address=" + ///
  "`r(percentencode)'" + ///
  "&key=50828d53749c02431a65192f7c092a77&city=" + ///
  "`=citygd[`i']'"
 copy "`url'" "json1/`=id[`i']'.json", replace
}

或者也可以使用 curl 代替 copy：

*- 或者也可以使用 curl 代替 copy：
*- !curl "`url'" -o "json1/`=id[`i']'.json"

运行的过程中如果出现了网络中断等问题，可以加个 if 语句重新运行上面的代码。运行结束之后也可以检查 json1 文件夹里面的文件，删去错误的结果再重新运行。

forval i = 1/`=_N' {
 if !fileexists("json1/`=id[`i']'.json") {
  local address = "`=address[`i']'"
  percentencode `address'
  local url = "https://restapi.amap.com/v3/geocode/geo?address=" + ///
   "`r(percentencode)'" + ///
   "&key=50828d53749c02431a65192f7c092a77&city=" + ///
   "`=citygd[`i']'"
  copy "`url'" "json1/`=id[`i']'.json", replace
 }
}

使用并行提高速度

Stata 中可以使用 parallel 命令进行并行运算，这里也可以使用 parallel 来提高下载速度：

*- 安装：ssc install parallel

use dfsim, clear
*- 设置进程数量
parallel initialize 8

*- 把循环部分的代码保存成 forval.do
parallel do forval.do

*- 这样效率会快很多！很快就会发现 json1 文件夹中的文件已经下载完了。

另外也有一种多线程的本方法，就是把循环拆分成几部分，然后打开多个 Stata 分别运行。Windows 电脑上可以对着 Stata 的图标右键选择打开新的 Stata，Mac 可以从终端打开多个 Stata 的命令行窗口。

使用 insheetjson 处理 json 文件

然后循环处理下这些 json 文件，运行前先可以手动检查下有没有错误的 json 文件，然后手动删除掉：

clear
gen str100 formatted_address = ""
gen str100 province = ""
gen str100 city = ""
gen str100 district = ""
gen str100 location = ""
gen str100 level = ""
gen str100 id = ""
*- 循环 json1 文件夹下的所有文件
local files: dir "json1" files "*.json"
local n: word count `files'

foreach i in `files' {
 insheetjson formatted_address province city ///
  district location level using "json1/`i'", ///
  columns("formatted_address" "province" ///
   "city" "district" "location" "level") ///
  tableselector("geocodes") offset(`=_N')
 replace id = "`i'" if mi(id)
}

compress

foreach i of varlist _all {
  cap format `i' %10s
}

order id
replace id = subinstr(id, ".json", "", .)
destring id, replace

和原始数据匹配：

merge m:1 id using dfsim
keep if _m == 3
drop _m
order id 机构名称 address citygd citybd

然后就可以进行筛选了，去除不符合要求的观测值：

*- 删去城市的前两个字不一致的（如果没有城市变量就算了）
drop if substr(citybd, 1, 6) != substr(city, 1, 6)

*- 删除 level 变量为 省市区县的，这说明结果范围过于宽泛
drop if inlist(level, "省", "市", "区县")

*- 其他大家自己觉得不合理的也可以进行删除。

*- 最后我们就可以解决结果重复的问题了，比较常见的一种做法就是保留每组的第一个观测值
bysort id: keep if _n == 1

*- 根据经纬度可以判断所处的省市区县，所以这里先删去省市区县信息
keep id location
save "解析结果1", replace

没有成功的：

use dfsim, clear
merge 1:1 id using 解析结果1
keep if _m == 1
drop _m

使用百度地图地理编码接口解析

剩下的这些使用百度地图地理编码接口解析：

use dfsim, clear
merge 1:1 id using 解析结果1
keep if _m == 1
drop _m
replace address = address + 机构名称
drop location
compress
save dfsim2, replace

use dfsim2, clear
cap mkdir "json2"
forval g = 1/`=_N' {
 percentencode `=address[`g']'
 local add = "`r(percentencode)'"
 percentencode `=citybd'
 local cty = "`r(percentencode)'"
 local url = "https://api.map.baidu.com/geocoding/v3/?address=`add'&output=json&ak=h4jN43d3ajnnXGzBHmIowERWwM2o6PIe&ret_coordtype=gcj02ll&city=`cty'"
 copy "`url'" "json2/`=id[`g']'.json"
}

*- 看起来都成功了
clear
gen str100 lng = ""
gen str100 lat = ""
gen str100 precise = ""
gen str100 confidence = ""
gen str100 comprehension = ""
gen str100 level = ""
gen str100 id = ""

*- 循环 json2 文件夹下的所有文件
local files: dir "json2" files "*.json"
local n: word count `files'
foreach i in `files' {
 insheetjson lng lat precise confidence ///
  comprehension level ///
  using "json2/`i'", table("result") ///
  columns("location:lng" "location:lat" "precise" ///
   "confidence" "comprehension" "level") offset(`=_N')
 replace id = "`i'" if mi(id)
}
compress
order id
replace id = subinstr(id, ".json", "", .)
destring id, replace
gen location = lng + "," + lat

*- 这里也可以根据数据中的变量判断解析情况。
keep id location
save "解析结果2", replace

合并所有的成功结果

use 解析结果1, clear
append using 解析结果2
gsort id
unique id
split location, parse(,)
drop location
ren location1 经度
ren location2 纬度
destring 经度 纬度, replace

GCJ02 坐标转换成 WGS84 坐标

附件中的 gcj02towgs84.ado 命令可以直接完成这个操作：

gcj02towgs84, lon(经度) lat(纬度)
save finalresult, replace

根据经纬度判断所处的省市区县

之前介绍过这个操作，不过不是 Stata 完成的，最近发现原来 Stata 中也有这种命令：geoinpoly。

安装：

ssc install geoinpoly

把 shp 文件转换成 dta 文件：

shp2dta using 2021行政区划/县.shp, database(county_db) coordinates(county_coord) genid(ID) replace

然后就可以使用该命令进行经纬度判断省市区县的操作了：

use finalresult, clear
geoinpoly 纬度_WGS84 经度_WGS84 using county_coord.dta
ren _ID ID
merge m:1 ID using county_db
keep if _m == 3
drop _m
foreach i of varlist _all {
 cap format `i' %10s
}
*- 这样我们就完成了判断所处省市区县的操作
save finalresult2, replace

效率非常高！

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县