Stata:空间计量之用-spmap-绘制地图.md

发布时间:2020-10-09 阅读 14662

Stata 连享会   主页 || 视频 || 推文

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

课程详情 https://gitee.com/arlionn/Course   |   lianxh.cn

课程主页 https://gitee.com/arlionn/Course

万莉 (北京航空航天大学,wanli_buaa@163.com)

左祥太 (武汉纺织大学,wanli_buaa@163.com


目录


引言

在绘制地图时,大家可能会推荐 ArcGIS 等专业软件。然而,对于习惯使用 Stata 的用户而言,若数据管理、数据分析、绘图等各方面都能在 Stata 实现,那就太棒了!Stata 的绘图功能日益强大,能满足大多数可视化需求。现在,Stata 能通过外部命令 spmap 绘制出各式各样的地图,如下图所示:

示例图.png
示例图.png

本文将详细介绍如何使用 spmap 命令绘制各种常用的地图。

1. 命令简介

外部命令 spmap (Pisati, 2007) 的前身为 tmap (Pisati, 2004)。Stata 15 版本的官方命令 grmap 改编自 spmap。相较 spmapgrmap 改动不大,但只用于空间数据 (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) 两个文件。 shp2dtaspshape2dta 命令能帮助我们分别将这两个文件转成 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

  • 上述方法用时较长,这里提供一种替代方案:

    • 首先打开 阿里云地图选择器 ,从页面左下角可以看见数据来源于 高德开放平台 ,在页面右方选择需要下载的地区,然后点击下载,可以下载 JSON 文件或者 SVG 文件
    image-20210913175122084
    image-20210913175122084
    • 下载完成之后导入 mapshaper 这个在线转化格式的网站,选择文件之后( JSON 和 SVG 均可)点击导入
    image-20210913175302288
    image-20210913175302288
    • 点击右上角的 Export ,选择 Shapefile 导出
    image-20210913175510089
    image-20210913175510089

第三步:将地图的矢量数据文件转成 Stata 格式数据文件

利用命令 shp2dtamif2dta 将地图的矢量数据文件转成 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 文件如下: Translate the files.png

第四步:将作图所需数据文件转换成 Stata 格式文件

具体而言,将我们需要作图的数据(attribute variable,比如,某年的人口数据)转换成 Stata 格式的 .dta 数据。

在绘制地图时,我们主要是通过标签数据中的 ID 变量连接(合并)地图数据和作图数据。我们需要保证不同数据中的 ID 变量对应的区域划分是一致的。

本例中,我们作图的数据为美国各州 1990 年的人口数据:stats.dta。 由于 stats.dta 中各州的 ID 变量为 scode ,与 usdb.dta 中的 ID 变量 id 编码不一致,我们生成一个中间文件 trans.dta,将 scodeid 对应起来。

*-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
美国地图.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)
image-20210913180706087
image-20210913180706087

连享会计量方法专题……

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
example1.png
example1.png
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
example2.png
example2.png
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
example3.png
example3.png
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
example4.png
example4.png

关于代码的一些说明:

  • 官方直接提供了 dta 文件,所以我们在绘制地图时省略 shp2dtamif2dta 转换地图矢量数据文件的步骤。
  • 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
example5.png
example5.png
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
example6.png
example6.png

关于代码的一些说明:

  • 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
example7.png
example7.png
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
example8.png
example8.png
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
example9.png
example9.png
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
example10.png
example10.png
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
example11.png
example11.png

关于代码的一些说明:

  • polygon(...) 在底图 (basemap) 上添加其他多边形。
  • point(...) 在底图 (basemap) 上添加点。
  • line(...) 在底图 (basemap) 上添加线条。
  • 绘制比较复杂的地图时,常用上述三个命令在底图 (basemap) 上叠加图层。

小结

用 Stata 绘制地图有两个关键点:一是获得地图的矢量数据文件,熟悉数据结构;二是熟悉命令 spmap 的语法结构。

其它相关命令

  • help tmap 外部命令安装:ssc install tmap, replace 该命令是 spmap 的前身。
  • help grmap 该命令源自 spmap,注意仅适用于 Stata 15 及更高版本。

相关链接

注:本推文相关数据,do file及资料点这里可获得→ 点这里(<ゝω・)☆

相关课程

连享会-直播课 上线了!
http://lianxh.duanshu.com

免费公开课:


课程一览

支持回看,所有课程可以随时购买观看。

专题 嘉宾 直播/回看视频
最新专题 DSGE, 因果推断, 空间计量等
Stata数据清洗 游万海 直播, 2 小时,已上线
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]

Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。


关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
  • 连享会-主页知乎专栏,300+ 推文,实证分析不再抓狂。
  • 公众号推文分类: 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类,主流方法介绍一目了然:DID, RDD, IV, GMM, FE, Probit 等。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD