使用 Stata 绘制中国地图与世界地图的网络流向图

今天给大家讲解下如何使用 Stata 绘制中国地图与世界地图的网络流向图的方法。

我准备了两个版本，一个是竖版：

另外一个是横版：

下面我们会以竖版的为例进行讲解。横版类似，感兴趣的小伙伴可以结合本课程自行学习。

本课程是在下面这些课程的基础上开展的，需要预先学习下面几个课程：

1. 使用 Stata 绘制历年中国省级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/e194d75fe1674f7c8921daec1ece7f7d
2. 使用 Stata 绘制历年中国市级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/0e9d78633ae64fefa684d6aad89633bc
3. 使用 Stata 绘制历年中国县级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/99c2e97b88ea401bbdb8da4511341cc8
4. 使用 Stata 绘制以太平洋为中心的世界地图（带指北针和比例尺）：https://rstata.duanshu.com/#/brief/course/f995d70c54094c968611eeeb73b25d76

使用 R 语言创建中国地图 + 世界地图的矢量数据

这部分需要使用 R 语言，不过不会使用也没关系，可以直接跳过这个部分。后面如果有需要特别的地图数据，可以联系李老师定制。

首先加载所需的 R 包：

library(tidyverse)
library(sf)

读取经过编辑后的世界地图：

read_sf("worldmap/worldmap.shp") -> worldmap
read_sf("worldmap/worldmap_line.shp") %>%
  filter(is.na(continent)) %>%
  select(name = country) -> worldmapline

st_crs(worldmap)$proj4string

注意这份世界地图的坐标系是“+proj=eck3 +lon_0=150 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs”。

再读取 aea 坐标系下的省级行政区划数据，这份数据来源上面的使用 Stata 绘制省级地图的课程里面：

mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

read_sf("chinaprov2021mini/chinaprov2021mini.shp") %>%
  select(-contains("类型")) %>%
  filter(省代码 != 0) -> prov

read_sf("chinaprov2021mini/chinaprov2021mini_line.shp") %>%
  select(-contains("类型")) %>%
  filter(省代码 == 0) -> provline

直接合并两份地图数据，后面再在 Stata 中进行地图数据的调整。

# 合并
bind_rows(
  worldmap %>% mutate(class = "世界地图"),
  worldmapline %>% mutate(class = "世界地图线条"),
  prov %>% mutate(class = "省级地图") %>% st_set_crs(st_crs(worldmap)),
  provline %>% mutate(class = "省级地图线条") %>% st_set_crs(st_crs(worldmap))
) -> dfall

# 保存为 shp 文件
dir.create("prov+worldmap_long")
dfall %>%
  st_make_valid() %>%
  st_collection_extract("POLYGON") %>%
  st_write("prov+worldmap_long/prov+worldmap_long.shp", layer_options = "ENCODING=UTF-8", delete_layer = TRUE, layer = "MULTIPOLYGON")

dfall %>%
  st_make_valid() %>%
  st_cast("MULTILINESTRING") %>%
  st_write("prov+worldmap_long/prov+worldmap_long_line.shp", layer_options = "ENCODING=UTF-8", delete_layer = TRUE, layer = "MULTIPOLYGON")

各国中心的坐标及各省的中心坐标经常会被用到，所以也保存下：

worldmap %>%
  filter(!is.na(continent)) %>%
  st_point_on_surface() -> worldmappoint

bind_cols(
  worldmappoint %>% st_drop_geometry(),
  worldmappoint %>% st_coordinates()
) %>%
  rename(value = iso3) %>%
  select(value, X, Y) %>%
  haven::write_dta("世界各国坐标.dta")

read_sf("chinaprov2021long/chinaprov2021long.shp") %>%
  select(-contains("类型")) %>%
  filter(省代码 != 0) -> prov

bind_cols(
  prov %>% st_drop_geometry(),
  prov %>% st_point_on_surface() %>%
    st_coordinates()
) %>%
  mutate(省 = str_sub(省, 1, 2)) %>%
  haven::write_dta("各省坐标.dta")

这样就完成了 R 语言部分的内容，简单来说就是直接把两种数据堆起来保存成 shp 文件。

使用 Stata 进行地图数据的调整

下面我们再使用 Stata 进行地图数据的调整。

首先把 shp 文件转换成 dta 文件：

