使用 Stata 绘制双变量填充中国地图
为了让大家更好理解本文内容，欢迎各位培训班小伙伴参加明晚 8 点的直播课：「使用 Stata 绘制双变量填充中国地图」
在之前「使用 Stata 绘制历年中国各省市区县地图（小地图版本 + 长版）」课程的基础上，我们今天再来学习下如何绘制双变量填充地图，之前的课程链接：
使用 Stata 绘制历年中国省级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/e194d75fe1674f7c8921daec1ece7f7d 使用 Stata 绘制历年中国市级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/0e9d78633ae64fefa684d6aad89633bc 使用 Stata 绘制历年中国县级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/99c2e97b88ea401bbdb8da4511341cc8
其中省市的有视频讲解，任意学习一个就会使用这套数据了。
使用 Stata 绘制双变量填充中国地图也同样包含两种版本：
小地图版本
在附件中我提供了 2x2、3x3、4x4、5x5 及 8x8 的双变量图例 shp 数据。课程中我们将以 4x4 的为例进行讲解。
首先我们需要读取和处理瞪羚企业数据，由于要绘制带小地图版本的中国地图，所以我们还需要为小地图部分准备数据：
*- 读取和处理瞪羚数据
import excel using "瞪羚、独角兽、创新型企业经纬度数据.xlsx", clear first
drop if missing(经度)
*- ssc install geo2xy
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
把 biscale_4x4 的 shp 文件转换成 dta 文件：
local name = "biscale4x4"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID2) gencentroids(centroids) replace
use biscale4x4_db.dta, clear
drop ID2
order ID
save biscale4x4_db.dta, replace
图例的标签位置数据也需要读取处理下：
import excel using "biscale4x4/biscale4x4_label.xlsx", clear first
replace label = "瞪羚企业 {&rarr}" in 1
replace label = "独角兽和创新企业 {&rarr}" in 2
ren x X
ren y Y
ren label cname
save "biscale4x4_label", replace
这里的变量重命名是为了等下和地图的标签数据合并。
下面我们就可以把这两个数据和绘图使用的地图数据合并起来了。
这里我准备绘制的是 市级填充地图 + 省级线条 + 散点图，需要下面的这些数据：
chinacity2020mini_coord.dta chinacity2020mini_db.dta chinacity2020mini_label2.dta chinaprov2020mini_line_coord2.dta polygon2.dta
把 biscale4x4 合并到 polygon2.dta 数据上面：
use polygon2, clear
gen class = "指北针和比例尺"
append using biscale4x4_coord
replace class = "双变量图例" if missing(class)
*- 如果使用的双变量图例阶数很高，可能会出现 _ID 重复的情况，所以这里我们最后给指北针和比例尺的 _ID 调大点
replace _ID = _ID + 100 if class == "指北针和比例尺"
save polygon2_with_biscale, replace
还需要把 biscale4x4_label 合并到 chinacity2020mini_label2.dta 文件中：
use chinacity2020mini_label2.dta, clear
append using biscale4x4_label
replace ename = cname if missing(ename)
replace angle = 0 if missing(angle)
save chinacity2020mini_label2_biscale, replace
为了绘制双变量填充地图，我们需要准备两个变量，这里我统计了每个城市的瞪羚企业数量和其他企业的数据量：
*- 统计每个城市瞪羚企业和其他企业的数量
import excel using "瞪羚、独角兽、创新型企业经纬度数据.xlsx", clear first
replace 认定分类 = "瞪羚企业" if index(认定分类, "瞪羚")
replace 认定分类 = "独角兽和创新企业" if 认定分类 != "瞪羚企业"
contract 认定分类 市 市代码
drop if missing(市)
spread 认定分类 _freq
replace 独角兽和创新企业 = 0 if missing(独角兽和创新企业)
replace 瞪羚企业 = 0 if missing(瞪羚企业)
然后需要对两个变量进行分组：
*- 对两个变量进行分段，都分成四段，这里我们使用分位数分段
egen group1 = cut(瞪羚企业), group(4)
replace group1 = group1 + 1
tab group1
egen group2 = cut(独角兽和创新企业), group(4)
replace group2 = group2 + 1
tab group2
*- 可以看到如果数据分布太不均匀，这种分组方法可能会失效，这个时候可以手动指定断点：
cap drop group2
egen group2 = cut(独角兽和创新企业), at(0 1 2 10 341) icodes
replace group2 = group2 + 1
tab group2
tostring group1 group2, replace
gen groupclass = group1 + " - " + group2
drop group1 group2
不过这个 groupclass 变量是个字符串，不能直接用来绘图，我们需要对其进行有序因子化，也就是使 1-1、1-2、1-3、1-4、2-1、2-1、...、4-4 恰好对应 1、2、3、4、...、16。
*- 对 groupclass 进行有序因子化
merge m:1 groupclass using biscale4x4_db
drop if _m == 2
drop _m
*- 安装 labmask: net install gr0034.pkg, from("http://www.stata-journal.com/software/sj8-2") replace
labmask ID, val(groupclass)
ren ID groupclass2
drop x_centroids - fill
*- 和 chinacity2020mini_db.dta 匹配
merge 1:1 市 市代码 using "chinacity2020mini_db.dta"
replace groupclass2 = -1 if missing(groupclass2)
save "mapdata", replace
然后就可以绘图了：
*- 绘制地图
use mapdata, clear
gsort 市代码
local colorlist = `""232 232 232" "189 222 222" "142 212 212" "90 200 200" "218 189 212" "189 189 212" "142 189 212" "90 189 200" "204 146 193" "189 146 193" "142 146 193" "90 146 193" "190 100 172" "189 100 172" "142 100 172" "90 100 172""'
local colorlist2 = `""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""'
grmap groupclass2 using chinacity2020mini_coord.dta, ///
id(ID) osize(vvvthin ...) ocolor(white ...) ///
clmethod(custom) clbreaks(-1 0 1(1)16) ///
fcolor("gs12" `colorlist') ///
leg(off) ///
graphr(margin(medium)) ///
line(data(chinaprov2020mini_line_coord2.dta) by(group) ///
size(vvthin *1 *0.5 *0.5 *0.5) pattern(solid ...) ///
select(drop if inlist(group, 4, 7)) ///
color(gs5 /// 省界颜色
"162 154 196" /// 国界线颜色
"0 85 170" /// 海岸线颜色
black /// 小地图框格颜色
black /// 比例尺和指北针颜色
)) ///
polygon(data(polygon2_with_biscale) by(_ID) ///
osize(vvvthin ...) ocolor(white ...) ///
fcolor(`colorlist' black ...)) ///
label(data(chinacity2020mini_label2_biscale) x(X) y(Y) ///
select(keep if inlist(cname, "N", "1000km") | index(cname, "{&rarr}")) ///
by(angle) angle(0 90) ///
label(cname) length(40 ...) size(*0.7 ...)) ///
point(data(pointdata) by(省代码) x(lon) y(lat) ///
fcolor(`colorlist2') shape(o ...) size(vtiny ...)) ///
ti("瞪羚、独角兽、创新型企业城市分布", size(*0.8)) ///
subti("数据爬取&绘制：微信公众号 RStata", size(*0.8)) ///
caption("数据来源：瞪羚云网站", size(*0.6))
gr export pic4x4.png, replace width(4800)
代码中的 colorlist 是需要 16 种颜色，为了便于生成这 16 种颜色，我设计了一个小工具：https://observablehq.com/d/539c972c4d48428b
使用这个工具就可以快速生成双变量填充色了。
colorlist2 是各个省份的散点颜色，可以使用我设计的这个工具：https://tidyfriday.cn/colors
相比于原始的 colorbrewer 网站，这个工具会对 Stata 和 R 语言用户更友好。
这样就可以直接绘制出来这种效果了，无需后期再修图。
长版地图
长版地图的绘制方法类似上面。具体可以参考附件中的 main_long.do 文件：
其他阶数
上面的内容是以 4x4 为例的，附件中也提供了 2x2、3x3、4x4、5x5 及 8x8 的数据，也可以绘制这些阶数的。
直播信息
欢迎各位培训班会员参加明晚 8 点的直播课：「使用 Stata 绘制双变量填充中国地图」！
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/36da18e9b108474b910559de78cca942