名师讲堂｜工具变量的计算：“八纵八横”光缆干线网络及其与网点城市各省市区县质心的距离

最近给大家分享了「经济学季刊」 2023 年 1 月 里面的「数字经济、家庭分工与性别平等」一文中的工具变量使用案例：

这个文献里面提到的两个工具变量都可以非常容易地使用 R 语言进行计算。

加载相关 R 包与读取矢量数据

首先我们加载所需的 R 包：

library(tidyverse)
library(sf)

其中 tidyverse 用于数据处理，sf 用于地理计算。

以 2019 年的省市区县行政区划为例：

read_sf("2019行政区划/省.shp") %>%
  st_make_valid() %>%
  st_transform(4326) -> prov

read_sf("2019行政区划/市.shp") %>%
  st_make_valid() %>%
  st_transform(4326) -> city

read_sf("2019行政区划/县.shp") %>%
  st_make_valid() %>%
  rename(县代码 = PAC, 县 = NAME) %>%
  st_transform(4326) -> county

# 计算质心
prov %>%
  st_centroid() -> prov_centroid
city %>%
  st_centroid() -> city_centroid
county %>%
  st_centroid() -> county_centroid

准备“八纵八横”光缆干线网络数据

关于“八纵八横”光缆干线，我找到了这样的一幅图：

然后我手工抄录了里面所有的连接线：line.csv 文件，下面我们将基于这个文件来构建矢量数据：

read_csv("line.csv") %>%
  mutate(ID = row_number()) -> linedf

linedf

#> # A tibble: 109 × 2
#>    line               ID
#>    <chr>           <int>
#>  1 齐齐哈尔-哈尔滨     1
#>  2 哈尔滨-牡丹江       2
#>  3 白城-齐齐哈尔       3
#>  4 白城-长春           4
#>  5 长春-延吉           5
#>  6 延吉-牡丹江         6
#>  7 白城-阜新           7
#>  8 阜新-沈阳           8
#>  9 沈阳-长春           9
#> 10 沈阳-丹东          10
#> # ℹ 99 more rows

linedf %>%
  separate_rows(line, sep = "-") %>%
  select(ID, everything()) %>%
  mutate(city2 = substr(line, 1, 2)) %>%
  left_join(
    bind_rows(
      city_centroid %>%
        filter(!str_detect(市, "张家界市")) %>%
        mutate(city2 = str_sub(市, 1, 2)),
      county_centroid %>%
        filter(县 %in% c("延吉市", "伊宁市", "格尔木市")) %>%
        mutate(city2 = str_sub(县, 1, 2)) %>%
        select(city2, 市 = 县)
    )
  ) %>%
  select(ID, city = 市, geometry) -> linepointdf

# 这里还可以编写一个函数用于 point 合成 linestring
st_point_to_linestring <- function(..., crs) {
  c(...) %>%
    st_coordinates() %>%
    st_linestring() %>%
    st_sfc(crs = crs)
}

linepointdf %>%
  group_by(ID) %>%
  summarise(geometry = st_point_to_linestring(geometry, crs = 4326)) %>%
  st_sf() -> linedfsf

mapview::mapview(linedfsf)

插图

然后可以把这个矢量数据保存成 shp 文件：

linedfsf %>%
  left_join(linedf) %>%
  select(ID, line) -> linedfsf

# 保存成 shp
dir.create("88network")
linedfsf %>%
  st_make_valid() %>%
  st_write("88network/88network.shp", layer_options = "ENCODING=UTF-8", delete_layer = TRUE, layer = "LINESTRING")
#> Deleting layer `LINESTRING' failed
#> Writing layer `LINESTRING' to data source
#>   `88network/88network.shp' using driver `ESRI Shapefile'
#> options:        ENCODING=UTF-8
#> Writing 109 features with 2 fields and geometry type Line String.

网络节点城市是有这些：

linepointdf %>%
  distinct(city, .keep_all = T) %>%
  st_sf() %>%
  select(-ID) -> linepointsf

使用 ggplot2 绘图展示

绘制中国地图通常使用这样的坐标参考系：

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"

这里我想绘制一个带邻国地区的中国地图，范围是这个：

# 绘图区域
st_bbox(c(xmin = -2725586,
          xmax = 2982768,
          ymax = 6000000,
          ymin = 1800655),
        crs = st_crs(mycrs)) %>%
  st_as_sfc() -> plotbbox

读取我编辑设计好的 shp 文件：

类似的文件都可以从平台上的课程：「使用 Stata 绘制中国市级地图」中获取。

read_sf("chinacity2019mini/chinacity2019mini.shp") %>%
  filter(!is.na(省代码)) -> citymap

