使用 Stata 绘制广东省市级双变量填充地图

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

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

cd "~/Desktop/使用 Stata 绘制广东省市级双变量填充地图"

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

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

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

import excel using "gd_biscale4x4/gd_biscale4x4_label.xlsx", clear first
replace label = "年逆温天数 {&rarr}" if label == "v1 →"
replace label = "年均PM2.5浓度 {&rarr}" if label == "v2 →"
ren x X
ren y Y
ren label cname
save "gd_biscale4x4_label", replace

gdmap_with_ns 文件夹是我准备的带指北针和比例尺的广东省市级地图：

local name = "gdmap_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 gdmap_with_ns_db.dta, clear
*- 删除比例尺和指北针
drop if index(市, "_")
save gdmap_with_ns_db.dta, replace

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

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

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

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

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

• gdmap_with_ns_db：底图数据库文件
• gdmap_with_ns_coord：底图坐标系文件
• gdmap_with_ns_line_coord：线条数据
• gd_polyon_with_biscale：多边形数据
• gd_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 gd_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

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

*- 线条数据
use gdmap_with_ns_line_db, clear
use gdmap_with_ns_line_coord, clear
keep if _ID > 21

然后就可以绘制地图了：

*- 绘制地图
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 gdmap_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(gdmap_with_ns_line_coord) by(_ID) ///
    select(keep if _ID > 21) ///
    size(*0.3 ...) pattern(solid ...) ///
    color(black ...)) ///
  polygon(data(gd_polyon_with_biscale) by(_ID) ///
    osize(vvvthin ...) ocolor(`colorlist' black ...) ///
    fcolor(`colorlist' black ...)) ///
  label(data(gd_biscale4x4_label) x(X) y(Y) ///
    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 绘制广东省市级双变量填充地图