使用 Stata 绘制中国乡镇地图及根据经纬度判断所处的乡镇

在之前的课程中，我给大家分享过使用 Stata 绘制省市区县地图的方法：

使用 Stata 绘制历年中国省级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/e194d75fe1674f7c8921daec1ece7f7d
使用 Stata 绘制历年中国市级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/0e9d78633ae64fefa684d6aad89633bc
使用 Stata 绘制历年中国县级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/99c2e97b88ea401bbdb8da4511341cc8

另外也给大家分享过使用 Stata 进行地理编码以及根据经纬度判断所处省市区县的方法：

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

乡镇地图的绘制

今天我们再来补充讲解下如何绘制乡镇地图以及如何根据经纬度判断所处的乡镇。

附件中的 town4326 文件夹是 WGS84 坐标系下的乡镇地图数据；town_mini 文件夹是经过编辑后的小地图版本的乡镇地图数据；town_long 则是长版本的。

town_mini 和 town_long 都是“+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 文件转换成 dta 格式：

*- 把 shp 文件转换成 dta 文件
local name = "town_mini"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) replace
local name = "town_long"
shp2dta using `name'/`name', database(`name'_db) coordinates(`name'_coord) genid(ID) replace

国界线、九段线、海岸线之类的数据都可以使用之前课程里面的。

这里使用 2019 年各乡镇的平均人口密度数据演示乡镇地图的绘制：

*- 各乡镇人口密度
use town_mini_db.dta, clear
ren townname 乡镇
ren towncode 乡镇代码
merge 1:1 乡镇 乡镇代码 using 2019年各乡镇平均人口密度.dta
drop _m
_pctile mean, n(8)
ret list
replace mean = -1 if mi(mean)
grmap mean using town_mini_coord.dta, id(ID) ///
 ocolor("black" ...) osize(vvthin ...) ///
 clmethod(custom) ///
 clbreaks(-1 0 68 125 197 312 458 663 1360 50000) ///
 fcolor("white" "252 251 253" "239 237 245" ///
  "218 218 235" "188 189 220" "158 154 200" ///
  "128 125 186" "106 81 163" "74 20 134") ///
 line(data(chinaprov2021mini_line_coord3.dta) ///
     /// 胡焕庸线（7）秦岭淮河线（4）
     select(keep if inlist(group, 1, 2, 3, 5, 6)) ///
     by(group) size(vvthin *1 *0.5 *0.5 *0.5) ///
     pattern(solid ...) ///
     color(white /// 省界颜色
         black /// 国界线颜色
         "0 85 170" /// 海岸线颜色
         black /// 小地图框格颜色
         black /// 比例尺和指北针颜色
         )) ///
 polygon(data(polygon3) fcolor(black) ///
  osize(vvthin)) ///
 label(data(chinacity2021mini_label3) x(X) y(Y) ///
  label(cname) length(20) size(*0.8)) ///
 leg(order(2 "无数据" 3 "<=68" 4 "68~125" 5 "125~197" ///
  6 "197~312" 7 "312~458" 8 "458~663" 9 "663~1360" ///
  10 ">1360") ti("人口密度(人/m{sup:2})", size(*0.5)) ///
  col(2)) ///
    ti("2019 年各乡镇平均人口密度", color(black)) ///
    subti("绘制：微信公众号 RStata") ///
    graphr(margin(medium)) ///
    caption("数据来源：中国科学院资源环境科学与数据中心", size(*0.8))

gr export 2019年各乡镇平均人口密度.png, width(4800) replace

绘图细节难以一一介绍，所以这里建议结合视频讲解学习。

长版地图的绘制方法类似：

*- 长版本
use town_long_db.dta, clear
ren townname 乡镇
ren towncode 乡镇代码
merge 1:1 乡镇 乡镇代码 using 2019年各乡镇平均人口密度.dta
drop _m
_pctile mean, n(8)
ret list
replace mean = -1 if mi(mean)
grmap mean using town_long_coord.dta, id(ID) ///
 ocolor("black" ...) osize(vvthin ...) ///
 clmethod(custom) ///
 clbreaks(-1 0 68 125 197 312 458 663 1360 50000) ///
 fcolor("white" "252 251 253" "239 237 245" ///
  "218 218 235" "188 189 220" "158 154 200" ///
  "128 125 186" "106 81 163" "74 20 134") ///
 line(data(chinaprov2021long_line_coord3.dta) ///
  by(group) size(vvthin *1 *0.5 *0.5) ///
  pattern(solid ...) ///
  select(drop if inlist(group, 4, 7)) ///
  color(white /// 省界颜色
     black /// 国界线颜色
     "0 85 170" /// 海岸线颜色
     black /// 比例尺和指北针颜色
     )) ///
 polygon(data(longpolygon3) fcolor(black) ///
  osize(vvthin)) ///
 label(data(chinacity2021long_label3) x(X) y(Y) ///
  label(cname) length(20) size(*0.8)) ///
 leg(order(2 "无数据" 3 "<=68" 4 "68~125" 5 "125~197" ///
  6 "197~312" 7 "312~458" 8 "458~663" 9 "663~1360" ///
  10 ">1360") ti("人口密度(人/m{sup:2})", size(*0.5)) ///
  col(1)) ///
    ti("2019 年各乡镇平均人口密度", color(black)) ///
    subti("绘制：微信公众号 RStata") ///
    graphr(margin(medium)) ///
    caption("数据来源：中国科学院资源环境科学与数据中心", size(*0.8))

gr export 2019年各乡镇平均人口密度_long.png, width(4800) replace

根据经纬度判断所处的乡镇

使用 geoinpoly 命令就可以进行这一操作了，首先把 town4326 转换成 dta 文件：

shp2dta using town4326/town4326.shp, database(town4326_db) coordinates(town4326_coord) genid(ID) replace

工企样本.dta 是大概 1 万条工企数据的样本，含有经纬度变量，使用下面的代码就可以判断每个坐标点所处的乡镇了：

use 工企样本.dta, clear
geoinpoly 纬度 经度 using town4326_coord.dta
ren _ID ID
merge m:1 ID using town4326_db
keep if _m == 3
drop _m
foreach i of varlist _all {
 cap format `i' %10s
}
drop ID
ren Name 乡镇
ren code 乡镇代码
*- 这样我们就完成了判断所处乡镇的操作
save finalresult, replace

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制中国乡镇地图及根据经纬度判断所处的乡镇

乡镇地图的绘制

根据经纬度判断所处的乡镇

相关文章

评论