為了讓大家更好的理解本文內(nèi)容，歡迎各位培訓(xùn)班會員參加明晚 8 點的直播課：「使用 Stata 進行地理編碼：地址解析經(jīng)緯度、坐標轉(zhuǎn)換 & 根據(jù)經(jīng)緯度判斷所處的省市區(qū)縣」

由于高德地圖接口的更新，之前的代碼不再適用了，主要是兩個更新：

1. 高德地圖的接口會根據(jù)一個地址返回多個結(jié)果；
2. 高德地圖不再建議一次解析多個地址了，會返回混亂的結(jié)果；
3. Stata 中可以使用 parallel 進行多線程解析，大大提高解析速度；
4. 如果需要進行大量的解析，可以聯(lián)系李老師幫忙，費用參考這個：https://rstata./#/course/6b44ce2712af47b6ba7b2a2fd74a4e68

之前給大家講解過根據(jù)地址解析經(jīng)緯度并根據(jù)經(jīng)緯度判斷所處的省市區(qū)縣的方法，不過那個方法效率較低，面對數(shù)百萬地址的工作就顯得力不從心了。今天再給大家介紹一套無比高效的解析流程。

這套工作流程包含下面幾個過程：

1. 使用高德地圖地理編碼接口解析（含多線程內(nèi)容）；
2. 反復(fù)重復(fù)上述過程，直到篩選出那些無法使用高德地圖地理編碼接口解析成功的地址；
3. 然后使用百度地圖地理編碼接口解析上面的操作沒有成功的部分；
4. 合并所有的解析結(jié)果；
5. GCJ02 坐標轉(zhuǎn)換成 WGS84 坐標；
6. 使用地理計算根據(jù)經(jīng)緯度判斷所處的省市區(qū)縣。

下面就讓我們以金融機構(gòu)網(wǎng)點數(shù)據(jù)的地理編碼為例進行講解。

讀取并處理數(shù)據(jù)

首先我們讀取 xlsx 文件并選擇前 10000 個觀測值：

申請高德地圖和百度地圖的密鑰

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

高德地圖的 city 參數(shù)既可以是城市名稱也可以是城市區(qū)劃代碼，百度地圖的 city 參數(shù)只能是城市名稱，所以我們需要準備兩個 city 變量。

gen address = 機構(gòu)地址
gen citygd = 地址代碼 
gen citybd = 機構(gòu)所在地 
format address city* %20s
replace citybd = ustrregexs(1) if ustrregexm(citybd, '-(.*)')
keep id 機構(gòu)名稱 address citybd citygd

*- 地址里面不能出現(xiàn) # 符號，替換掉
replace address = subinstr(address, '#', '號', .)
replace address = subinstr(address, ' ', '', .)
save dfsim, replace 

使用高德地圖解析

在這套工作流程中我們會優(yōu)先使用高德地圖地理編碼接口進行解析，因為高德地圖支持更高的 QPS（每秒請求的次數(shù)）。不過不幸運的是現(xiàn)在的高德和百度都只提供每日 5000 次的免費額度，難以進行大量的解析。

這里插播一條廣告。如果大家需要大量的解析，可以考慮聯(lián)系李老師幫忙解析。費用可以參考這個：https://rstata./#/course/6b44ce2712af47b6ba7b2a2fd74a4e68

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

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

*- 解析 json 文件
*- 查看 json 文件的結(jié)構(gòu) 
*- 安裝
*- 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 = 北京市豐臺區(qū)中國農(nóng)業(yè)發(fā)展銀行(北京市分行)
*>         geocodes:1:country = 中國
*>         geocodes:1:province = 北京市
*>         geocodes:1:citycode = 010
*>         geocodes:1:city = 北京市
*>         geocodes:1:district = 豐臺區(qū)
*>         geocodes:1:adcode = 110106
*>         geocodes:1:location = 116.329972,39.865275
*>         geocodes:1:level = 興趣點
*>         geocodes:2:formatted_address = 北京市豐臺區(qū)南四環(huán)西路186號院
*>         geocodes:2:country = 中國
*>         geocodes:2:province = 北京市
*>         geocodes:2:citycode = 010
*>         geocodes:2:city = 北京市
*>         geocodes:2:district = 豐臺區(qū)
*>         geocodes:2:adcode = 110106
*>         geocodes:2:street = 南四環(huán)西路
*>         geocodes:2:number = 186號院
*>         geocodes:2:location = 116.295491,39.825012
*>         geocodes:2:level = 門址
*>         geocodes:3:formatted_address = 北京市豐臺區(qū)南四環(huán)西路186號
*>         geocodes:3:country = 中國
*>         geocodes:3:province = 北京市
*>         geocodes:3:citycode = 010
*>         geocodes:3:city = 北京市
*>         geocodes:3:district = 豐臺區(qū)
*>         geocodes:3:adcode = 110106
*>         geocodes:3:street = 南四環(huán)西路
*>         geocodes:3:location = 116.283150,39.823781
*>         geocodes:3:level = 門牌號

