其他
使用 Stata 爬取全国宗庙祠堂数据并绘图展示
由于之前的思路不太好，这里又重新更新了下～
该课程已经讲解完，感兴趣的小伙伴可以点击文末的阅读原文跳转到平台上观看学习！视频讲解中使用的是 utm 坐标系，这里改用了 aea 坐标系，这样由于坐标系转换产生的距离变化会更小，课程里面的检索结果更多是因为 utm 转换会 wgs84 后的格点实际上更细密了。
最近有个小伙伴想要获取全国宗庙祠堂数据，还附带给我转发了一篇推文：绘制全国祠堂密度地图。
这篇推文是使用 Stata 爬取的，具体思路是先计算各区县的经纬度，然后使用百度地图的“地点检索”接口搜索附近 50km 范围内的宗庙祠堂。不过显然这种方法会遗漏很多，主要是：
这个范围会遗漏中国的很多地区； 50km 的范围太大，百度地图只会给出该范围的一些热门的搜索结果。
因此我们需要生成更细密的检索中心，然后使用更小的搜索半径。
百度地图的地点检索接口文档在这里：https://lbsyun.baidu.com/faq/api?title=webapi/guide/webservice-placeapi/circle ，虽然多边形检索看起来更高效，不过这个接口需要单独申请，所以我们还是只能使用圆形区域检索。
在使用地点检索之前，我们需要先构建一个经纬度中心点数据。
为了便于演示，在这里我们依然使用 50km 的半径进行检索，由于检索是圆形检索，所以为了不遗漏，我们需要建立分辨率与 50√2km x 50√2km 的经纬度网格，这样以每个网格的中心为检索中心、范围是 50km 的圆形区域恰好可以覆盖这一区域，如下图所示：
由于 Stata 中地理计算的工具不完善，所以我们需要借助 R 语言来生成这样的一个网格。
不会使用 R 语言的小伙伴也不必担心，附件中提供了最终的结果数据。
首先加载所需的 R 包：
library(tidyverse)
library(sf)
合并各省份的数据得到中国范围的矢量数据：
# 生成中国范围的矢量数据
read_sf("2020行政区划/省.shp") %>%
st_union() -> cn
nngeo::st_remove_holes(cn[2]) -> cn
这个中国的范围是这样的：
url <- "http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}"
library(leaflet)
leaflet() %>%
addTiles(url, attribution = "微信公众号 RStata") -> map
mapview::mapview(cn, map = map, layer.name = "中国范围")
创建一些边长为 50√2 km的格网，这样每个方格正好可以使用一个半径 50km 的圆覆盖。转换成以 m 为单位的坐标系更容易创建：
# 这里以讲义材料的这个代码为准
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"
cn %>%
st_transform(mycrs) -> cnaea
cnaea %>%
st_make_grid(cellsize = units::set_units(c(50 * sqrt(2), 50 * sqrt(2)), "km")) -> cngrid
这个网格的效果是这样的：
# 效果是这样的
read_sf("九段线.geojson") %>%
st_transform(mycrs) -> jdxaea
plot(cngrid)
plot(cnaea, add = T, col = "red")
plot(jdxaea, add = T, col = "black")
提取与中国相交的格点：
cngrid %>%
st_intersection(cnaea) -> cngrid
plot(cngrid)
计算这些格点的质心：
cngrid %>%
st_sf() %>%
st_centroid() %>%
st_transform(4326) %>%
st_coordinates() %>%
as_tibble() %>%
mutate(id = row_number()) -> coorddf
# 保存
coorddf %>%
haven::write_dta("50km 网格中心.dta")
这样我们就得到了检索中心坐标。另外附件中还提供了 10km 分辨率的网格中心数据。
再回到 Stata 中。
首先我们看一下百度地图地理检索接口如何使用，例如搜索 (39.915,116.404) 附近 2km 的 宗祠、家庙 和 祠堂 可以使用下面的链接进行请求：
clear all
*- "a8niO5h8icau3IRNG2Wv2k5CduVyrktL" 是我的百度地图密钥，需要替换成自己的
copy "https://api.map.baidu.com/place/v2/search?query=宗祠$家庙$祠堂&location=39.915,116.404&radius=2000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=0" "temp.json", replace
这里使用了这些参数：
query=宗祠祠堂（检索关键词） location=39.915,116.404（检索中心） radius=2000（检索半径 2km） output=json（返回结果为 json 格式） ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL（用户密钥，注意要替换成自己的） coord_tyle=wgs84ll（检索中心经纬度的坐标系） ret_coordtype=gcj02ll（返回结果中经纬度的坐标系） page_size=20（每次请求的最大数量，最大只能为 20） page_num=0（0 表示第一页，如果检索结果超过 20 ，需要循环多页爬取）。
有时候也会需要对中文进行 URL 转码，附件中提供了一个 perccentencode.ado 命令，可以进行 URL 转码：
*- 也可能需要对 query 部分进行 URL 转码
percentencode 宗祠
local a = r(percentencode)
percentencode 家庙
local b = r(percentencode)
percentencode 祠堂
local c = r(percentencode)
copy "https://api.map.baidu.com/place/v2/search?query=`a'$`b'$`c'&location=39.915,116.404&radius=2000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=0" "temp.json", replace
这是个 json 文件，在课程「Stata 网络数据爬取：JSON篇」(https://rstata.duanshu.com/#/brief/course/c6de9fae65df4eeb814d2275545e224d) 中我系统的介绍过如何在 Stata 中处理 json 格式的数据，这里我们可以这样解析：
*- 查看响应
insheetjson using "temp.json", showr flatten
*- 解析数据
clear
gen str100 name = ""
gen str100 lat = ""
gen str100 lng = ""
gen str100 address = ""
gen str100 province = ""
gen str100 city = ""
gen str100 area = ""
gen str100 street_id = ""
gen str100 detail = ""
gen str100 uid = ""
insheetjson name lat lng address province city ///
area street_id detail uid using "temp.json", ///
table("results") col("name" "location:lat" ///
"location:lng" "address" "province" "city" ///
"area" "street_id" "detail" "uid")
compress
destring, replace
这里仅仅处理了 page_num=0 页，如果结果超过 20（total）的话，就需要爬取多页了：
注意看接口文档中的 total 介绍：POI 检索总数，开发者请求中设置了 page_num 字段才会出现 total 字段。出于数据保护目的，单次请求 total 最多为150。
*- 第 2 页（如果有的话）
copy "https://api.map.baidu.com/place/v2/search?query=宗祠$家庙$祠堂&location=39.915,116.404&radius=2000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=1" "temp2.json", replace
insheetjson name lat lng address ///
province city area detail uid ///
using temp2.json, table(results) ///
col("name" "location:lat" "location:lng" "address" "province" ///
"city" "area" "detail" "uid") offset(20)
*- 第 3 页（如果有的话）
copy "https://api.map.baidu.com/place/v2/search?query=宗祠$家庙$祠堂&location=39.915,116.404&radius=2000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=2" "temp3.json", replace
insheetjson name lat lng address ///
province city area detail uid ///
using temp3.json, table(results) ///
col("name" "location:lat" "location:lng" "address" "province" ///
"city" "area" "detail" "uid") offset(20)
然后就可以循环爬取处理了。
对于复杂的网络数据爬取任务，建议先下载、下载完成后再集中解析。
首先循环下载首页 json 文件到工作目录下的 res 文件夹：
use "50km 网格中心.dta", clear
*- 循环每个观测值
cap mkdir "res"
*- 爬取首页
forval i = 1/`=_N' {
di "`i'"
if !fileexists("res/`=id[`i']'_0.json") {
copy "https://api.map.baidu.com/place/v2/search?query=宗祠$家庙$祠堂&location=`=Y[`i']',`=X[`i']'&radius=50000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=0" "res/`=id[`i']'_0.json", replace
}
}
从这些首页 json 文件里面提取 total 值：
*- 提取 total 结果
clear all
gen str10 total = ""
insheetjson total using "temp.json", table("total") offset(0)
use "50km 网格中心.dta", clear
tostring id, replace
gen file = "res/" + id + "_0.json"
gen str10 total = ""
forval i = 1/`=_N' {
insheetjson total using "`=file[`i']'", table("total") offset(`=`i' - 1')
}
destring, replace
*- total <= 20 的都是只有一页的
drop if total <= 20
gen pagenum = int(total/20)
*- 循环下载其余页面
forval i = 1/`=_N' {
forval j = 1/`=pagenum[`i']' {
if !fileexists("res/`=id[`i']'_`j'.json") {
copy "https://api.map.baidu.com/place/v2/search?query=宗祠$家庙$祠堂&location=`=Y[`i']',`=X[`i']'&radius=50000&output=json&ak=a8niO5h8icau3IRNG2Wv2k5CduVyrktL&coord_tyle=wgs84ll&ret_coordtype=gcj02ll&page_size=20&page_num=`j'" "res/`=id[`i']'_`j'.json", replace
}
}
}
循环处理下载得到的文件：
local files: dir "res" files "*.json"
di `"`files'"'
clear
gen str100 name = ""
gen str100 lat = ""
gen str100 lng = ""
gen str100 address = ""
gen str100 province = ""
gen str100 city = ""
gen str100 area = ""
gen str100 detail = ""
gen str100 uid = ""
foreach i in `files' {
local N = `=_N'
di "`i'"
insheetjson name lat lng address province city ///
area detail uid using "res/`i'", ///
table("results") col("name" "location:lat" ///
"location:lng" "address" "province" "city" ///
"area" "detail" "uid") offset(`N')
}
compress
destring, replace
foreach i of varlist _all {
cap format `i' %10s
}
*- 删除重复的
duplicates drop _all, force
另外因为得到结果里面的经纬度是 GCJ02 坐标系的，为了后续的地理计算，我们需要把它转换成 WGS84 坐标系的，附件中提供了一个转换命令：
gcj02towgs84, lon(lng) lat(lat)
save "全国宗庙祠堂分布(50km范围的检索)", replace
50km 的检索大概只得到了一万多个宗庙祠堂（实际上每次解析也可能得到不同的结果），在附件中我还提供了一份更全的结果，是使用 10km 范围进行检索的（使用的网格也是 10√2km 的）。下面我们绘图比较下两个结果：
*- 绘制地图（使用更全的一个版本）
*- 这个里面的 lat lng 是我转换成 WGS84 过的，不需要再转换了。
import excel using "全国宗庙祠堂分布（10km 范围的检索）.xlsx", clear first
foreach i of varlist _all {
cap format `i' %10s
}
gen id = _n
encode province, gen(prov)
*- ssc install geo2xy
ren lat 纬度
ren lng 经度
geo2xy 纬度 经度, gen(lat lon) projection(albers, 6378137 298.257223563 25 47 0 105) replace
gen class = "main"
save mainpoint, replace
*- 提取小地图框格的
use mainpoint, clear
keep if inrange(lon, 120000, 1766004.1) & inrange(lat, 320000, 2557786.0)
replace lon = lon * 0.5 + 2100000
replace lat = lat * 0.5 + 1665139
replace class = "smallbox"
append using mainpoint
save pointdata, replace
use chinacity2020mini_db.dta, clear
drop if 省代码 == 0
spmap using chinacity2020mini_coord.dta, id(ID) ///
ocolor("black" ...) osize(vvthin ...) ///
line(data(chinacity2020mini_line_coord2.dta) ///
/// 胡焕庸线（7）秦岭淮河下（4）
select(keep if inlist(group, 1, 2, 3, 5, 6)) ///
by(group) size(vvthin *1 *0.5 *0.5 *0.5) ///
pattern(solid ...) ///
color(white /// 市界颜色
black /// 国界线颜色
"0 85 170" /// 海岸线颜色
black /// 小地图框格颜色
black /// 比例尺和指北针颜色
)) ///
polygon(data(polygon2) fcolor(black) ///
osize(vvthin)) ///
label(data(chinacity2020mini_label2) x(X) y(Y) label(cname) length(20) size(*0.8)) ///
point(data(pointdata) by(prov) ///
fcolor("80 80 255" "206 61 50" "116 155 88" ///
"240 230 133" "70 105 131" "186 99 56" "93 177 221" ///
"128 34 104" "107 215 107" "213 149 167" "146 72 34" ///
"131 123 141" "199 81 39" "213 143 92" "122 101 165" ///
"228 175 105" "59 27 83" "205 222 183" "97 42 121" ///
"174 31 99" "231 199 111" "90 101 94" "204 153 0" ///
"153 204 0" "169 169 169" "204 153 0" "153 204 0" ///
"51 204 0" "0 204 51" "0 204 153" "0 153 204" ///
"10 71 255" "71 117 255" "255 194 10") ///
x(lon) y(lat) ///
size(*0.1 ...) legenda(off)) ///
leg(off) ///
ti("全国宗庙祠堂分布（10km 范围的检索）", color(black)) ///
subti("数据处理 & 绘制：微信公众号 RStata") ///
graphr(margin(medium)) ///
caption("数据来源：百度地图地点检索接口", size(*0.8))
gr export "全国宗庙祠堂分布（10km 范围的检索）.png", replace width(4800)
*- 50km 的
use "全国宗庙祠堂分布(50km范围的检索).dta", clear
gen id = _n
encode province, gen(prov)
*- ssc install geo2xy
drop 纬度 经度
ren 纬度_WGS84 纬度
ren 经度_WGS84 经度
geo2xy 纬度 经度, gen(lat lon) projection(albers, 6378137 298.257223563 25 47 0 105) replace
gen class = "main"
save mainpoint, replace
*- 提取小地图框格的
use mainpoint, clear
keep if inrange(lon, 120000, 1766004.1) & inrange(lat, 320000, 2557786.0)
replace lon = lon * 0.5 + 2100000
replace lat = lat * 0.5 + 1665139
replace class = "smallbox"
append using mainpoint
save pointdata, replace
use chinacity2020mini_db.dta, clear
drop if 省代码 == 0
spmap using chinacity2020mini_coord.dta, id(ID) ///
ocolor("black" ...) osize(vvthin ...) ///
line(data(chinacity2020mini_line_coord2.dta) ///
/// 胡焕庸线（7）秦岭淮河下（4）
select(keep if inlist(group, 1, 2, 3, 5, 6)) ///
by(group) size(vvthin *1 *0.5 *0.5 *0.5) ///
pattern(solid ...) ///
color(white /// 市界颜色
black /// 国界线颜色
"0 85 170" /// 海岸线颜色
black /// 小地图框格颜色
black /// 比例尺和指北针颜色
)) ///
polygon(data(polygon2) fcolor(black) ///
osize(vvthin)) ///
label(data(chinacity2020mini_label2) x(X) y(Y) label(cname) length(20) size(*0.8)) ///
point(data(pointdata) by(prov) ///
fcolor("80 80 255" "206 61 50" "116 155 88" ///
"240 230 133" "70 105 131" "186 99 56" "93 177 221" ///
"128 34 104" "107 215 107" "213 149 167" "146 72 34" ///
"131 123 141" "199 81 39" "213 143 92" "122 101 165" ///
"228 175 105" "59 27 83" "205 222 183" "97 42 121" ///
"174 31 99" "231 199 111" "90 101 94" "204 153 0" ///
"153 204 0" "169 169 169" "204 153 0" "153 204 0" ///
"51 204 0" "0 204 51" "0 204 153" "0 153 204" ///
"10 71 255" "71 117 255" "255 194 10") ///
x(lon) y(lat) ///
size(*0.1 ...) legenda(off)) ///
leg(off) ///
ti("全国宗庙祠堂分布（50km 范围的检索）", color(black)) ///
subti("数据处理 & 绘制：微信公众号 RStata") ///
graphr(margin(medium)) ///
caption("数据来源：百度地图地点检索接口", size(*0.8))
gr export "全国宗庙祠堂分布（50km 范围的检索）.png", replace width(4800)
可以看到两个结果虽然数量上有差异，但是地理分布趋势相近。因此我猜测使用不全的数据结果很可能也可能得到靠谱的结果。所以如果只是准备采用这个数据作为某个变量的代理变量，不用太担心全不全的事情。
最后我们再计算下各个城市的宗庙祠堂密度，首先使用 R 语言计算各个城市的行政区划面积：
# 计算各市的面积
read_sf("2020行政区划/市.shp") -> city
city %>%
mutate(面积 = st_area(.)) %>%
st_drop_geometry() %>%
mutate(面积 = units::set_units(面积, km2),
面积 = as.numeric(面积)) %>%
select(-contains("类型")) %>%
haven::write_dta("2020年各城市行政区划面积.dta")
然后再统计、计算、绘图即可：
*- 各地级市的宗庙祠堂密度
*- 首先根据经纬度判断所处的省市区县
shp2dta using 2020行政区划/县.shp, database(county_db) coordinates(county_coord) genid(ID) replace
import excel using "全国宗庙祠堂分布（10km 范围的检索）.xlsx", clear first
geoinpoly lat lng using county_coord.dta
ren _ID ID
merge m:1 ID using county_db
drop if missing(lng)
drop 省 省代码 市 县 县代码 县类型 市类型 省类型 _merge
gen id = _n
*- 计算每个城市的宗庙祠堂密度
collapse (count) count = id, by(市代码)
*- 匹配面积数据
merge 1:1 市代码 using 2020年各城市行政区划面积.dta
replace count = 0 if missing(count)
keep count 市 市代码 面积
*- 面积的单位是平方公里
gen density = count / (面积 / 10000)
label var density "宗庙祠堂密度（个/万平方公里）"
save "各城市宗庙祠堂密度", replace
*- 绘图
use chinacity2020mini_db.dta, clear
merge 1:1 市 市代码 using "各城市宗庙祠堂密度"
keep if _m == 3
replace density = 0 if missing(density)
grmap density using chinacity2020mini_coord.dta, ///
id(ID) osize(vvthin ...) ocolor(white ...) ///
clmethod(custom) clbreaks(0 3 9 26 48 80 686) ///
fcolor("255 255 229" "247 252 185" "217 240 163" "173 221 142" "120 198 121" "65 171 93" "35 132 67" "0 90 50") ///
graphr(margin(medium)) ///
leg(order(2 "<= 3个/万平方公里" 3 "3 ~ 9个/万平方公里" ///
4 "9 ~ 26个/万平方公里" 5 "26 ~ 48个/万平方公里" ///
6 "48 ~ 80个/万平方公里" 7 "> 80个/万平方公里")) ///
line(data(chinaprov2010mini_line_coord2.dta) by(group) ///
size(vvthin *1 *0.5 *0.5 *0.5) pattern(solid ...) ///
select(drop if inlist(group, 4, 7)) ///
color(gs20 /// 省界颜色
black /// 国界线颜色
"0 85 170" /// 海岸线颜色
black /// 小地图框格颜色
black /// 比例尺和指北针颜色
)) ///
polygon(data(polygon2) fcolor(black) ///
osize(vvthin)) ///
label(data(chinacity2020mini_label2) x(X) y(Y) label(cname) length(20) size(*0.8)) ///
ti("各城市宗庙祠堂密度") ///
subti("数据整理 & 绘制：微信公众号 RStata") ///
caption("数据来源：百度地图地点检索接口", size(*0.8))
gr export "各城市宗庙祠堂密度（10km 范围的检索）.png", replace width(4800)
这样我们就解决了这个问题。
课程信息
该课程已经讲解完，感兴趣的小伙伴可以点击文末的阅读原文跳转到平台上观看学习！
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要报名 RStata 培训班，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/1a86473ee3c14a1c8f4a1ac55d35d05d