read_sf("chinacity2019mini/chinacity2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> citylinemap

read_sf("china_neighboring/china_neighboring.shp") -> china_neighboring

然后就可以绘制地图了：

library(ggspatial)
ggplot(citymap) +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(fill = "white") +
  geom_sf(data = china_neighboring, fill = "#EDEDED") +
  stat_sf_coordinates(data = china_neighboring,
                      geom = "text", color = "gray",
                      aes(label = country_cn), family = cnfont,
                      size = 2) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  geom_sf(data = linedfsf, color = "#00c19b", linewidth = 1) +
  geom_sf(data = linepointsf, color = "#c40003", size = 2) +
  scale_color_manual(
    values = c("九段线" = "black",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.8,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "“八纵八横”光缆骨干网络与骨干城市分布",
       caption = "绘制：微信公众号 RStata") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p

ggsave("pic1.png", width = 10, height = 8.5, device = png)

计算各省市区县质心到“八纵八横”光缆骨干城市的距离

首先我们把所有的点集合成一个 MULTIPOINT 对象：

linepointsf$geometry %>%
  st_coordinates() %>%
  st_multipoint() %>%
  st_sfc(crs = 4326) -> multipoint

st_distance() 返回的结果就是最小距离（大圆距离）：

prov_centroid %>%
  mutate(与光缆骨干城市的最小距离_km = prov_centroid %>%
           st_distance(multipoint) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与光缆骨干城市的最小距离_km = if_else(
    与光缆骨干城市的最小距离_km < 0.001, 0,
    与光缆骨干城市的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各省份质心与光缆骨干城市的最小距离(单位km).xlsx")

city_centroid %>%
  mutate(与光缆骨干城市的最小距离_km = city_centroid %>%
           st_distance(multipoint) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与光缆骨干城市的最小距离_km = if_else(
    与光缆骨干城市的最小距离_km < 0.001, 0,
    与光缆骨干城市的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各城市质心与光缆骨干城市的最小距离(单位km).xlsx")

county_centroid %>%
  mutate(与光缆骨干城市的最小距离_km = county_centroid %>%
           st_distance(multipoint) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与光缆骨干城市的最小距离_km = if_else(
    与光缆骨干城市的最小距离_km < 0.001, 0,
    与光缆骨干城市的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各区县质心与光缆骨干城市的最小距离(单位km).xlsx")

计算各省市区县质心到杭州的距离

这个就非常容易了：

city_centroid %>%
  filter(市 == "杭州市") -> hz

prov_centroid %>%
  mutate(与杭州的最小距离_km = prov_centroid %>%
           st_distance(hz) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与杭州的最小距离_km = if_else(
    与杭州的最小距离_km < 0.001, 0,
    与杭州的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各省份质心与杭州的最小距离(单位km).xlsx")

city_centroid %>%
  mutate(与杭州的最小距离_km = city_centroid %>%
           st_distance(hz) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与杭州的最小距离_km = if_else(
    与杭州的最小距离_km < 0.001, 0,
    与杭州的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各城市质心与杭州的最小距离(单位km).xlsx")

county_centroid %>%
  mutate(与杭州的最小距离_km = county_centroid %>%
           st_distance(hz) %>%
           .[,1] %>%
           units::set_units("km") %>%
           as.numeric()) %>%
  mutate(与杭州的最小距离_km = if_else(
    与杭州的最小距离_km < 0.001, 0,
    与杭州的最小距离_km
  )) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  filter(省 != "中朝共有") %>%
  writexl::write_xlsx("各区县质心与杭州的最小距离(单位km).xlsx")

绘制各城市到杭州的连接线网络

有时候我们会需要绘制某个城市到其他城市的辐射网络，这里以杭州为例进行演示：

# 各城市与杭州的连接线
city_centroid %>%
  filter(市 != "杭州市") %>%
  select(city = 市) -> cityfrom

city_centroid %>%
  filter(市 == "杭州市") %>%
  select(city = 市) -> cityto

cityfrom %>%
  mutate(ID = row_number()) -> cityfrom

cityto %>%
  crossing(ID = 1:nrow(cityfrom)) %>%
  bind_rows(cityfrom) %>%
  group_by(ID) %>%
  summarise(geometry = st_point_to_linestring(geometry, crs = 4326)) %>%
  st_sf() -> linedfsf2

# 绘图展示
ggplot(citymap) +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(fill = "white") +
  geom_sf(data = china_neighboring, fill = "#EDEDED") +
  stat_sf_coordinates(data = china_neighboring,
                      geom = "text", color = "gray",
                      aes(label = country_cn), family = cnfont,
                      size = 2) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  geom_sf(data = linedfsf2, color = "#00c19b", linewidth = 0.5) +
  geom_sf(data = city_centroid, color = "#c40003", size = 1) +
  scale_color_manual(
    values = c("九段线" = "black",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.8,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "所有城市到杭州的连接线",
       caption = "绘制：微信公众号 RStata") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p

ggsave("pic2.png", width = 10, height = 8.5, device = png)

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜工具变量的计算：“八纵八横”光缆干线网络及其与网点城市各省市区县质心的距离