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

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

该数据的使用方法类似之前讲解的 使用 Stata 绘制双变量填充中国地图，感兴趣的小伙伴可以观看视频讲解：https://rstata.duanshu.com/#/brief/course/36da18e9b108474b910559de78cca942

首先需要准备我特别为湖南地图设置的双变量填充地图图例：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

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