使用 Stata 绘制湖南省双变量填充地图

今天给大家分享一份使用 Stata 绘制湖南省双变量填充地图的数据和方法。

首先需要准备我特别为湖南地图设置的双变量填充地图图例：hn_biscale4x4 文件夹。使用 shp2dta 命令转换成 dta 文件再处理下：

cd "~/Desktop/使用 Stata 绘制湖南省双变量填充地图/"

*- 把 shp 文件转换成 dta 文件
local name = "hn_biscale4x4"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID2) gencentroids(centroids) replace

*- 绘图
use hn_biscale4x4_db.dta, clear
drop ID2
order ID
save hn_biscale4x4_db.dta, replace

再准备图例的标签位置数据：

import excel using "hn_biscale4x4/hn_biscale4x4_label.xlsx", clear first
replace label = "年逆温天数 {&rarr}" in 1
replace label = "年均PM2.5浓度 {&rarr}" in 2
ren x X
ren y Y
ren label cname
save "hn_biscale4x4_label", replace

hnmap_with_ns 文件夹是我准备的带指北针和比例尺的湖南省县级地图：

local name = "hnmap_with_ns"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) gencentroids(centroids) replace
shp2dta using `name'/`name'_line, database(`name'_line_db) coordinates(`name'_line_coord) genid(ID) gencentroids(centroids) replace

简单处理下：

use hnmap_with_ns_db.dta, clear
*- 删除比例尺和指北针
drop if index(县, "_")
save hnmap_with_ns_db.dta, replace

use hnmap_with_ns_coord.dta, clear
*- 保留指北针和比例尺对应的观测值
keep if _ID > 122
gen value = 1
save hn_polyon, replace

由于 Stata 绘制地图只支持添加一个面数据，所以需要把 hn_biscale4x4 合并到 hn_polyon.dta 数据上面：

use hn_polyon, clear
gen class = "指北针和比例尺"
append using hn_biscale4x4_coord
replace class = "双变量图例" if missing(class)

*- 如果使用的双变量图例阶数很高，可能会出现 _ID 重复的情况，所以这里我们最后给指北针和比例尺的 _ID 调大点
replace _ID = _ID + 100 if class == "指北针和比例尺"
save hn_polyon_with_biscale, replace

这样我们就准备好了所有的数据文件：

• hnmap_with_ns_db：底图数据库文件
• hnmap_with_ns_coord：底图坐标系文件
• hnmap_with_ns_line_coord：线条数据
• hn_polyon_with_biscale：多边形数据
• hn_biscale4x4_label：文本标签数据

再准备下绘制地图的数据，这里我们以湖南省各区县的年均 PM2.5 浓度和逆温天数为例：

*- 统计 2021 年每个区县的逆温天数和年均 PM2.5 浓度
import excel using "1998~2021年中国各区县PM2_5面板数据.xlsx", clear first
keep if 省 == "湖南省" & 年份 == 2021
drop 年份 sum min max sd
ren mean PM2_5
save "PM2_5", replace

import excel using "1980～2022年各区县年逆温天数.xlsx", clear first
keep if 省 == "湖南省" & year == 2021
keep 县代码 days1
ren days1 逆温天数
merge 1:1 县代码 using PM2_5
drop _m

*- 对两个变量进行分段，都分成四段，这里我们使用分位数分段
egen group1 = cut(PM2_5), group(4)
replace group1 = group1 + 1
tab group1

egen group2 = cut(逆温天数), group(4)
replace group2 = group2 + 1
tab group2

tostring group1 group2, replace
gen groupclass = group1 + " - " + group2
drop group1 group2

*- 对 groupclass 进行有序因子化
merge m:1 groupclass using hn_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

*- 和 hnmap_with_ns_db 匹配
drop 省 省代码 市 市代码 县
merge 1:1 县代码 using "hnmap_with_ns_db"
replace groupclass2 = -1 if missing(groupclass2)
save "mapdata", replace

然后就可以绘制地图了：

*- 绘制地图
use mapdata, clear

*- 16 种颜色
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""'

grmap groupclass2 using hnmap_with_ns_coord, ///
  id(ID) osize(vvvthin ...) ocolor(white ...) ///
  clmethod(custom) clbreaks(-1 0 1(1)16) ///
    fcolor("gs12" `colorlist') ///
    leg(off) ///
  graphr(margin(medium)) ///
  line(data(hnmap_with_ns_line_coord) by(_ID) ///
    select(keep if _ID > 122) ///
    size(*0.3 ...) pattern(solid ...) ///
    color(black ...)) ///
  polygon(data(hn_polyon_with_biscale) by(_ID) ///
    osize(vvvthin ...) ocolor(`colorlist' black ...) ///
    fcolor(`colorlist' black ...)) ///
  label(data(hn_biscale4x4_label) x(X) y(Y) ///
    select(keep if inlist(cname, "N", "100km") | index(cname, "{&rarr}")) ///
    by(angle) angle(0 90) ///
    label(cname) length(40 ...) size(*0.7 ...)) ///
  ti("2021 年湖南各区县年均 PM2.5 浓度与年总逆温天数", size(*0.8)) ///
  subti("数据处理&绘制：微信公众号 RStata", size(*0.8)) ///
  caption("数据来源：华盛顿大学路易斯分校、NASA", size(*0.6))

gr export pic4x4.png, replace width(4800)
gr export pic4x4.pdf, replace

代码中的 colorlist 需要 16 种颜色，为了便于生成这 16 种颜色，我设计了一个小工具：https://observablehq.com/d/539c972c4d48428b

使用这个工具就可以快速生成双变量填充色了。

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制湖南省双变量填充地图