twoway 系列命令绘制地图

在谈到用stata绘制地图时，很多人首先想到的是spmap命令。其实，最常见的twoway系列命令就可以完成我们的大多数绘图需求。

对于spmap命令的使用，在线资源已经有很多了，这主要介绍使用最常用的twoway系列命令如何进行地图可视化。spmap命令的使用只做简单介绍。(后台回复20200412获取相关数据)

首先，使用shp2dta命令将 shp 文件转为 stata 能够识别的 .dta 文件。shp2dta语法格式如下：

 shp2dta using shpfilename, database(filename) coordinates(filename) [options]

    options                    Description
    -----------------------------------------------------------------------------------------------------------------------------------------
    Main
      database(filename)       Specifies filename of new dBase dataset
      coordinates(filename)    Specifies filename of new coordinates dataset
      replace                  Overwrites existing database and coordinates filenames
      genid(newvarname)        Creates unique id variable for database.dta
      gencentroids(stub)       Creates centroid variables
    -----------------------------------------------------------------------------------------------------------------------------------------

database(filename)用于将shp文件的属性表导出为 .dta 格式，名称为 filename；
coordinates(filename)则用于导出地理边界坐标点信息，是绘制底图所用的关键数据。
gencentroids(stub)用与生成中心点坐标值

下面是具体演示案例：

导出数据

cd E:/微信/stata地图可视化
* 导入面状文件
shp2dta using 14sheng_matchWGS.shp, database(data_shp) coor(coor0) genid(id) genc(centroid) replace
* 导入线状文件
shp2dta using 三级分界线WGS.shp, database(data_shp3) coor(coor3) genid(id) genc(centroid) replace

绘图

# 调用数据
use coor0, clear
rename _Y lan
rename _X lon
# 将面和线状文件的边界坐标数据进行合并
merge using coor3
drop _merge
merge using data_shp
# 绘制图形
twoway area lan lon, cmiss(n) col(gray) || scatter y_centroid x_centroid, mcolor(dark) msymbol("0") msize(vsmall) || line _Y _X, legend(off)

image

区域上色

如果想要分级显示，则需要分别进行绘制，这是这样方式绘图的不足之处。但是优点在于，对于各个图层的处理有很大的自由度。演示如下：

# 导入属性表
use data_shp, clear
# 与我们的其他经济变量数据进行合并
merge 1:1 xian_code using data2008
rename id _ID
drop _merge
save shuxingbiao.dta, replace
# 代开地理边界数据
use coor0.dta, clear
# 与以及合并的属性表再次合并
merge m:1 _ID using shuxingbiao.dta
# 重命名变量，放置与后面所有合并的表冲突
rename _Y lan
rename _X lon
drop _merge
# 合并分界线的地理边界数据
merge using coor3
# 正式绘图
twoway area lan lon, cmiss(n) col(gray) || area lan lon if crank3 == 1, cmiss(n) col(khaki) || scatter y_centroid x_centroid, mcolor(dark) msymbol(circle_hollow) msize(vtiny) || line _Y _X, legend(off)

image

spmap 绘图

* 使用 shp2dta 转换数据格式
shp2dta using nycwrpczb_201601, database(nydb) coordinates(nycd) genid(id) genc(centroid) replace //行政边界 底图
shp2dta using nycwrpsnwa_201601, database(wadb) coordinates(wacd) genid(id) genc(centroid) replace  //面状图
shp2dta using nycwrprec_201601, database(recdb) coordinates(reccd) genid(id) genc(centroid) replace  //点图
* 地图绘制
use "reccd", clear 
spmap using "nycd.dta", id(_ID) fcolor(purple) ocolor(gray) mosize(medium) osize(medium) /// 
    polygon(data("wacd.dta") fcolor(orange%40) legenda(on) leglabel(special areas)) /// 
    point(data("reccd.dta") x(_X) y(_Y) size(*1.2) fcolor(red) ocolor(black) osize(*0.5) /// 
    legenda(on) leglabel(event pints)) ///
    title("Mapping with Stata", size(*0.8)) /// 
    legend(region(lcolor(gray%0) fcolor(white)) ///  
    title("Legend", size(*0.4) bexpand justification(center)) position (4))

image
不知是否为本人个人问题，发现spmap只允许每种类型的子图叠加一次，希望后面有所改进。

往期回顾

1. 如何使用R、matlab、arcgis进行趋势面分析
2. 使用R进行空间自相关检验
3. 空间权重矩阵的那些事(八)-球面距离权重矩阵
4. R中空间权重矩阵的处理
5. 使用R进行地图可视化
6. 使用QGIS进行地图可视化
7. 南海九段线可以这样处理
8. Stata 绘图也可以很好看

image