由此我們可以編寫下面的代碼處理 json 文件：

這里得到了 3 個結(jié)果，通常我們會選擇第一條作為最優(yōu)結(jié)果，另外也可以根據(jù) level 篩選有較大可能性的結(jié)果。

處理策略可以是先把所有的 json 文件下載下來，然后再處理：

use dfsim, clear 

*- 創(chuàng)建文件夾保存 json 文件
*- 如果數(shù)據(jù)里面沒有城市變量的話，這里去掉 '&city=...' 參數(shù)就可以了
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：

運行的過程中如果出現(xiàn)了網(wǎng)絡(luò)中斷等問題，可以加個 if 語句重新運行上面的代碼。運行結(jié)束之后也可以檢查 json1 文件夾里面的文件，刪去錯誤的結(jié)果再重新運行。

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 來提高下載速度：

另外也有一種多線程的本方法，就是把循環(huán)拆分成幾部分，然后打開多個 Stata 分別運行。Windows 電腦上可以對著 Stata 的圖標右鍵選擇打開新的 Stata，Mac 可以從終端打開多個 Stata 的命令行窗口。

使用 insheetjson 處理 json 文件

然后循環(huán)處理下這些 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 = ''
*- 循環(huán) 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 

和原始數(shù)據(jù)匹配：

然后就可以進行篩選了，去除不符合要求的觀測值：

*- 刪去城市的前兩個字不一致的（如果沒有城市變量就算了）
drop if substr(citybd, 1, 6) != substr(city, 1, 6)

*- 刪除 level 變量為 省市區(qū)縣的，這說明結(jié)果范圍過于寬泛 
drop if inlist(level, '省', '市', '區(qū)縣')

*- 其他大家自己覺得不合理的也可以進行刪除。

*- 最后我們就可以解決結(jié)果重復(fù)的問題了，比較常見的一種做法就是保留每組的第一個觀測值
bysort id: keep if _n == 1 

*- 根據(jù)經(jīng)緯度可以判斷所處的省市區(qū)縣，所以這里先刪去省市區(qū)縣信息 
keep id location 
save '解析結(jié)果1', replace 

沒有成功的：

使用百度地圖地理編碼接口解析

剩下的這些使用百度地圖地理編碼接口解析：

use dfsim, clear 
merge 1:1 id using 解析結(jié)果1 
keep if _m == 1 
drop _m 
replace address = address + 機構(gòu)名稱 
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 = ''

*- 循環(huán) 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

*- 這里也可以根據(jù)數(shù)據(jù)中的變量判斷解析情況。
keep id location 
save '解析結(jié)果2', replace 

合并所有的成功結(jié)果

GCJ02 坐標轉(zhuǎn)換成 WGS84 坐標

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

gcj02towgs84, lon(經(jīng)度) lat(緯度)
save finalresult, replace 

根據(jù)經(jīng)緯度判斷所處的省市區(qū)縣

之前介紹過這個操作，不過不是 Stata 完成的，最近發(fā)現(xiàn)原來 Stata 中也有這種命令：geoinpoly。

安裝：

把 shp 文件轉(zhuǎn)換成 dta 文件：

shp2dta using 2021行政區(qū)劃/縣.shp, database(county_db) coordinates(county_coord) genid(ID) replace 

然后就可以使用該命令進行經(jīng)緯度判斷省市區(qū)縣的操作了：

效率非常高！

直播信息

為了讓大家更好的理解本文內(nèi)容，歡迎各位培訓(xùn)班會員參加明晚 8 點的直播課：「使用 Stata 進行地理編碼：地址解析經(jīng)緯度、坐標轉(zhuǎn)換 & 根據(jù)經(jīng)緯度判斷所處的省市區(qū)縣」

1. 直播地址：騰訊會議(需要報名 RStata 培訓(xùn)班參加)
2. 講義材料：需要購買 RStata 名師講堂會員，詳情可閱讀：一起來學(xué)習(xí) R 語言和 Stata 啦！學(xué)習(xí)過程中遇到的問題也可以隨時提問！