使用 Stata 绘制中国地图+空间网络图

前不久给大家分享过「上市公司前5大供应商和客户工商注册数据匹配结果（含经纬度和所处的省市区县）」数据，在数据介绍中我展示了「2022 年上市公司前 5 大客户与供应商地理分布」：

这幅图是使用 Stata 绘制的，今天的课程中我们将一起学习下该如何在 Stata 中绘制。

在学习该课程前，需要预先学习下面两个课程：

×

今天的课程也将在之前的这两个课程的基础上进行讲解。

绘制地图和在地图上添加散点图的方法之前我们都已经讲解过了，因此我们今天就直接来看如何在地图上添加连接线。

运行 help spmap 可以看到连接线可以使用 arrow() 选项添加：

然后在下面可以看到 arrow data 的格式是这样的：

_ID 是每个连接线的编号、两组 _X 和 _Y 显然是连接线首尾的坐标，byvar_ar 则是连接线的分组变量。

下面我们就来构造这个数据：

*- 供应商连接线
use 2022年上市公司供应商信息.dta, clear

keep 经度 纬度 省 股票代码
ren 经度 供应商经度
ren 纬度 供应商纬度
ren 省 供应商所在省
merge m:1 股票代码 using "上市公司注册地址经纬度数据.dta"
keep if _m == 3
keep 股票代码 供应商经度 供应商纬度 供应商所在省 注册地址经度 注册地址纬度 注册地址所在省份
ren 注册地址经度 上市公司经度
ren 注册地址纬度 上市公司纬度
ren 注册地址所在省份 上市公司所在省
gen _ID = _n
order _ID
geo2xy 供应商纬度 供应商经度, gen(_Y1 _X1) projection(albers, 6378137 298.257223563 25 47 0 105) replace
geo2xy 上市公司纬度 上市公司经度, gen(_Y2 _X2) projection(albers, 6378137 298.257223563 25 47 0 105) replace
encode 上市公司所在省, gen(byvar)

keep _* byvar 股票代码
gsort 股票代码

save arrowdata, replace

为了避免图不是过于单调，我们可以再把散点图层加上去，下面构造散点图所需的数据：

use arrowdata, clear
keep _X1 _Y1
ren _X1 lon
ren _Y1 lat
gen class = "供应商"
save tempdata1, replace

use arrowdata, clear
keep _X2 _Y2
ren _X2 lon
ren _Y2 lat
duplicates drop _all, force
gen class = "上市公司"
append using tempdata1
encode class, gen(class_group)
save "mainpoint3", replace

*- 提取小地图框格的
use mainpoint3, clear
keep if inrange(lon, 120000, 1766004.1) & inrange(lat, 320000, 2557786.0)
replace lon = lon * 0.5 + 2100000
replace lat = lat * 0.5 + 1665139
replace class = "smallbox"
append using mainpoint3
save pointdata3, replace

然后就可以绘制地图了，由于细节较多，这里就不再一步步演示了，详情可以观看视频讲解学习：

use chinaprov2021mini_db, clear
local colorlist = `""80 80 255" "206 61 50" "116 155 88" "240 230 133" "70 105 131" "186 99 56" "93 177 221" "128 34 104" "107 215 107" "213 149 167" "146 72 34" "131 123 141" "199 81 39" "213 143 92" "122 101 165" "228 175 105" "59 27 83" "205 222 183" "97 42 121" "174 31 99" "231 199 111" "90 101 94" "204 153 0" "153 204 0" "169 169 169" "204 153 0" "153 204 0" "51 204 0" "0 204 51" "0 204 153" "0 153 204" "10 71 255" "71 117 255" "255 194 10""'
grmap using chinaprov2021mini_coord.dta, ///
  id(ID) osize(vthin ...) ocolor(gray ...) ///
  graphr(margin(medium)) ///
  fcolor("white" ...) ///
  line(data(line_with_nhb) by(group) ///
    select(drop if inlist(group, 4, 7)) ///
    size(vvthin *1.2 *0.5 *0.5 *0.5 *1.2) pattern(solid ...) ///
    color(white /// 省界颜色
        "162 154 196" /// 国界线颜色
        "0 85 170" /// 海岸线颜色
        black /// 小地图框格颜色
        black /// 比例尺和指北针颜色
        )) ///
  polygon(data(polygon_with_nhb) by(class_group) ///
    fcolor(black black "237 237 237") ///
    osize(vvthin ...) ocolor(black black black)) ///
  label(data(label_with_nhb) x(X) y(Y) label(cname) ///
    length(20) size(*0.8) ///
    select(drop if inlist(cname, "乌兹别克斯坦", "塔吉克斯坦", ///
      "阿富汗", "巴基斯坦", "锡亚琛冰川", "马来西亚", "柬埔寨", ///
      "印度尼西亚", "文莱"))) ///
  ti("前五大供应商") ///
  plotr(fcolor("187 209 235") margin(-0.95 -0.95 -0.95 -0.95)) ///
  point(data(pointdata3) by(class_group) ///
      fcolor("196 0 3" "0 193 155") ///
        x(lon) y(lat) ///
        size(tiny ...) legenda(off)) leg(off) ///
  arrow(data(arrowdata) by(byvar) ///
        lcolor(`colorlist') ///
        lpattern(solid ...) lsize(vvthin ...) ///
        hfcolor(`colorlist') hocolor(`colorlist') ///
        hosize(vvthin ...) hsize(small ...) hbarbsize(small ...)) ///
  name(c, replace) nodraw

