Stata：空间计量之用-spmap-绘制地图.md
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
万莉 (北京航空航天大学，wanli_buaa@163.com)
左祥太 (武汉纺织大学，wanli_buaa@163.com）
目录
引言
在绘制地图时，大家可能会推荐 ArcGIS 等专业软件。然而，对于习惯使用 Stata 的用户而言，若数据管理、数据分析、绘图等各方面都能在 Stata 实现，那就太棒了！Stata 的绘图功能日益强大，能满足大多数可视化需求。现在，Stata 能通过外部命令 spmap 绘制出各式各样的地图，如下图所示：
本文将详细介绍如何使用 spmap 命令绘制各种常用的地图。
1. 命令简介
外部命令 spmap (Pisati, 2007) 的前身为 tmap (Pisati, 2004)。Stata 15 版本的官方命令 grmap 改编自 spmap。相较 spmap， grmap 改动不大，但只用于空间数据 (spset data)。
由于篇幅有限，本推文只介绍 spmap。由于其命令的选项 (options) 较复杂，我们先熟悉绘制地图的一般步骤；再进一步了解其语法结构。
1.1 绘制步骤
参考 「Stata 官网提供的 FAQ」，绘制地图的一般步骤如下：
第一步：下载三个外部命令：spmap | shp2dta | mif2dta
值得注意的是，在 Stata 15 之后，Stata 官方引入了 spshape2dta 命令，可以直接实现对 shape file的 转换
*-Step 1: Obtain and install the spmap, shp2dta, and mif2dta commands
// 画地图的外部命令
. ssc install spmap, replace
// 将外部的 dbf 文件和 shp 文件转化为 dta 文件
. ssc install shp2dta, replace
// 将外部的 mif 文件和 mid 文件转化为 dta 文件
. ssc install mif2dta, replace
* 获得帮助文件
help spmap
help shp2dta
help mif2dta
第二步：下载描绘地图的矢量数据文件
绘制地图时，一般需要两份数据：标签数据 (比如 ID 变量) 和地理坐标数据。对于后者，则有两种常用的矢量数据存储格式： ESRI shapefile （简称 shapefile）或 MapInfo Interchange Format。具体说明如下：
ESRI shapefile 是美国环境系统研究所公司 (ESRI) 开发的一种空间数据开放格式，用于存储地理要素的位置、形状和属性。Shapefile 通常有一个 .zip 压缩文件，包含多个文件，其中必须包括以下三个基本文件：
主文件 (*.shp)： 存储几何要素的的空间信息，也就是 XY 坐标。 索引文件 (*.shx)： 存储有关 *.shp 存储的索引信息。它记录了在 *.shp 中，空间数据是如何存储的，XY 坐标的输入点在哪里，有多少 XY 坐标对等信息。 表文件 (*.dbf)： 存储地理数据的属性信息的 dBase 表。
用 Stata 绘制地图时，我们只需要 shape (*.shp) 和 dBase (*.dbf) 两个文件。 shp2dta 与 spshape2dta 命令能帮助我们分别将这两个文件转成 Stata 专属的 .dta 格式的数据。（最近在使用
shp2dta命令时发现容易报错，报错内容为 “invalid dbase data type” ，所以这里推荐使用 spshape2dta 命令进行文件格式转换操作）MapInfo Interchange Format 是美国公司 MapInfo 开发的通用数据交换格式。这种格式是 ASCⅡ 码，可以编辑，易于生成，且可以工作在 MapInfo 支持的所有平台上。它将 MapInfo 数据保存于以下两个文件：
*.mif 文件： 存储图形数据，可理解为 XY 坐标。 *.mid 文件： 存储文本（属性）数据。
用 Stata 绘制地图时，
mif2dta命令能帮助我们分别将 *.mif 和 *.mid 文件转成 Stata 专属的 .dta 格式的数据。
操作步骤：
可以直接通过 全国地理信息资源目录服务系统 (webmap.cn) 下载 shape file
上述方法用时较长，这里提供一种替代方案：
image-20210913175122084 下载完成之后导入 mapshaper 这个在线转化格式的网站，选择文件之后（ JSON 和 SVG 均可）点击导入
image-20210913175302288 点击右上角的 Export ，选择 Shapefile 导出
image-20210913175510089
第三步：将地图的矢量数据文件转成 Stata 格式数据文件
利用命令 shp2dta 或 mif2dta 将地图的矢量数据文件转成 Stata .dta 格式的数据。
检查 .dbf (.mid) 文件经转换后对应的 .dta 数据。该数据为标签数据，包含 ID 变量（命名为 id），不同 id 对应不同区域划分。我们应该了解地图矢量数据的提供者是如何划分区域的。比如，在一份数据中， id=1 对应 Alaska，id=2 对应 Alabama；而在另一份数据中，id=1 对应 Albania，id=2 对应 Argentina。
*-Step 3: Translate the files
* note:将数据放在当前工作路径
shp2dta using s_11au16, database(usdb) ///
coordinates(uscoord) genid(id)
/* 如果上述命令存在报错，则使用以下命令 */
spshape2dta s_11au16
关于代码的一些说明：
database(usdb)将 database file 转换成 usdb.dta；coordinates(uscoord)将 coordinate file 转换成 uscoord.dta；genid(id)将 usdb.dta 中的 ID 变量命名为 id;转换 MapInfo files 时，语法一致，只需将 shp2dta 换成 mif2dta。 转换后的两个 .dta 文件如下：
第四步：将作图所需数据文件转换成 Stata 格式文件
具体而言，将我们需要作图的数据（attribute variable，比如，某年的人口数据）转换成 Stata 格式的 .dta 数据。
在绘制地图时，我们主要是通过标签数据中的 ID 变量连接（合并）地图数据和作图数据。我们需要保证不同数据中的 ID 变量对应的区域划分是一致的。
本例中，我们作图的数据为美国各州 1990 年的人口数据：stats.dta。 由于 stats.dta 中各州的 ID 变量为 scode ，与 usdb.dta 中的 ID 变量 id 编码不一致，我们生成一个中间文件 trans.dta，将 scode 与 id 对应起来。
*-Step 4: Determine the coding used by the map
* note:将数据放在当前工作路径
. use stats, clear
. merge 1:1 scode using trans
Result # of obs.
-----------------------------------------
not matched 0
matched 51 (_merge==3)
-----------------------------------------
. drop _merge
*仔细检查合并的结果，确认无误后可删掉 _merge
*- _merge 变量的含义：
* _merge==1 obs. from master data
* _merge==2 obs. from only one using dataset
* _merge==3 obs. from at least two datasets, master or using
第五步：数据合并
将 .dbf (.mid) 文件经转换后对应的 .dta 数据与我们需要作图的 .dta 数据进行合并。
*-Step 5: Merge datasets
* note:将数据放在当前工作路径
. merge 1:1 id using usdb
Result # of obs.
-----------------------------------------
not matched 6
from master 0 (_merge==1)
from using 6 (_merge==2)
matched 51 (_merge==3)
-----------------------------------------
. drop if _merge!=3
(6 observations deleted)
. drop _merge
*仔细检查合并的结果，确定无误后可删掉多余观测值
第六步：用 spmap 命令绘图
最后一步啦！可以用 spmap 命令进行画图！show your data!
*-Step 6: Draw the graph
* 图1：不完美示范
spmap pop1990 using uscoord, ///
id(id) fcolor(Blues)
* 图2：进一步美化
replace pop1990 = pop1990/1e+6
* 改变单位，避免图例里显示的数字过大
format pop1990 %4.2f
* 变量数据格式（整数与小数部分长度）
spmap pop1990 using uscoord ///
if id !=1 & id !=54 & ///
id !=39 & id !=56, ///
id(id) fcolor(Blues)
关于代码的一些说明：
pop1990 是我们作图的数据（1990 年美国各州人口数据）。 using uscoord 是必填项。绘制地图地图时，必需坐标数据（本例子中坐标数据为 uscoord.dta）。 id(id) 是必填项。 绘制地图地图时，必需不同区域划分的 ID 变量（本例子中命名为 id）。 fcolor(Blues) 是可选项，设置专题地图 (Choropleth map) 的填充颜色。本例中选用 Blues 主题，即以蓝色为主题色，不同区域会根据作图数据（本例中为人口数据）大小填充深浅不同的颜色。 clnumber(#) 是可选项，设置分组数。 spmap默认基于分位数，将作图数据（本例中为人口数据）分为四组；若要自定义组数，则可加上该选项。if 命令的使用：本例中使用 if去掉四个州，是因为这四个州分布过于分散，画出的图不美观。两幅图的对比如下： 美国地图.png
连享会计量方法专题……
1.2 命令的语法
该命令的基本语法如下：
spmap [attribute] [if] [in] using basemap [,
basemap_options
polygon(polygon_suboptions)
line(line_suboptions)
point(point_suboptions)
diagram(diagram_suboptions)
arrow(arrow_suboptions)
label(label_suboptions)
scalebar(scalebar_suboptions)
graph_options]
主要可选项说明如下：
basemap_options设置地图底图。其中最重要的选项为 id(idvar)。该选项必须加上，用来识别不同区域划分的 ID 变量 idvar。 其他选项用来控制底图的四个方面： Cartogram（比如背景色）, Choropleth map（比如将作图数据分为几组）, Format（比如填充颜色）, Legend（比如是否显示图例）。
polygon(polygon_suboptions)在底图基础上添加其他多边形。line(line_suboptions)在底图基础上添加线条。point(point_suboptions)在底图基础上添加点。diagram(diagram_suboptions)在底图基础上添加条形/扇形统计图。arrow(arrow_suboptions)在底图基础上添加箭头。label(label_suboptions)在底图基础上添加文字标签。scalebar(scalebar_suboptions)设置比例尺。graph_options与我们常规画图的选项 (twoway_options) 基本一致。help spmap获得更多详细信息。
2. Stata 范例
目前，我们已经熟悉了绘制地图的一般步骤和语法结构。接下来我们将结合具体的例子进一步掌握 spmap 的使用方法。本节 2.2-2.4 的例子均来自 help spmap 帮助文件。官方提供的数据集可通过下述方式获得：
net get spmap.pkg, from(http://fmwww.bc.edu/RePEc/bocode/s)
2.1 如何获得地图的矢量数据文件？
用 Stata 绘制地图至少需要两份数据：带有标签 (比如 ID 变量) 的数据和地理坐标数据。在 spmap 中，标签数据被称为 "master" dataset，地理坐标数据被称为 "basemap" dataset。
因此，我们绘图前必须先下载地图的矢量数据文件。若能直接下载 .dta 格式文件，则无需转换格式；若只能下载 shapefile （或 MapInfo Interchange Format） 格式的文件，则需用 shp2dta （或 mif2dta ）将其转成 .dta 格式。
获得途径总结如下：
在搜索引擎上搜索关键词，比如 "中国地图 shapefile"。 在 GADM 网站上下载 shapefile 文件。该网站提供世界及世界各国的地图数据。 利用 Stata 命令，如下：
*-获得中国地图数据
*-注意：为了保证图形直接可以用于论文（保证领土完整性），这里建议通过 1.1 当中第二步推荐的方式在高德地图开放平台进行下载
*-获得世界地图数据
use "http://fmwww.bc.edu/repec/bocode/w/world-c.dta",clear
save world-c.dta
use "http://fmwww.bc.edu/repec/bocode/w/world-d.dta",clear
save world-d.dta
在获得 JSON / SVG 格式的数据后，先经由 mapshaper 转换为 Shapefile，再使用 spshape2dta 转换为 dta 文件后，在现有的工作路径下会得到两份 dta 文件，其中带有 _shp 后缀的是我们需要作图基础文件，可以理解为地图的轮廓，没有后缀的文件则需要与我们的变量文件（人口、GDP等）进行合并。我们可利用 scatter 检查数据，看地图轮廓是否满足要求。
use "china_shp.dta", clear
scatter _Y _X, msize(vtiny)
连享会计量方法专题……
2.2 Stata 范例 1: Choropleth maps
Choropleth map 可翻译为 分级统计地图。在这种地图中，每个区域以不同深浅度的颜色表示数据变量。此法因常用色级表示，所以也叫色级统计图法。我们可简单理解为 "地图+热力：颜色渐变常见"。
use "Italy-RegionsData.dta", clear
#d ;
spmap relig1 using "Italy-RegionsCoordinates.dta", id(id)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8));
#d cr
use "Italy-RegionsData.dta", clear
#d ;
spmap relig1m using "Italy-RegionsCoordinates.dta", id(id)
ndfcolor(red)
clmethod(eqint) clnumber(5) eirange(20 70)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
legstyle(2) legend(region(lcolor(black)));
#d cr
use "Italy-RegionsData.dta", clear
#d ;
spmap relig1 using "Italy-RegionsCoordinates.dta", id(id)
clnumber(20) fcolor(Reds2) ocolor(none ..)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
legstyle(3) legend(ring(1) position(3));
#d cr
use "Italy-RegionsData.dta", clear
#d ;
spmap relig1 using "Italy-RegionsCoordinates.dta", id(id)
clnumber(20) fcolor(Greens2) ocolor(white ..) osize(thin ..)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
legstyle(3) legend(ring(1) position(3))
plotregion(icolor(stone)) graphregion(icolor(stone))
polygon(data("Italy-Highlights.dta") ocolor(white)
osize(medthick))
scalebar(units(500) scale(1/1000) xpos(-100) label(Kilometers));
#d cr
关于代码的一些说明：
官方直接提供了 dta 文件，所以我们在绘制地图时省略 shp2dta或mif2dta转换地图矢量数据文件的步骤。ndfcolor(red)用红色表示缺失数据区域。clmethod(...)设置作图数据分组方式，默认是根据分位数分为 4 组。括号里可选填 quantile, boxplot, eqint, stdev, kmeans, custom, unique。 clmethod(eqint) clnumber(5) eirange(20 70)：将作图数据按 20-70 范围等分为 5 组。
clunmber(#)设置分组数。legstyle(#)设置图例显示模式，默认为 1，可选 2,3。1 显示为 (num,num]；2 显示为 num-num；3 为色卡。
legend(...)设置图例是否带框线、摆放位置等。fcolor(...)设置填充颜色。icolor(...)设置背景色。ocolor(...)设置区域边界颜色。osize(...)设置区域边界线宽度。polygon(...)在底图 (basemap) 上添加其他多边形。scalebar(...)设置比例尺。
2.3 Stata 范例 2: Proportional symbol maps
Proportional symbol map 可理解为 带气泡的地图 (Bubble Map)。在这种地图中，气泡的面积代表了这个数据的大小。我们可简单理解为 "地图+气泡：圆越大数越大"。
我们可以通过添加选项 point(... proportional(propvar_pn) ...) 完成带气泡的地图。若我们想将气泡（圆）变成其他符号（比如为正方形），可在里面添加子选项 shape(s)。
若在该选项中去掉子选项 proportional(propvar_pn)，则该图可退化为 点描法地图 (dot map)，即气泡变为大小相等的点。我们可简单理解为 "地图+点：哪里有数点哪里"。
use "Italy-OutlineData.dta", clear
#d;
spmap using "Italy-OutlineCoordinates.dta", id(id)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
point(data("Italy-RegionsData.dta") xcoord(xcoord)
ycoord(ycoord) deviation(relig1) fcolor(red) dmax(30)
legenda(on) leglabel(Deviation from the mean));
#d cr
use "Italy-OutlineData.dta", clear
#d;
spmap using "Italy-OutlineCoordinates.dta", id(id)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
point(data("Italy-RegionsData.dta") xcoord(xcoord)
ycoord(ycoord) proportional(relig1) fcolor(green)
ocolor(white) size(*3))
label(data("Italy-RegionsData.dta") xcoord(xcoord)
ycoord(ycoord) label(relig1) color(white) size(*0.7));
#d cr
关于代码的一些说明：
point(...)在底图 (basemap) 上添加点。label(...)在底图 (basemap) 上添加标签（比如地名，数字等）。
2.4 Stata 范例 3: Other maps
我们还可以通过 spmap 画其他类型的图，比如在地图上添加条形图/扇形图。
use "Italy-RegionsData.dta", clear
#d ;
spmap using "Italy-RegionsCoordinates.dta", id(id) fcolor(stone)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Italy, 1994-98" " ", size(*0.8))
diagram(variable(relig1) range(0 100) refweight(pop98)
xcoord(xcoord) ycoord(ycoord) fcolor(red));
#d cr
use "Italy-RegionsData.dta", clear
#d ;
spmap using "Italy-RegionsCoordinates.dta", id(id) fcolor(stone)
diagram(variable(relig1 relig2 relig3) proportional(fortell)
xcoord(xcoord) ycoord(ycoord) legenda(on))
legend(title("Religious orientation", size(*0.5) bexpand
justification(left)))
note(" "
"NOTE: Chart size proportional to number of fortune tellers per million population",
size(*0.75));
#d cr
use "Italy-OutlineData.dta", clear
#d ;
spmap using "Italy-OutlineCoordinates.dta", id(id) fc(bluishgray)
ocolor(none)
title("Provincial capitals" " ", size(*0.9) color(white))
point(data("Italy-Capitals.dta") xcoord(xcoord)
ycoord(ycoord) by(size) fcolor(orange red maroon) shape(s ..)
legenda(on))
legend(title("Population 1998", size(*0.5) bexpand
justification(left)) region(lcolor(black) fcolor(white))
position(2))
plotregion(margin(medium) icolor(dknavy) color(dknavy))
graphregion(icolor(dknavy) color(dknavy));
#d cr
use "Italy-RegionsData.dta", clear
#d ;
spmap relig1 using "Italy-RegionsCoordinates.dta" if zone==1,
id(id) fcolor(Blues2) ocolor(white ..) osize(medthin ..)
title("Pct. Catholics without reservations", size(*0.8))
subtitle("Northern Italy, 1994-98" " ", size(*0.8))
polygon(data("Italy-OutlineCoordinates.dta") fcolor(gs12)
ocolor(white) osize(medthin)) polyfirst;
#d cr
use "Italy-OutlineData.dta", clear
#d ;
spmap using "Italy-OutlineCoordinates.dta", id(id) fc(sand)
title("Main lakes and rivers" " ", size(*0.9))
polygon(data("Italy-Lakes.dta") fcolor(blue) ocolor(blue))
line(data("Italy-Rivers.dta") color(blue) )
freestyle aspect(1.4) xlab(400000 900000 1400000, grid);
#d cr
关于代码的一些说明：
polygon(...)在底图 (basemap) 上添加其他多边形。point(...)在底图 (basemap) 上添加点。line(...)在底图 (basemap) 上添加线条。绘制比较复杂的地图时，常用上述三个命令在底图 (basemap) 上叠加图层。
小结
用 Stata 绘制地图有两个关键点：一是获得地图的矢量数据文件，熟悉数据结构；二是熟悉命令 spmap 的语法结构。
其它相关命令
help tmap外部命令安装：ssc install tmap, replace 该命令是spmap的前身。help grmap该命令源自spmap，注意仅适用于 Stata 15 及更高版本。
相关链接
How do I graph data onto a map with spmap? How do I graph data onto a map with tmap? Pisati, M. 2004. Simple thematic mapping. Stata Journal 4: 361-378. Pisati, M. 2007. spmap: Stata module to visualize spatial data. Statistical Software Components S456812, Department of Economics, Boston College. ShapeFile 格式说明 【地理信息】GIS 之 shapefile 文件趣解（上） 一文看尽 Stata 绘图 [Stata] 地图绘制（一）：中国地图与世界地图 [Stata] 地图绘制（二）: spkde/spgrid/spmap 【盘点】这些经典地图可视化你一定要知道！
注：本推文相关数据，do file及资料点这里可获得→ 点这里(<ゝω・)☆
相关课程
连享会-直播课 上线了！
http://lianxh.duanshu.com
免费公开课：
直击面板数据模型 - 连玉君，时长：1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT，dofiles等)
课程一览
支持回看，所有课程可以随时购买观看。
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
| ⭕ Stata数据清洗 | 游万海 | 直播, 2 小时，已上线 |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。 连享会-主页 和 知乎专栏，300+ 推文，实证分析不再抓狂。 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD
热门资讯
