在談到用stata繪制地圖時(shí)，很多人首先想到的是spmap命令。其實(shí)，最常見的twoway系列命令就可以完成我們的大多數(shù)繪圖需求。

對(duì)于spmap命令的使用，在線資源已經(jīng)有很多了，這主要介紹使用最常用的twoway系列命令如何進(jìn)行地圖可視化。spmap命令的使用只做簡單介紹。(后臺(tái)回復(fù)20200412獲取相關(guān)數(shù)據(jù))

首先，使用shp2dta命令將 shp 文件轉(zhuǎn)為 stata 能夠識(shí)別的 .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文件的屬性表導(dǎo)出為 .dta 格式，名稱為 filename；
coordinates(filename)則用于導(dǎo)出地理邊界坐標(biāo)點(diǎn)信息，是繪制底圖所用的關(guān)鍵數(shù)據(jù)。
gencentroids(stub)用與生成中心點(diǎn)坐標(biāo)值

下面是具體演示案例：

導(dǎo)出數(shù)據(jù)

cd E:/微信/stata地圖可視化
* 導(dǎo)入面狀文件
shp2dta using 14sheng_matchWGS.shp, database(data_shp) coor(coor0) genid(id) genc(centroid) replace
* 導(dǎo)入線狀文件
shp2dta using 三級(jí)分界線WGS.shp, database(data_shp3) coor(coor3) genid(id) genc(centroid) replace

繪圖

# 調(diào)用數(shù)據(jù)
use coor0, clear
rename _Y lan
rename _X lon
# 將面和線狀文件的邊界坐標(biāo)數(shù)據(jù)進(jìn)行合并
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

區(qū)域上色

如果想要分級(jí)顯示，則需要分別進(jìn)行繪制，這是這樣方式繪圖的不足之處。但是優(yōu)點(diǎn)在于，對(duì)于各個(gè)圖層的處理有很大的自由度。演示如下：

# 導(dǎo)入屬性表
use data_shp, clear
# 與我們的其他經(jīng)濟(jì)變量數(shù)據(jù)進(jìn)行合并
merge 1:1 xian_code using data2008
rename id _ID
drop _merge
save shuxingbiao.dta, replace
# 代開地理邊界數(shù)據(jù)
use coor0.dta, clear
# 與以及合并的屬性表再次合并
merge m:1 _ID using shuxingbiao.dta
# 重命名變量，放置與后面所有合并的表沖突
rename _Y lan
rename _X lon
drop _merge
# 合并分界線的地理邊界數(shù)據(jù)
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 轉(zhuǎn)換數(shù)據(jù)格式
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  //點(diǎn)圖
* 地圖繪制
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

往期回顧

1. 如何使用R、matlab、arcgis進(jìn)行趨勢面分析
2. 使用R進(jìn)行空間自相關(guān)檢驗(yàn)
3. 空間權(quán)重矩陣的那些事(八)-球面距離權(quán)重矩陣
4. R中空間權(quán)重矩陣的處理
5. 使用R進(jìn)行地圖可視化
6. 使用QGIS進(jìn)行地圖可視化
7. 南海九段線可以這樣處理
8. Stata 繪圖也可以很好看

image