gr display c

上市公司与客户的连接线绘制方法类似：

use 2022年上市公司客户信息.dta, clear

keep 经度 纬度 省 股票代码
ren 经度 客户经度
ren 纬度 客户纬度
ren 省 客户所在省
merge m:1 股票代码 using "上市公司注册地址经纬度数据.dta"
keep if _m == 3
keep 股票代码 客户经度 客户纬度 客户所在省 注册地址经度 注册地址纬度 注册地址所在省份
ren 注册地址经度 上市公司经度
ren 注册地址纬度 上市公司纬度
ren 注册地址所在省份 上市公司所在省
gen _ID = _n
order _ID
geo2xy 客户纬度 客户经度, gen(_Y2 _X2) projection(albers, 6378137 298.257223563 25 47 0 105) replace
geo2xy 上市公司纬度 上市公司经度, gen(_Y1 _X1) projection(albers, 6378137 298.257223563 25 47 0 105) replace
encode 上市公司所在省, gen(byvar)

keep _* byvar 股票代码
gsort 股票代码

save arrowdata2, replace

use arrowdata2, clear

* 2021 年的上市公司和客户
use arrowdata2, clear
keep _X2 _Y2
ren _X2 lon
ren _Y2 lat
gen class = "客户"
save tempdata2, replace

use arrowdata2, clear
keep _X1 _Y1
ren _X1 lon
ren _Y1 lat
duplicates drop _all, force
gen class = "上市公司"
append using tempdata1
encode class, gen(class_group)
save "mainpoint4", replace

* 提取小地图框格的
use mainpoint4, clear
keep if inrange(lon, 120000, 1766004.1) & inrange(lat, 320000, 2557786.0)
replace lon = lon * 0.5 + 2100000
replace lat = lat * 0.5 + 1665139
replace class = "smallbox"
append using mainpoint4
save pointdata4, replace

use chinaprov2021mini_db, clear
local colorlist = `""80 80 255" "206 61 50" "116 155 88" "240 230 133" "70 105 131" "186 99 56" "93 177 221" "128 34 104" "107 215 107" "213 149 167" "146 72 34" "131 123 141" "199 81 39" "213 143 92" "122 101 165" "228 175 105" "59 27 83" "205 222 183" "97 42 121" "174 31 99" "231 199 111" "90 101 94" "204 153 0" "153 204 0" "169 169 169" "204 153 0" "153 204 0" "51 204 0" "0 204 51" "0 204 153" "0 153 204" "10 71 255" "71 117 255" "255 194 10""'
grmap using chinaprov2021mini_coord.dta, ///
  id(ID) osize(vthin ...) ocolor(gray ...) ///
  graphr(margin(medium)) ///
  fcolor("white" ...) ///
  line(data(line_with_nhb) by(group) ///
    select(drop if inlist(group, 4, 7)) ///
    size(vvthin *1.2 *0.5 *0.5 *0.5 *1.2) pattern(solid ...) ///
    color(white /// 省界颜色
        "162 154 196" /// 国界线颜色
        "0 85 170" /// 海岸线颜色
        black /// 小地图框格颜色
        black /// 比例尺和指北针颜色
        )) ///
  polygon(data(polygon_with_nhb) by(class_group) ///
    fcolor(black black "237 237 237") ///
    osize(vvthin ...) ocolor(black black black)) ///
  label(data(label_with_nhb) x(X) y(Y) label(cname) ///
    length(20) size(*0.8) ///
    select(drop if inlist(cname, "乌兹别克斯坦", "塔吉克斯坦", ///
      "阿富汗", "巴基斯坦", "锡亚琛冰川", "马来西亚", "柬埔寨", ///
      "印度尼西亚", "文莱"))) ///
  ti("前五大客户") ///
  plotr(fcolor("187 209 235") margin(-0.95 -0.95 -0.95 -0.95)) ///
  point(data(pointdata4) by(class_group) ///
      fcolor("196 0 3" "0 193 155") ///
        x(lon) y(lat) ///
        size(tiny ...) legenda(on)) ///
  arrow(data(arrowdata2) by(byvar) ///
        lcolor(`colorlist') ///
        lpattern(solid ...) lsize(vvthin ...) ///
        hfcolor(`colorlist') hocolor(`colorlist') ///
        hosize(vvthin ...) hsize(small ...) hbarbsize(small ...)) ///
  leg(order(11 "上市公司" 12 "客户/供应商") row(1) size(*1.5)) ///
  name(d, replace) nodraw

最后我们再使用 grc1leg2 把两个图合并起来并使用统一的图例：

* net install grc1leg2.pkg, from(http://digital.cgdev.org/doc/stata/MO/Misc)

grc1leg2 c d, row(1) legendfrom(d) xsize(10) ysize(6) ///
  ti("2022 年上市公司前 5 大客户与供应商地理分布", bexpand justification(center) size(*1.2)) ///
  subti("数据处理 & 绘制：微信公众号 RStata", bexpand justification(center) size(*1.2)) ///
  caption("数据来源：国泰安数据库、天眼查，使用高德地图地理编码解析经纬度")
gr export pic.png, width(4800) replace

这样就绘制好了这幅地图。

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制中国地图+空间网络图