温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
万莉 (北京航空航天大学,wanli_buaa@163.com)
左祥太 (武汉纺织大学,wanli_buaa@163.com)
目录
在绘制地图时,大家可能会推荐 ArcGIS 等专业软件。然而,对于习惯使用 Stata 的用户而言,若数据管理、数据分析、绘图等各方面都能在 Stata 实现,那就太棒了!Stata 的绘图功能日益强大,能满足大多数可视化需求。现在,Stata 能通过外部命令 spmap
绘制出各式各样的地图,如下图所示:
本文将详细介绍如何使用 spmap
命令绘制各种常用的地图。
外部命令 spmap
(Pisati, 2007) 的前身为 tmap
(Pisati, 2004)。Stata 15 版本的官方命令 grmap
改编自 spmap
。相较 spmap
, grmap
改动不大,但只用于空间数据 (spset data)。
由于篇幅有限,本推文只介绍 spmap
。由于其命令的选项 (options) 较复杂,我们先熟悉绘制地图的一般步骤;再进一步了解其语法结构。
参考 「Stata 官网提供的 FAQ」,绘制地图的一般步骤如下:
值得注意的是,在 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 压缩文件,包含多个文件,其中必须包括以下三个基本文件:
用 Stata 绘制地图时,我们只需要 shape (*.shp) 和 dBase (*.dbf) 两个文件。 shp2dta 与 spshape2dta 命令能帮助我们分别将这两个文件转成 Stata 专属的 .dta 格式的数据。(最近在使用 shp2dta
命令时发现容易报错,报错内容为 “invalid dbase data type” ,所以这里推荐使用 spshape2dta 命令进行文件格式转换操作)
MapInfo Interchange Format 是美国公司 MapInfo 开发的通用数据交换格式。这种格式是 ASCⅡ 码,可以编辑,易于生成,且可以工作在 MapInfo 支持的所有平台上。它将 MapInfo 数据保存于以下两个文件:
用 Stata 绘制地图时, mif2dta
命令能帮助我们分别将 *.mif 和 *.mid 文件转成 Stata 专属的 .dta 格式的数据。
操作步骤:
可以直接通过 全国地理信息资源目录服务系统 (webmap.cn) 下载 shape file
上述方法用时较长,这里提供一种替代方案:
利用命令 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 文件如下:
具体而言,将我们需要作图的数据(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
命令进行画图!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
连享会计量方法专题……
该命令的基本语法如下:
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
设置地图底图。
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
获得更多详细信息。
目前,我们已经熟悉了绘制地图的一般步骤和语法结构。接下来我们将结合具体的例子进一步掌握 spmap
的使用方法。本节 2.2-2.4 的例子均来自 help spmap
帮助文件。官方提供的数据集可通过下述方式获得:
net get spmap.pkg, from(http://fmwww.bc.edu/RePEc/bocode/s)
用 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)
连享会计量方法专题……
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
关于代码的一些说明:
shp2dta
或 mif2dta
转换地图矢量数据文件的步骤。ndfcolor(red)
用红色表示缺失数据区域。clmethod(...)
设置作图数据分组方式,默认是根据分位数分为 4 组。
clmethod(eqint) clnumber(5) eirange(20 70)
:将作图数据按 20-70 范围等分为 5 组。clunmber(#)
设置分组数。legstyle(#)
设置图例显示模式,默认为 1,可选 2,3。
legend(...)
设置图例是否带框线、摆放位置等。fcolor(...)
设置填充颜色。icolor(...)
设置背景色。ocolor(...)
设置区域边界颜色。osize(...)
设置区域边界线宽度。polygon(...)
在底图 (basemap) 上添加其他多边形。scalebar(...)
设置比例尺。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) 上添加标签(比如地名,数字等)。我们还可以通过 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) 上添加线条。
用 Stata 绘制地图有两个关键点:一是获得地图的矢量数据文件,熟悉数据结构;二是熟悉命令 spmap
的语法结构。
help tmap
外部命令安装:ssc install tmap, replace
该命令是 spmap
的前身。help grmap
该命令源自 spmap
,注意仅适用于 Stata 15 及更高版本。
注:本推文相关数据,do file及资料点这里可获得→ 点这里(<ゝω・)☆
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看,所有课程可以随时购买观看。
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
⭕ Stata数据清洗 | 游万海 | 直播, 2 小时,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD