在学习本课程前需要预先学习使用 Stata 绘制省级地图的课程：「使用 Stata 绘制历年中国省级行政区划（小地图版本 + 长版）」
之前给大家分享过使用 Stata 绘制中国地图的教程，例如省级的：
最近有小伙伴想要绘制带邻国地区的地图，也就是类似这样的：
于是我就想是不是可以设计一份辅助数据，这样就可以直接在之前方法的基础上直接添加邻国了。
这里以带小地图的版本为例进行讲解。
使用 R 语言设计数据
不会 R 语言的小伙伴可以直接跳过这部分，设计的结果可以直接在 Stata 中使用的。
加载所需的 R 包：
library( tidyverse) library( sf)
之前提供的 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"
读取之前设计的 shp 文件：
read_sf( "chinaprov2021mini/chinaprov2021mini.shp" ) %>% st_transform( mycrs) -> prov prov
这个数据的经纬度范围是：
然后我们使用这个范围（最好再扩大点）从世界地图上截取部分：
st_bbox( c ( xmin = - 2725586 , xmax = 2962768 , ymax = 6000000 , ymin = 1808655 ) , crs = st_crs( mycrs) ) %>% st_as_sfc( ) -> provbbox read_sf( "worldmap0/worldmap0.shp" ) %>% st_transform( mycrs) -> world world %>% st_intersection( provbbox) -> asia asia %>% st_simplify( dTolerance = 2000 ) -> asia_sim plot( asia_sim[ 1 ] )
然后再把中国区域的数据去除：
prov %>% filter( ! is.na ( 省) ) -> prov2 asia %>% st_difference( st_union( prov2) ) -> asia2 asia2 %>% write_rds( "asia2.rds" ) read_rds( "asia2.rds" ) -> asia2 asia2 %>% st_simplify( dTolerance = 2000 ) -> asia2_sim plot( asia2_sim[ 1 ] )
添加中文国家名称变量：
asia2 %>% mutate( country_cn = c ( "阿富汗" , "孟加拉国" , "不丹" , "中国" , "印度" , "日本" , "哈萨克斯坦" , "吉尔吉斯斯坦" , "老挝" , "蒙古" , "缅甸" , "尼泊尔" , "朝鲜" , "巴基斯坦" , "菲律宾" , "俄罗斯" , "锡亚琛冰川" , "韩国" , "塔吉克斯坦" , "泰国" , "乌兹别克斯坦" , "越南" ) ) %>% select( - continent) -> asia2
由于带小地图版本的数据里面还有单独的小地图部分，所以再提取小地图范围的：
small_bbox <- st_bbox( c ( xmin = 120000 , xmax = 1766004.1 , ymax = 2557786.0 , ymin = 320000 ) , crs = st_crs( mycrs) ) %>% st_as_sfc( ) world %>% st_intersection( small_bbox) %>% mutate( geometry = geometry * 0.5 + c ( 2100000 , 1665139 ) ) %>% sf:: st_set_crs( mycrs) %>% mutate( country_cn = c ( "文莱" , "柬埔寨" , "中国" , "印度尼西亚" , "老挝" , "马来西亚" , "菲律宾" , "越南" ) ) -> asia3 asia3 %>% st_difference( st_union( prov2) ) -> asia3a
然后合并两部分的数据：
bind_rows( asia2, asia3a) %>% select( - continent) %>% group_by( country, iso3, country_cn) %>% summarise( ) %>% ungroup( ) -> asia4
保存为 shp 格式的数据：
fs:: dir_delete( "china_neighboring" ) dir.create( "china_neighboring" ) asia4 %>% filter( country_cn != "中国" ) %>% st_make_valid( ) %>% st_collection_extract( "POLYGON" ) %>% st_write( "china_neighboring/china_neighboring.shp" , layer_options = "ENCODING=UTF-8" , delete_layer = TRUE , layer = "MULTIPOLYGON" )
所以大家对上面的代码不能理解的话也没关系，之后直接使用这个处理好的 shp 数据即可。
在附件中 main2.R 文件中也有长版地图数据的生成代码，更简单。
在 Stata 绘图中使用 在 Stata 中首先我们需要把 shp 文件转换成 dta 文件：
local name = "china_neighboring" shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) gencentroids(centroids) replace
之前绘制的不带邻国的地图的代码这样的：
use chinaprov2021mini_db.dta, clear encode 省, gen (prov)codebook provgrmap prov using chinaprov2021mini_coord.dta, id(ID) osize(vvthin ...) ocolor(white ...) clmethod(unique) fcolor("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" "255 209 71" "153 0 51" "153 26 0" "153 102 0" "128 153 0" "51 153 0" "0 153 26" "0 153 102" "0 128 153" "0 51 153" "26 0 153" "102 0 153" "153 0 128" "214 0 71" "255 20 99" "0 214 143" "20 255 177" ) leg(off) graphr(margin(medium)) line (data(chinaprov2021mini_line_coord.dta) by (group) size(vvthin *1 *0.5 *1.2 *0.5 *0.5 *1.2) pattern(solid ...) color(white black "0 85 170" "24 188 156" black black "227 26 28" )) polygon(data(polygon) fcolor(black) osize(vvthin)) label (data(chinaprov2021mini_label) x(X) y (Y) label (cname) length (20) size(*0.8)) ti("使用 Stata 绘制 2021 年中国省级行政区划" ) subti("绘制：微信公众号 RStata" ) caption("版本：使用 Stata 绘制中国省级地图数据包2021" , size(*0.8)) gr export pic1c.png, replace width(4800)
由于 china_neighboring 主要是多边形数据，所以我们只要把它和 polygon.dta 数据合并即可：
use polygon, clear use china_neighboring_coord.dta, clear gen class = "邻国" gen value = 2append using polygonencode class , gen (class_group)codebook class_groupreplace _X = _X - 1000000 if index (class , "比例尺" )replace _X = _X + 3620000 if index (class , "指北针" )replace _Y = _Y + 3400000 if index (class , "指北针" )save polygon_with_nhb, replace
为了给这些邻国添加上文本标签，我们需要把标签数据也和 chinaprov2021mini_label.dta 合并起来：
use chinaprov2021mini_label.dta, clear use china_neighboring_db, clear keep x_ y_ country country_cnren country_cn cnameren country enameren x_ Xren y_ Yappend using chinaprov2021mini_label.dtareplace X = X - 1000000 if index (cname, "1000km" )replace X = X + 3620000 if index (cname, "N" )replace Y = Y + 3400000 if index (cname, "N" )replace X = X - 100000 if index (cname, "菲律宾" )save label_with_nhb, replace
比例尺向左平移、指北针移动到右上方（上面的代码中也有平移的代码）：
use chinaprov2021mini_line_coord.dta, clear replace _X = _X - 1000000 if index (class , "比例尺" )replace _X = _X + 3620000 if index (class , "指北针" )replace _Y = _Y + 3400000 if index (class , "指北针" )save line_with_nhb, replace
然后就可以绘图了：
use chinaprov2021mini_db.dta, clear encode 省, gen (prov)codebook provgrmap prov using chinaprov2021mini_coord.dta, id(ID) osize(vvthin ...) ocolor(white ...) clmethod(unique) fcolor("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" "255 209 71" "153 0 51" "153 26 0" "153 102 0" "128 153 0" "51 153 0" "0 153 26" "0 153 102" "0 128 153" "0 51 153" "26 0 153" "102 0 153" "153 0 128" "214 0 71" "255 20 99" "0 214 143" "20 255 177" ) leg(off) graphr(margin(medium)) line (data(line_with_nhb) by (group) size(vvthin *1.2 *0.5 *1.2 *0.5 *0.5 *1.2) pattern(solid ...) color(white "162 154 196" "0 85 170" "24 188 156" black black "227 26 28" )) 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("使用 Stata 绘制 2021 年中国省级行政区划" ) subti("绘制：微信公众号 RStata" ) caption("版本：使用 Stata 绘制中国省级地图数据包2021" , size(*0.8)) plotr(fcolor("187 209 235" ) margin(-0.5 -0.5 -0.5 -0.5)) gr export pic2c.png, replace width(4800)
在附件中 main2.do 文件中也有绘制长版地图的代码，方法类似，这里就不再赘述了。
当然因为时间关系，本文仅仅演示了 2021 年省级地图的绘制，实际上其他年份的、市级和县级的也可以使用类似的方法绘图。
点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制带邻国的中国地图（含使用 R 语言操作矢量数据的内容）
×
评论