之前给大家分享了使用 R 语言绘制城市间专利合作申请数量网络图的方法：
使用 R 语言绘制城市间专利合作申请数量网络图(一)：https://rstata.duanshu.com/#/brief/course/4ef23aa8632046e283b8e7175f6a6f7f
使用 R 语言绘制城市间专利合作申请数量网络图(二)：https://rstata.duanshu.com/#/brief/course/cc8ff6330fd944039f1af4c2d6622342
考虑到有小伙伴没有使用过 R 语言，今天我们再来讲解下如何使用 Stata 绘制这样的图表：
实际上之前也有讲解过类似的课程：
使用 Stata 绘制中国地图+空间网络图：https://rstata.duanshu.com/#/brief/course/06ef5cad0db54403ba93916e3b11c404
使用 Stata 绘制地图是一门精细的操作，因此在学习本课程之前，需要预先学习下面的课程：
使用 Stata 绘制历年中国市级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/0e9d78633ae64fefa684d6aad89633bc
构造 arrow 数据
网络图里面的每条边都是从起点到终点的一条线，但是我们的数据里面并没有起点和终点的坐标：
use 2020年城市间各类型专利合作数量统计.dta, clear
list in 1/10
|
因此第一步我们就需要在数据里面添加每条边的起点和终点坐标，也就是各个城市的质心：
shp2dta using 2021行政区划/市.shp, database(city_db) coordinates(city_coord) genid(ID) replace gencentroids(centroids)
use city_db, clear
keep 市 *cen*
geo2xy y_centroids x_centroids, gen(_Y1 _X1) projection(albers, 6378137 298.257223563 25 47 0 105) replace
drop *cen*
ren 市 城市1
save 各城市质心坐标, replace
|
由于提供的地图数据是 aea 坐标系的，所以所有的坐标也都需要转换成 aea 坐标系的。
然后就可以分别把起点和终点的坐标添加到数据里面了：
use 2020年城市间各类型专利合作数量统计.dta, clear
drop if 城市1 == "三沙市" | 城市2 == "三沙市"
merge m:1 城市1 using 各城市质心坐标
drop _m
ren (_X1 _Y1) (_X2 _Y2)
ren 城市1 to
ren 城市2 城市1
merge m:1 城市1 using 各城市质心坐标
drop _m
ren 城市1 from
keep from to _* 合作申请专利数量
ren 合作申请专利数量 value
drop if mi(_X1) | mi(_X2)
bysort from: egen sum = sum(value)
egen rank = rank(-sum)
gsort rank
egen group1 = group(rank)
replace group1 = 11 if group1 > 10
drop rank
egen group2 = cut(value), group(6)
tab group2
replace group2 = 5 if group1 <= 10
gsort group1 group2
egen group = group(group1 group2)
gsort group1
gen _ID = _n
save arrowdata, replace
list in 1/10
|
这里需要注意 arrow 数据的格式要求，运行 help spmap##spatdata：
_ID _X1 _Y1 _X2 _Y2 byvar_ar
---------------------------------------------------------
1 11 30 18 30 1
2 15 40 15 45 1
3 15 40 25 40 1
4 20 35 28 45 2
5 17 20 20 11 2
---------------------------------------------------------
|
这也是为什么我们要把 arrowdata 构造成前面的样子。
构造 point 数据
point 数据是用来绘制一些关键节点城市：
此处代码需下载讲义材料查看～
可以看到 point 数据的格式也要按照要求。
构造 label 数据
除了指北针和比例尺之外，我们还希望添加 value 和最大的 10 个城市的标签：
此处代码需下载讲义材料查看～
最后我们就可以绘图了，由于细节较多，这里就不在一句句解释了，建议结合视频讲解学习：
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""'
local rawcolor = `""0 115 194" "239 192 0" "134 134 134" "205 83 76" "122 166 220" "0 60 103" "143 119 0" "59 59 59" "167 48 48" "74 105 144" gray"'
local colorlist2 = ""
forval i = 1/11 {
local tempcolor: word `i' of `rawcolor'
forval j = 1/5 {
local colorlist2 = `"`colorlist2' "`tempcolor'""'
}
}
local rawsize = `"vvvthin vvvthin vvvthin vvvthin vvthin"'
local sizelist2 = ""
forval i = 1/11 {
local sizelist2 = `"`sizelist2' `rawsize'"'
}
spmap 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(mylabel) x(X) y(Y) label(cname)
length(20) size(*0.8))
ti("2020 年城市间专利合作网络")
subti("数据处理 & 绘制：微信公众号 RStata", size(*0.8))
caption("数据来源：国家知识产权局", size(*0.8))
plotr(fcolor("187 209 235") margin(-0.95 -0.95 -0.95 -0.95))
arrow(data(arrowdata) by(group)
lcolor(`colorlist2')
lpattern(solid ...) lsize(`sizelist2')
hfcolor(`colorlist2') hocolor(`colorlist2')
hosize(vvvthin ...) hsize(tiny ...) hbarbsize(tiny ...))
point(data(first100df) by(byvar)
fcolor("black" ...)
x(_X) y(_Y)
size(vtiny tiny vsmall small medsmall) legenda(off))
leg(off)
gr export pic1.png, width(4800) replace
|
把 point 图层放置在上层
由于没有办法解决散点被压在下面的问题，所以我们可以考虑另外一种绘制线条的方案：
use arrowdata, clear
expand 3
gsort _ID
bysort _ID: replace _X1 = _X2 if _n == 2
bysort _ID: replace _Y1 = _Y2 if _n == 2
drop _X2 _Y2
ren (_X1 _Y1) (_X _Y)
bysort _ID: replace _X = . if _n == _N
bysort _ID: replace _Y = . if _n == _N
replace _ID = _ID + 43
replace group = group + 7
append using line_with_nhb
save myline, replace
> 此处代码需下载讲义材料查看～
|
去除邻国地区
邻国区域主要是包括 polygon 和 label 数据，因此从二者里面删除相应的内容即可：
use mylabel, clear
drop if !mi(ename) & !inlist(ename, "N", "1000km")
replace Y = Y + 80000 if mi(ename)
save mylabel2, replace
use polygon_with_nhb, clear
drop if class == "邻国"
save polygon_with_nhb2, 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""'
local rawcolor = `""0 115 194" "239 192 0" "134 134 134" "205 83 76" "122 166 220" "0 60 103" "143 119 0" "59 59 59" "167 48 48" "74 105 144" gray"'
local colorlist2 = ""
forval i = 1/11 {
local tempcolor: word `i' of `rawcolor'
forval j = 1/5 {
local colorlist2 = `"`colorlist2' "`tempcolor'""'
}
}
local rawsize = `"vvvthin vvvthin vvvthin vvvthin vvthin"'
local sizelist2 = ""
forval i = 1/11 {
local sizelist2 = `"`sizelist2' `rawsize'"'
}
spmap using chinaprov2021mini_coord.dta,
id(ID) osize(vthin ...) ocolor(gray ...)
graphr(margin(medium))
fcolor("white" ...)
line(data(myline) by(group)
select(drop if inlist(group, 4, 7))
size(vvthin *1.2 *0.5 *0.5 *0.5 `sizelist2')
pattern(solid ...)
color(white
"162 154 196"
"0 85 170"
black
black
`colorlist2'))
polygon(data(polygon_with_nhb2) by(class_group)
fcolor(black black "237 237 237")
osize(vvthin ...) ocolor(black black black))
label(data(mylabel2) x(X) y(Y) label(cname)
length(20) size(*0.8))
ti("2020 年城市间专利合作网络")
subti("数据处理 & 绘制：微信公众号 RStata", size(*0.8))
caption("数据来源：国家知识产权局", size(*0.8))
plotr(margin(-0.95 -0.95 -0.95 -0.95))
point(data(first100df) by(byvar)
fcolor("black%80" ...)
x(_X) y(_Y)
size(vtiny tiny vsmall small medsmall) legenda(off))
leg(off)
gr export pic3.png, width(4800) replace
|
点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制城市间专利合作申请数量网络图
评论