本文内容较为简单，就不再进行视频讲解了～
今天给大家分享一份使用 Stata 绘制江苏省区县地图+湖泊河流分布的数据和方法。
江苏省区县地图
首先需要准备我特别设计的江苏地图数据（包含指北针和比例尺）：
jsmap 文件夹：地理矢量数据； labeldf.dta 数据：文本标签位置数据；
这份数据是使用 R 语言设计的，感兴趣的小伙伴可以参考附件中的 main.R 文件。
下面我们看一下如何在 Stata 中绘制江苏区县地图。
首先把 shp 文件转换成 dta 文件：
*- 把 shp 文件转换成 dta 文件
local name = "jsmap"
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
这样就会得到四个文件：
jsmap_db.dta jsmap_coord.dta jsmap_line_coord.dta jsmap_line_db.dta
其中前两个是底图的属性数据和坐标系数据，后面两个是线条的。
下面还需要对数据进行下调整。
首先是删除底图上的比例尺和指北针数据，因为这两个元素计划使用线条和多边形选项绘制：
use jsmap_db.dta, clear
*- 删除比例尺和指北针
drop if index(县, "_")
save jsmap_db.dta, replace
再创建比例尺和指北针的多边形数据：
use jsmap_coord.dta, clear
*- 保留指北针和比例尺对应的观测值
keep if _ID > 96
gen value = 1
save jsmap_polyon, replace
再对线条进行分组：
*- 线条
use jsmap_line_db.dta, clear
*- 分组
gen group = 1
replace group = 2 if index(市, "边界")
replace group = 3 if index(省, "边界")
replace group = 4 if index(省, "_")
keep ID group
ren ID _ID
save temp, replace
use jsmap_line_coord.dta, clear
merge m:1 _ID using temp
drop _m
save jsmap_line_coord, replace
这里我把线条分为四类：
区县边界； 城市边界； 省份边界； 比例尺和指北针的线条。
后面的代码中也会根据不同的线条类型使用不同的颜色和粗细。
统计 2021 年每个区县的年均 PM2.5 浓度：
import excel using "1998~2021年中国各区县PM2_5面板数据.xlsx", clear first
keep if inlist(省, "江苏省") & 年份 == 2021
drop 年份 sum min max sd
ren mean PM2_5
keep 县代码 PM2_5
save "PM2_5", replace
然后就可以绘制地图了：
use jsmap_db.dta, clear
merge 1:1 县代码 using PM2_5
drop if _m == 2
drop _m
*- 缺失值替换成 -1
replace PM2_5 = -1 if mi(PM2_5)
hist PM2_5
grmap PM2_5 using jsmap_coord.dta, ///
id(ID) osize(vvvthin ...) ocolor(white ...) ///
clmethod(custom) clbreaks(-1 0 10 20 30 40 50 60) ///
fcolor("gs12" "237 231 246" "209 196 233" "179 157 219" "149 117 205" "126 87 194" "103 58 183") ///
leg(off) ///
graphr(margin(medium)) ///
line(data(jsmap_line_coord.dta) by(group) ///
size(*0.2 *0.7 *1.2 *0.5) pattern(solid ...) ///
color("247 247 247" "150 150 150" black black)) ///
polygon(data(jsmap_polyon.dta) by(value) ///
osize(vvvthin ...) ocolor(black ...) ///
fcolor(black ...)) ///
label(data(labeldf.dta) x(X) y(Y) ///
select(keep if inlist(cname, "N", "150km") | class == "市级标签") ///
label(cname) length(40 ...) size(*1 ...)) ///
ti("2021 年江苏各区县年均 PM2.5 浓度", size(*1)) ///
subti("数据处理&绘制：微信公众号 RStata", size(*1)) ///
caption("数据来源：华盛顿大学路易斯分校", size(*0.7))
添加湖泊河流分布
下面我们看一下如何在江苏区县地图的基础上添加湖泊河流分布。
首先需要准备我特别设计的江苏省湖泊与河流数据：
江苏湖泊文件夹：地理矢量数据； 江苏河流文件夹：地理矢量数据；
这些数据是使用 R 语言设计的，感兴趣的小伙伴可以参考附件中的 提取江苏省湖泊与河流.R 文件。
接下来处理湖泊数据：
local name = "江苏湖泊"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) replace
*- 把这个数据转换坐标系后再和 jsmap_polyon.dta 合并起来即可
use jsmap_polyon, clear
use 江苏湖泊_db, clear
use 江苏湖泊_coord, clear
gen value = 2
geo2xy _Y _X, projection(albers, 6378137 298.257223563 25 47 0 105) replace
append using jsmap_polyon
save jsmap_polyon2, replace
处理河流数据：
local name = "江苏河流"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) replace
*- 把这个数据转换坐标系后再和 jsmap_line_coord.dta 合并起来即可
use jsmap_line_coord, clear
use 江苏河流_db, clear
use 江苏河流_coord, clear
gen group = 5
geo2xy _Y _X, projection(albers, 6378137 298.257223563 25 47 0 105) replace
append using jsmap_line_coord
save jsmap_line_coord2, replace
最后就可以江苏省区县地图上添加湖泊河流分布了：
use jsmap_db.dta, clear
merge 1:1 县代码 using PM2_5
drop if _m == 2
drop _m
*- 缺失值替换成 -1
replace PM2_5 = -1 if mi(PM2_5)
hist PM2_5
grmap PM2_5 using jsmap_coord.dta, ///
id(ID) osize(vvvthin ...) ocolor(white ...) ///
clmethod(custom) clbreaks(-1 0 10 20 30 40 50 60) ///
fcolor("gs12" "237 231 246" "209 196 233" "179 157 219" "149 117 205" "126 87 194" "103 58 183") ///
leg(off) ///
graphr(margin(medium)) ///
line(data(jsmap_line_coord2.dta) by(group) ///
size(*0.2 *0.7 *1.2 *0.5 *0.5) pattern(solid ...) ///
color("247 247 247" "150 150 150" black black "0 85 170")) ///
polygon(data(jsmap_polyon2.dta) by(value) ///
osize(vvvthin ...) ocolor(black "0 85 170") ///
fcolor(black "0 85 170")) ///
label(data(labeldf.dta) x(X) y(Y) ///
select(keep if inlist(cname, "N", "150km") | class == "市级标签") ///
label(cname) length(40 ...) size(*1 ...)) ///
ti("2021 年江苏各区县年均 PM2.5 浓度", size(*1)) ///
subti("数据处理&绘制：微信公众号 RStata", size(*1)) ///
caption("数据来源：华盛顿大学路易斯分校", size(*0.7))
获取附件
是不是感觉很硬核！欢迎报名 RStata 培训班获取全部课程和以会员价获取数据资料（10元/份）详情可阅读这篇推文：数据处理、图表绘制、效率分析与计量经济学如何学习～
详情可点击阅读原文进入 RStata 学院了解（从首页的会员卡专区即可查看和购买会员卡）。
更多关于 RStata 培训班的信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/bf15d66088aa4fd0b7909a4332582f04
