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

教育   2024-10-25 16:02   安徽  


本文内容较为简单,就不再进行视频讲解了~


今天给大家分享一份使用 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 培训班获取全部课程和以会员价获取数据资料(10元/份)详情可阅读这篇推文:数据处理、图表绘制、效率分析与计量经济学如何学习~

详情可点击阅读原文进入 RStata 学院了解(从首页的会员卡专区即可查看和购买会员卡)。

更多关于 RStata 培训班的信息可添加微信号 r_stata 咨询:

附件下载(点击文末的阅读原文即可跳转):

https://rstata.duanshu.com/#/brief/course/d77667cb17b34a0bb76d60fdaa746898


RStata
一起学习 R 语言和 Stata 吧!
 最新文章