之前给大家分享过两个课程：
不过课程里面的南北方区域是手动绘制的，而且是使用 R，对于 Stata 用户 不是很方便，今天再给大家介绍一种新的方案。
Stata 绘制地图中的点线面 Stata 中主要是使用 spmap/grmap 绘制地图，两个命令的用法几乎一样，grmap 是 Stata 15 版本之后的官方命令。
绘制地图需要使用两个 dta 文件，一个是数据表文件、一个是坐标系文件，例如绘制秦岭淮河线：
shp2dta using 2021行政区划/省.shp, database(provdb) coordinates(provcoord) replace genid(ID) use provdb, clear keep if 省 == "湖北省" grmap using provcoord, id(ID)
注意观察 provcoord 的结构：
use provcoord, clear keep if _ID == 13list in 1/5list in `=_N-4'/`=_N'
如果我们直接使用 tw area _Y _X 绘图：
会得到这样的奇奇怪怪的结果。再仔细观察数据，其实里面有一些空行进行分隔：
list in 24/28 *> +-----------------------------+ *> | _ID _X _Y | *> |-----------------------------| *> 24. | 13 109.10203 30.57451 | *> 25. | 13 109.10314 30.574022 | *> 26. | 13 109.10373 30.573526 | *> 27. | 13 . . | *> 28. | 13 109.10532 30.578827 | *> +-----------------------------+
另外可以看到 26 行和第 2 行的坐标值相同：
list if inlist(_n, 2, 26) *> +-----------------------------+ *> | _ID _X _Y | *> |-----------------------------| *> 2. | 13 109.10373 30.573526 | *> 26. | 13 109.10373 30.573526 | *> +-----------------------------+
也就是说 1～26 行的坐标正好可以围成一个闭合的多边形：
所以湖北省的地图是由多个闭合的多边形得到的，而闭合的多边形数据要求第一个坐标和最后一个坐标相同 。这些多边形使用空行分隔，拥有共同的 _ID。
gen class = _n if missing (_X)carryforward class , replace tab class tw area _Y _X if class == 1, color("254 212 57" ) || area _Y _X if class == 27, color("112 154 225" ) || area _Y _X if class == 59, color("138 145 151" ) || area _Y _X if class == 117, color("210 175 129" ) || area _Y _X if class == 17293, color("253 116 70" ) || area _Y _X if class == 17321, color("213 228 162" ) || area _Y _X if class == 17336, color("25 126 192" )
点和线条的数据结构和多边形的类似，不过不要求首尾坐标相同。
处理工企地理位置信息 首先简单处理下工企地理位置信息，去除无关的变量：
use 2014年工企数据库地理位置数据.dta, clear keep gqid 经度 纬度drop if missing (经度)save df, replace list in 1/10
转换 shp 数据生成 dta 文件 再把需要用到的 shp 文件转换成 dta 文件：
shp2dta using 2021行政区划/省.shp, database(provdb) coordinates(provcoord) replace genid(ID) shp2dta using 秦岭淮河线/秦岭淮河线.shp, database(qhdb) coordinates(qhcoord) replace genid(ID) shp2dta using 胡焕庸线/胡焕庸线.shp, database(hhydb) coordinates(hhycoord) replace genid(ID) shp2dta using 九段线/九段线.shp, database(jdxdb) coordinates(jdxcoord) replace genid(ID) shp2dta using 海岸线/海岸线.shp, database(haxdb) coordinates(haxcoord) replace genid(ID) use qhcoord, clear gen class = "秦岭-淮河线" append using hhycoordreplace class = "胡焕庸线" if missing (class )append using jdxcoordreplace class = "九段线" if missing (class )append using haxcoordreplace class = "海岸线" if missing (class )egen newid = group (_ID class )drop _IDren newid _IDencode class , gen (classgroup)codebook classgroupsave "linestring" , replace
绘图展示下：
use provdb, clear grmap using provcoord, id(ID) line (data(linestring) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" ) size(*1.2 *1.2 *2 *2)) graphr(margin(medium))
东西部、南北方的划分标准 通常我们使用胡焕庸线划分东西部，使用秦岭-淮河线划分南北方，不过这两个都是线条元素，不能用于经纬度的判别。所以我们下面需要生成东西部区域、南北方区域的多边形。一种思路就是先生成中国范围的方框，然后延长胡焕庸线、延长秦岭-淮河线把这个方框切分成两个区域。
生成中国范围的方框：
use provcoord, clear drop _IDsum _Xdi "xmin = `r(min)'" di "xmax = `r(max)'" sum _Ydi "ymin = `r(min)'" di "ymax = `r(max)'" clear input _ID _X _Y str40 label 1 . . "" 1 73.502355 53.563624 "topleft" 1 135.09567 53.563624 "topright" 1 135.09567 3.83703 "bottomright" 1 73.502355 3.83703 "bottomleft" 1 73.502355 53.563624 "" end compress save "cnbbox" , replace use provdb, clear grmap using provcoord, id(ID) line (data(linestring) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" ) size(*1.2 *1.2 *2 *2)) graphr(margin(vlarge)) polygon(data(cnbbox) fcolor("none" )) label (data(cnbbox) x(_X) y (_Y) l (label ) length (40))
使用胡焕庸线分东西 选择 cnbbox 的 bottomleft 和 topright 延伸胡焕庸线：
use cnbbox, clear keep if inlist (label , "bottomleft" , "topright" )append using hhycoordgsort _X _Ydrop label gen sort = _nreplace sort = 0 if _n == _Ngsort sort drop sort replace _ID = 5save "longhhyline" , replace use linestring, clear tab _IDappend using longhhylinereplace class = "延长的胡焕庸线" if missing (class )drop classgroupencode class , gen (classgroup)save linestring2, replace use provdb, clear grmap using provcoord, id(ID) line (data(linestring2) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" "196 0 3" ) size(*1.2 *1.2 *2 *2 *2)) graphr(margin(vlarge)) polygon(data(cnbbox) fcolor("none" )) label (data(cnbbox) x(_X) y (_Y) l (label ) length (40))
这样延长的胡焕庸线就可以和 topleft 点构成西部的区域多边形，可以和 bottomright 构成东部的区域多边形：
use cnbbox, clear keep if label == "topleft" save topleft, replace append using longhhylineappend using topleftdrop label replace _ID = 1gen class = "west" drop if missing (_X)save westpart, replace clear input _ID _X _Y str1 class . . . "" end save emptydata, replace use cnbbox, clear keep if label == "bottomright" save bottomright, replace append using longhhylineappend using bottomrightdrop label replace _ID = 2gen class = "east" drop if missing (_X)save eastpart, replace use emptydata, clear replace class = "west" replace _ID = 1append using westpartappend using emptydatareplace _ID = 2 if missing (class )replace class = "east" if missing (class )append using eastpartsave "west_eastdf" , replace use provdb, clear grmap using provcoord, id(ID) line (data(linestring2) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" "196 0 3" ) size(*1.2 *1.2 *2 *2 *2)) graphr(margin(vlarge)) polygon(data(west_eastdf) by (class ) fcolor("234 200 98%60" "127 210 255%60" )) gr_edit .plotregion1.AddTextBox added_text editor 23 80gr_edit .plotregion1.added_text[1].text = {}gr_edit .plotregion1.added_text[1].text.Arrpush 西部地区范围gr_edit .plotregion1.AddTextBox added_text editor 10 90gr_edit .plotregion1.added_text[2].text = {}gr_edit .plotregion1.added_text[2].text.Arrpush 东部地区范围
然后就可以使用 west_eastdf 进行东西分类了：
use df, clear geoinpoly 纬度 经度 using west_eastdf gen class = "东部" replace class = "西部" if _ID == 1drop _IDsave "东西部处理结果" , replace list in 1/10
使用秦岭-淮河线分南北 方法类似，不过需要注意秦岭-淮河线的数据方向是从右边开始到左边的，所以构建多边形的时候要从 rightpoint -> qh -> leftpoint -> bottomleft -> bottomright -> rightpoint。
再看一下图：
use provdb, clear grmap using provcoord, id(ID) line (data(linestring) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" ) size(*1.2 *1.2 *2 *2)) graphr(margin(vlarge)) polygon(data(cnbbox) fcolor("none" )) label (data(cnbbox) x(_X) y (_Y) l (label ) length (40))
可以选择 cnbbox 左右框线的中点作为延伸点：
use cnbbox, clear clear input _X _Y73.50236 28.70033 end save leftpoint, replace clear input _X _Y135.0957 28.70033 end save rightpoint, replace
构造南北方区域数据：
use cnbbox, clear keep if label == "topright" save topright, replace use cnbbox, clear keep if label == "bottomleft" save bottomleft, replace use emptydata, clear append using rightpointappend using qhcoordappend using leftpointappend using topleftappend using toprightappend using rightpointdrop _ID class label gen _ID = 1gen class = "north" save northpart, replace use emptydata, clear append using rightpointappend using qhcoordappend using leftpointappend using bottomleftappend using bottomrightappend using rightpointdrop _ID class label gen _ID = 2gen class = "south" save southpart, replace use northpart, clear append using southpartsave north_southdf, replace use provdb, clear grmap using provcoord, id(ID) line (data(linestring) by (classgroup) color("black" "0 85 170" "0 193 155" "196 0 3" ) size(*1.2 *1.2 *2 *2)) graphr(margin(vlarge)) polygon(data(north_southdf) by (class ) fcolor("234 200 98%60" "127 210 255%60" )) gr_edit .plotregion1.AddTextBox added_text editor 49 99gr_edit .plotregion1.added_text[1].text = {}gr_edit .plotregion1.added_text[1].text.Arrpush 北方区域范围gr_edit .plotregion1.AddTextBox added_text editor 14 80gr_edit .plotregion1.added_text[2].text = {}gr_edit .plotregion1.added_text[2].text.Arrpush 南方区域范围
然后就可以使用 north_southdf 进行南北方分类：
use df, clear geoinpoly 纬度 经度 using north_southdf gen class = "北方" replace class = "南方" if _ID == 2drop _IDsave "南北方处理结果" , replace
在另外的课程中我们也绘制过工企南北方的分布图：
感兴趣的小伙伴可以继续学习该课程掌握 Stata 中绘制精美地图的方法。
点击这里跳转到 RStata 短书平台获取附件：Stata：如何根据工企经纬度判断所处的南北方、东西部
×
×
×
评论