local name = "prov+worldmap_long"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) replace
shp2dta using `name'/`name'_line, database(`name'_line_db) coordinates(`name'_line_coord) genid(ID) replace

*- 概览
use prov+worldmap_long_db.dta, clear
grmap using prov+worldmap_long_coord.dta, id(ID) freestyle

世界地图的比例尺和指北针使用 polygon 绘制：

use prov+worldmap_long_coord.dta, clear
keep if inlist(_ID, 239, 240)
gen value = 1
save worldmap_long_polygon, replace

平移中国地图：

use prov+worldmap_long_db.dta, clear
*- 可以看到中国省级地图的 ID 范围是 >= 240
use prov+worldmap_long_coord.dta, clear
replace _X = _X + 0 if _ID >= 241
replace _Y = _Y + 800000 if _ID >= 241
replace _X = _X * 4 if _ID >= 241
replace _Y = _Y * 4 if _ID >= 241

*- 删除比例尺和指北针
drop if inlist(_ID, 239, 240)
save prov+worldmap_long_coord2.dta, replace

再次概览：

use prov+worldmap_long_db.dta, clear
grmap using prov+worldmap_long_coord2.dta, id(ID)

线条数据也要进行同步修改：

use prov+worldmap_long_line_db.dta, clear
keep if index(class, "线条")
*- 删除中国地图的比例尺和指北针
keep if mi(省) | inlist(省, "国界线", "海岸线", "小地图框格")
replace name = 省 if mi(name)
keep class name ID
save prov+worldmap_long_line_db2, replace

use prov+worldmap_long_line_coord.dta, clear
keep if inlist(_ID, 241, 242, 243, 244, 245, 280, 281, 283)

*- 国界线和海岸线作相同的变换：
replace _X = _X + 0 if _ID >= 280
replace _Y = _Y + 800000 if _ID >= 280
replace _X = _X * 4 if _ID >= 280
replace _Y = _Y * 4 if _ID >= 280

*- 生成类别变量
gen classgroup = 1 if _ID == 241 // 世界地图上的九段线和国界线
replace classgroup = 2 if inlist(_ID, 242, 243, 244, 245) // 世界地图的指北针和比例尺
replace classgroup = 3 if inlist(_ID, 280) // 中国地图的国界线和九段线
replace classgroup = 4 if inlist(_ID, 281) // 中国地图的海岸线
replace classgroup = 5 if inlist(_ID, 283) // 中国地图的小地图框格
save prov+worldmap_long_line_coord2.dta, replace

再次概览：

use prov+worldmap_long_db.dta, clear
grmap using prov+worldmap_long_coord2.dta, id(ID) ///
 line(data(prov+worldmap_long_line_coord2.dta) by(classgroup) ///
  color(black black "162 154 196" "0 85 170" black) ///
  size(thin thin *1.2 *0.5 *0.5)) ///
 polygon(data(worldmap_long_polygon) fcolor(black ...) osize(thin)) ///
 graphr(margin(medium))

5000km 和 N 向上移动一些：

use worldmap_label2, clear
replace Y = Y + 200000 if inlist(country_cn, "N", "5000km")
save worldmap_label3, replace

这样就调整好了。

处理绘图数据

这里我们使用 2006年海关进出口数据汇总结果.dta 演示。

首先我们需要获取每个消费地或生产地的经纬度，可以使用地理编码：

这部分内容可以详细学习：

使用 Stata 进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县：https://rstata.duanshu.com/#/brief/course/537300af1a9947edb758789785c600f3

use 2006年海关进出口数据汇总结果.dta, clear

*- 消费地或生产地使用地理编码获得经纬度
keep 消费地或生产地
duplicates drop 消费地或生产地, force
gen address = subinstr(消费地或生产地, "其它", "", .)
gen addid = _n
save "待解析经纬度", replace

*- 解析经纬度
use 待解析经纬度, clear
*- 创建一个文件夹保存
*- 这里的密钥需要替换成自己的
cap mkdir "json1"
forval i = 1/`=_N' {
 di "`i'"
 qui {
  if !fileexists("json1/`=addid[`i']'.json") {
   percentencode `=address[`i']'
   local a = r(percentencode)

   local url = "https://restapi.amap.com/v3/geocode/geo?address=`a'&output=JSON&key=ad13679a4baac443f5a63268ce0d7b37"
   copy "`url'" "json1/`=addid[`i']'.json", replace
  }
 }
}

