「Place-based policies, state-led industrialisation, and regional development: Evidence from China’s Great Western Development Programme」一文中使用到了西部大开发边界相关的距离数据：
今天我们一起来看下如何使用 R 语言提取西部大开发边界并计算各城市质心与该边界的距离。
西部大开发主要针对 12 个省份：新疆维吾尔自治区、甘肃省、内蒙古自治区、西藏自治区、青海省、宁夏回族自治区、陕西省、四川省、重庆市、云南省、贵州省和广西壮族自治区，以及在此之外的 3 个地级市：湘西土家族苗族自治州、延边朝鲜族自治州、恩施土家族苗族自治州。
首先我们读取中国市级行政区划：
library( tidyverse) library( sf) read_sf( "2019行政区划/市.shp" ) %>% select( - contains( "类型" ) ) %>% st_transform( 4326 ) -> city city
#> Simple feature collection with 371 features and 4 fields #> Geometry type: MULTIPOLYGON #> Dimension: XY #> Bounding box: xmin: 73.50114 ymin: 6.323421 xmax: 135.0885 ymax: 53.5609 #> Geodetic CRS: WGS 84 #> # A tibble: 371 × 5 #> 省代码 省 市代码 市 geometry #> * <dbl> <chr> <dbl> <chr> <MULTIPOLYGON [°]> #> 1 110000 北京市 110000 北京市 (((116.6753 41.0401, 116.6762 41.04006, 116.67… #> 2 120000 天津市 120000 天津市 (((117.4438 40.25101, 117.4561 40.24615, 117.4… #> 3 130000 河北省 130100 石家庄市 (((113.8242 38.75805, 113.8312 38.74815, 113.8… #> 4 130000 河北省 130200 唐山市 (((118.8539 39.10692, 118.8493 39.10679, 118.8… #> 5 130000 河北省 130300 秦皇岛市 (((119.1521 40.6128, 119.1517 40.60917, 119.15… #> 6 130000 河北省 130400 邯郸市 (((113.8711 37.01219, 113.8724 37.01182, 113.8… #> 7 130000 河北省 130500 邢台市 (((115.1259 37.79847, 115.1287 37.79843, 115.1… #> 8 130000 河北省 130600 保定市 (((115.4378 39.95016, 115.4435 39.9472, 115.44… #> 9 130000 河北省 130700 张家口市 (((114.8005 42.14749, 114.8045 42.14733, 114.8… #> 10 130000 河北省 130800 承德市 (((117.7998 42.6137, 117.8 42.61273, 117.7997 … #> # ℹ 361 more rows
然后我们筛选出西部大开发政策的地区：
city %>% filter( 省 %in% c ( "新疆维吾尔自治区" , "甘肃省" , "内蒙古自治区" , "西藏自治区" , "青海省" , "宁夏回族自治区" , "陕西省" , "四川省" , "重庆市" , "云南省" , "贵州省" , "广西壮族自治区" ) | 市 %in% c ( "湘西土家族苗族自治州" , "延边朝鲜族自治州" , "恩施土家族苗族自治州" ) ) -> city1 mapview:: mapview( city1)
然后我们把这部分区域融合起来：
city1 %>% st_make_valid( ) %>% nngeo:: st_remove_holes( ) %>% st_union( ) %>% st_sf( ) -> city1_union mapview:: mapview( city1_union)
提取该部分的边界保存成 geojson 文件：
city1_union %>% st_cast( "MULTILINESTRING" ) %>% st_simplify( dTolerance = 5000 ) %>% write_sf( "line1.geojson" , overwrite = T )
下面我们再借助一个在线应用提取边界：https://tidyfriday.cn/geojson/
删除 line1 的线条：
然后把右侧的内容保存到 polygon.geojson 文件里面。
这样我们就可以根据这个多边形提取西部大开发的边界了：
read_sf( "polygon.geojson" ) -> polygon st_union( city %>% filter( 市 == "延边朝鲜族自治州" ) %>% st_cast( "MULTILINESTRING" ) %>% pull( geometry) , city1_union %>% st_cast( "MULTILINESTRING" ) %>% st_intersection( polygon) ) %>% st_sf( ) -> lineall mapview:: mapview( lineall)
lineall %>% mutate( name = "西部大开发边界" ) %>% write_sf( "西部大开发边界.geojson" , overwrite = T )
最后我们再计算各城市质心到西部大开发边界的距离：
read_sf( "西部大开发边界.geojson" ) -> lineall read_sf( "2019行政区划/市.shp" ) %>% select( - contains( "类型" ) ) %>% st_transform( 4326 ) -> city city %>% st_make_valid( ) %>% st_centroid( ) %>% st_distance( lineall) -> distmat ( distmat/ 1000 ) %>% as.numeric ( ) %>% as_tibble( ) %>% bind_cols( city %>% st_drop_geometry( ) , . ) %>% rename( 各城市质心到西部大开发边界的距离 = value) %>% writexl:: write_xlsx( "各城市质心到西部大开发边界的距离_km.xlsx" )
这样就完成了这项工作。
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言提取西部大开发边界及计算各城市与该边界的距离
评论