*- 处理所有的 json 文件
clear
gen str100 location = ""
gen str100 addid = ""

*- 循环
local files: dir "json1" files "*.json"
di `"`files'"'
local n: word count `files'
di "`n'"

foreach i in `files' {
 insheetjson location using "json1/`i'", table("geocodes") ///
  columns("location") ///
  offset(`=_N')
 replace addid = "`i'" if mi(addid)
}

compress
foreach i of varlist _all {
 cap format `i' %10s
}
order addid
replace addid = subinstr(addid, ".json", "", .)
destring addid, replace
merge m:1 addid using 待解析经纬度.dta
keep if _m == 3
drop _m
keep 消费地或生产地 location
split location, parse(,)
drop location
ren location1 经度
ren location2 纬度

destring, replace

*- 转 WGS84 坐标系
gcj02towgs84, lon(经度) lat(纬度)
drop 经度 纬度
ren 经度_WGS84 经度
ren 纬度_WGS84 纬度
duplicates drop 消费地或生产地, force
save 消费地或生产地的经纬度, replace

起运国或目的国的经纬度就没办法使用地理编码获取了，只能先和 worldmap_label2.dta 匹配，匹配不上的再一一手动调整：

*- 起运国或目的国的经纬度
use 2006年海关进出口数据汇总结果.dta, clear
keep 起运国或目的国
duplicates drop 起运国或目的国, force
drop if inlist(起运国或目的国, "国别（地区）不详", "中华人民共和国")
drop if index(起运国或目的国, "其他国家")
save 待获取起运国或目的国的经纬度, replace

*- 和标签数据匹配
use worldmap_label2.dta, clear

use 待获取起运国或目的国的经纬度, clear
gen country_cn = 起运国或目的国
replace country_cn = subinstr(country_cn, "共和国", "", .)
replace country_cn = subinstr(country_cn, "联邦", "", .)
replace country_cn = "几内亚比绍" if country_cn == "几内亚(比绍)"
replace country_cn = "前南马其顿" if country_cn == "前南斯拉夫马其顿"
replace country_cn = "-17.5,28.5" if country_cn == "加那利群岛"
replace country_cn = "20.37,44.50" if country_cn == "南斯拉夫联盟"
replace country_cn = "-68.2385339,12.1783611" if country_cn == "博内尔"
replace country_cn = "博茨瓦纳" if country_cn == "博茨瓦那"
replace country_cn = "121.509064,25.044333" if country_cn == "台湾省"
replace country_cn = "吉尔吉斯斯坦" if country_cn == "吉尔吉斯"
replace country_cn = "179.13,-8.31" if country_cn == "图瓦卢"
replace country_cn = "-149.5,-23.5" if country_cn == "土布艾群岛"
replace country_cn = "圣基茨和尼维斯" if country_cn == "圣其茨--尼维斯"
replace country_cn = "法属圣马丁" if country_cn == "圣马丁岛"
replace country_cn = "摩洛哥" if country_cn == "塞卜泰(休达)"
replace country_cn = "塞舌尔群岛" if country_cn == "塞舌尔"
replace country_cn = "密克罗尼西亚联邦" if country_cn == "密克罗尼西亚"
replace country_cn = "-68.990021,12.169570" if country_cn == "库腊索岛"
replace country_cn = "27.5,-11.12" if country_cn == "扎伊尔"
replace country_cn = "格陵兰（丹）" if country_cn == "格陵兰"
replace country_cn = "35.2922775,-2.9380973" if country_cn == "梅利利亚"
replace country_cn = "梵蒂冈" if country_cn == "梵蒂冈城国"
replace country_cn = "圭亚那" if country_cn == "法属圭亚那"
replace country_cn = "-149.5,-17.5" if country_cn == "法属波利尼西亚"
replace country_cn = "17.807894,43.34" if country_cn == "波斯尼亚-黑塞哥维那"
replace country_cn = "17.807894,43.34" if country_cn == "波斯尼亚 黑塞哥维那"
replace country_cn = "113.5,22.2" if country_cn == "澳门"
replace country_cn = "-61.5,16.23" if country_cn == "瓜德罗普岛"
replace country_cn = "瓦利斯群岛和富图纳群岛" if country_cn == "瓦利斯和浮图纳"
replace country_cn = "55.5,-21.1" if country_cn == "留尼汪"
replace country_cn = "-134.9,-23.1" if country_cn == "盖比群岛"
replace country_cn = "-5.22,36.09" if country_cn == "直布罗陀"
replace country_cn = "-149.5,-17.5" if country_cn == "社会群岛"
replace country_cn = "-69.9,12.43" if country_cn == "荷属安地列斯群岛"
replace country_cn = "-4.13,50.36" if country_cn == "蒙特塞拉特"
replace country_cn = "撒哈拉西部" if country_cn == "西撒哈拉"
replace country_cn = "帕劳" if country_cn == "贝劳"
replace country_cn = "阿联酋" if country_cn == "阿拉伯联合酋长国"
replace country_cn = "阿鲁巴" if country_cn == "阿鲁巴岛"
replace country_cn = "114.25,22.25" if country_cn == "香港"
replace country_cn = "-149.5,-17.5" if country_cn == "马克萨斯群岛"
replace country_cn = "-61.64,14.6" if country_cn == "马提尼克岛"
replace country_cn = "45.167,-12.98" if country_cn == "马约特岛"
replace country_cn = "-162,61.9" if country_cn == "马绍尔群岛"

merge m:1 country_cn using worldmap_label2.dta
drop if _m == 2
keep 起运国或目的国 country_cn X Y
split country_cn, parse(,)
destring country_cn*, replace force
save temp, replace

keep country_cn*
ren country_cn1 lon
ren country_cn2 lat
drop if mi(lon)
export delimited using 待上传.csv, replace
*- 使用 https://czxb.shinyapps.io/crs-trans/ 转换坐标系，crs 要设置成 +proj=eck3 +lat_0=0 +lon_0=150 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

import delimited using "转换后的数据.csv", clear
gen id = _n
save temp2, replace

use temp, clear
keep country_cn*
ren country_cn1 lon
ren country_cn2 lat
drop if mi(lon)
gen id = _n
merge 1:1 id using temp2
drop _m
duplicates drop lon lat, force
save 待补充的坐标信息, replace

use temp, clear
ren country_cn1 lon
ren country_cn2 lat
merge m:1 lon lat using 待补充的坐标信息
drop _m
replace X = x if mi(X)
replace Y = y if mi(Y)
keep 起运国或目的国 X Y
save 起运国或目的国的经纬度结果, replace

然后就可以把坐标都匹配上，然后处理成所需要的样子了：

use 2006年海关进出口数据汇总结果.dta, clear
*- 出口部分
keep if 进口或出口 == "出口"
drop if 消费地或生产地 == "NULL"
merge m:1 消费地或生产地 using 消费地或生产地的经纬度.dta

*- 根据地图设计，省份的经纬度需要进行如下的变换：
geo2xy 纬度 经度, gen(_Y1 _X1) projection(albers, 6378137 298.257223563 25 47 0 105) replace
replace _X1 = _X1 + 0
replace _Y1 = _Y1 + 800000
replace _X1 = _X1 * 4
replace _Y1 = _Y1 * 4
drop *度 _m
merge m:1 起运国或目的国 using 起运国或目的国的经纬度结果.dta
drop _m
drop if mi(起运国或目的国)
ren X _X2
ren Y _Y2
gen _ID = _n
drop if mi(_X1)
save arrowdata, replace

*- 保留最大的 300 条、数据分组
use arrowdata, clear
gsort -金额
keep in 1/300
replace 金额 = 金额/100000000
egen group1 = cut(金额), group(8) label
save arrowdata2, replace

这就是 arrowdata 需要的样子。

然后就可以绘制地图+网络图了：

* 绘制地图
use prov+worldmap_long_db.dta, clear
local colorlist = `""128 203 196" "128 203 196" "128 203 196" "77 182 172" "38 166 154" "0 150 136" "0 137 123" "0 121 107""'
grmap using prov+worldmap_long_coord2.dta, id(ID) fcolor(gs13 ...) ///
 osize(vthin ...) ///
 line(data(prov+worldmap_long_line_coord2.dta) by(classgroup) ///
  color(black black "162 154 196" "0 85 170") ///
  size(vvthin vthin *1.2 *0.4)) ///
 polygon(data(worldmap_long_polygon) fcolor(black ...) osize(vthin)) ///
 graphr(margin(medium)) ///
 arrow(data(arrowdata2) by(group1) ///
  lcolor(`colorlist') ///
  lpattern(solid ...) ///
  lsize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hfcolor(`colorlist') hocolor(`colorlist') ///
  hosize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hsize(medium ...) hbarbsize(medium ...) legenda(on)) ///
 label(data(worldmap_label3) x(X) y(Y) ///
  color(black) size(*0.8) ///
  label(country_cn) length(80 ...) ///
  select(keep if inlist(country_cn, "N", "5000km"))) ///
 leg(order(9 "<=5.08" 10 "5.08~5.95" 11 "5.95~6.92" 12 "6.92~8.08" ///
  13 "8.08~10.78" 14 "10.78~14.07" 15 "14.07~24.59" 16 ">24.59") ///
  pos(6) row(3) ring(1) size(*1) ti("出口金额（亿美元）", size(*0.4) pos(9))) ///
 ti("2006 年中国各地部分出口方向及金额", size(*0.7)) ///
 subti("数据计算 & 绘制：微信公众号 RStata，仅展示总金额排名前 300 的流向", size(*0.7)) ///
 caption("数据来源：2006 年中国海关数据库", size(*0.7)) ///
 xsize(10) ysize(12)

gr export pic1.png, width(4800) replace

进口数据绘制方法类似：

*- 再处理进口的
use 2006年海关进出口数据汇总结果.dta, clear
*- 进口部分
keep if 进口或出口 == "进口"
drop if 消费地或生产地 == "NULL"
merge m:1 消费地或生产地 using 消费地或生产地的经纬度.dta

*- 根据地图设计，省份的经纬度需要进行如下的变换：
geo2xy 纬度 经度, gen(_Y2 _X2) projection(albers, 6378137 298.257223563 25 47 0 105) replace
replace _X2 = _X2 + 0
replace _Y2 = _Y2 + 800000
replace _X2 = _X2 * 4
replace _Y2 = _Y2 * 4
drop *度 _m
merge m:1 起运国或目的国 using 起运国或目的国的经纬度结果.dta
drop _m
drop if mi(起运国或目的国)
ren X _X1
ren Y _Y1
gen _ID = _n
drop if mi(_X1)
save arrowdata_b, replace

*- 保留最大的 300 条、数据分组
use arrowdata_b, clear
gsort -金额
keep in 1/300
replace 金额 = 金额/100000000
egen group1 = cut(金额), group(8) label
save arrowdata2b, replace

tab group1

*- 绘制地图
use prov+worldmap_long_db.dta, clear
local colorlist = `""255 235 238" "255 205 210" "239 154 154" "229 115 115" "239 83 80" "244 67 54" "229 57 53" "211 47 47""'
grmap using prov+worldmap_long_coord2.dta, id(ID) fcolor(gs13 ...) ///
 osize(vthin ...) ///
 line(data(prov+worldmap_long_line_coord2.dta) by(classgroup) ///
  color(black black "162 154 196" "0 85 170") ///
  size(vvthin vthin *1.2 *0.4)) ///
 polygon(data(worldmap_long_polygon) fcolor(black ...) osize(vthin)) ///
 graphr(margin(medium)) ///
 arrow(data(arrowdata2b) by(group1) ///
  lcolor(`colorlist') ///
  lpattern(solid ...) ///
  lsize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hfcolor(`colorlist') hocolor(`colorlist') ///
  hosize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hsize(medium ...) hbarbsize(medium ...) legenda(on)) ///
 label(data(worldmap_label3) x(X) y(Y) ///
  color(black) size(*0.8) ///
  label(country_cn) length(80 ...) ///
  select(keep if inlist(country_cn, "N", "5000km"))) ///
 leg(order(9 "<=5.08" 10 "5.08~5.95" 11 "5.95~6.92" 12 "6.92~8.08" ///
  13 "8.08~10.78" 14 "10.78~14.07" 15 "14.07~24.59" 16 ">24.59") ///
  pos(6) row(3) ring(1) size(*1) ti("进口金额（亿美元）", size(*0.4) pos(9))) ///
 ti("2006 年中国各地部分进口方向及金额", size(*0.7)) ///
 subti("数据计算 & 绘制：微信公众号 RStata，仅展示总金额排名前 300 的流向", size(*0.7)) ///
 caption("数据来源：2006 年中国海关数据库", size(*0.7)) ///
 xsize(10) ysize(12)

也可以把两幅图合并起来：

use prov+worldmap_long_db.dta, clear
local colorlist = `""128 203 196" "128 203 196" "128 203 196" "77 182 172" "38 166 154" "0 150 136" "0 137 123" "0 121 107""'
grmap using prov+worldmap_long_coord2.dta, id(ID) fcolor(gs13 ...) ///
 osize(vthin ...) ///
 line(data(prov+worldmap_long_line_coord2.dta) by(classgroup) ///
  color(black black "162 154 196" "0 85 170") ///
  size(vvthin vthin *1.2 *0.4)) ///
 polygon(data(worldmap_long_polygon) fcolor(black ...) osize(vthin)) ///
 graphr(margin(medium)) ///
 arrow(data(arrowdata2) by(group1) ///
  lcolor(`colorlist') ///
  lpattern(solid ...) ///
  lsize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hfcolor(`colorlist') hocolor(`colorlist') ///
  hosize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hsize(medium ...) hbarbsize(medium ...) legenda(on)) ///
 label(data(worldmap_label3) x(X) y(Y) ///
  color(black) size(*0.8) ///
  label(country_cn) length(80 ...) ///
  select(keep if inlist(country_cn, "N", "5000km"))) ///
 leg(order(9 "<=5.08" 10 "5.08~5.95" 11 "5.95~6.92" 12 "6.92~8.08" ///
  13 "8.08~10.78" 14 "10.78~14.07" 15 "14.07~24.59" 16 ">24.59") ///
  pos(6) row(3) ring(1) size(*1) ti("出口金额（亿美元）", size(*0.4) pos(9))) ///
 ti("出口", size(*0.7)) ///
 xsize(10) ysize(12) name(a, replace) nodraw

use prov+worldmap_long_db.dta, clear
local colorlist = `""255 235 238" "255 205 210" "239 154 154" "229 115 115" "239 83 80" "244 67 54" "229 57 53" "211 47 47""'
grmap using prov+worldmap_long_coord2.dta, id(ID) fcolor(gs13 ...) ///
 osize(vthin ...) ///
 line(data(prov+worldmap_long_line_coord2.dta) by(classgroup) ///
  color(black black "162 154 196" "0 85 170") ///
  size(vvthin vthin *1.2 *0.4)) ///
 polygon(data(worldmap_long_polygon) fcolor(black ...) osize(vthin)) ///
 graphr(margin(medium)) ///
 arrow(data(arrowdata2b) by(group1) ///
  lcolor(`colorlist') ///
  lpattern(solid ...) ///
  lsize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hfcolor(`colorlist') hocolor(`colorlist') ///
  hosize(vvvthin vvvthin vvthin vthin thin medthin medium medthick) ///
  hsize(medium ...) hbarbsize(medium ...) legenda(on)) ///
 label(data(worldmap_label3) x(X) y(Y) ///
  color(black) size(*0.8) ///
  label(country_cn) length(80 ...) ///
  select(keep if inlist(country_cn, "N", "5000km"))) ///
 leg(order(9 "<=5.08" 10 "5.08~5.95" 11 "5.95~6.92" 12 "6.92~8.08" ///
  13 "8.08~10.78" 14 "10.78~14.07" 15 "14.07~24.59" 16 ">24.59") ///
  pos(6) row(3) ring(1) size(*1) ti("进口金额（亿美元）", size(*0.4) pos(9))) ///
 ti("进口", size(*0.7)) ///
 xsize(10) ysize(12) name(b, replace) nodraw

gr combine a b, ti("2006 年中国各地部分进出口方向及金额", ///
  size(*0.9) bexpand justification(center)) ///
 subti("数据计算 & 绘制：微信公众号 RStata，仅展示总金额排名前 300 的流向", size(*0.9) bexpand justification(center)) ///
 caption("数据来源：2006 年中国海关数据库", size(*0.8))

gr export pic3.png, width(4800) replace

这样就解决了这个问题。

中国市级地图 + 世界地图、县级地图 + 世界地图的绘制方法也类似。

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制中国地图与世界地图的网络流向图