如何使用 R 语言计算各省份质心到秦岭淮河线的纬度差距离

最近有个小伙伴想要计算一些点距离秦岭淮河线的纬度差距离（也就是沿着纬度线的距离）。今天我们就来一起学习下计算这个距离的方法。

这里我们以每个省的质心为例进行讲解。

首先加载所需的 R 包：

library(tidyverse)
library(sf)

读取秦岭淮河线和省份的矢量数据：

read_sf("秦岭-淮河线/秦岭-淮河线.shp") -> qh
read_sf("2020行政区划/省.shp") -> prov
read_sf("九段线/九段线.shp") -> jdx

计算各个省的质心：

prov %>%
  st_centroid() -> prov_centroid
prov_centroid

#> Simple feature collection with 34 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 85.19159 ymin: 17.86559 xmax: 127.7673 ymax: 47.86168
#> Geodetic CRS:  WGS 84
#> # A tibble: 34 x 4
#>    省             省代码 省类型                geometry
#>  * <chr>           <dbl> <chr>              <POINT [°]>
#>  1 安徽省         340000 省          (117.2261 31.8259)
#>  2 澳门特别行政区 820000 特别行政区 (113.5622 22.16291)
#>  3 北京市         110000 直辖市     (116.4126 40.18563)
#>  4 福建省         350000 省         (118.0047 26.05492)
#>  5 甘肃省         620000 省         (100.6046 37.88112)
#>  6 广东省         440000 省         (113.4132 23.30776)
#>  7 广西壮族自治区 450000 自治区     (108.7883 23.82111)
#>  8 贵州省         520000 省         (106.8746 26.81563)
#>  9 海南省         460000 省         (110.3052 17.86559)
#> 10 河北省         130000 省         (116.1432 39.54576)
#> # … with 24 more rows

注意并不是每个省份的质心都可以计算距离秦岭淮河线的纬度差距离（因为他们所处的经线并不和秦岭淮河线相交），可以计算的有这些：

# 秦岭淮河线的范围
st_bbox(qh) -> qhrange
qhrange

#>      xmin      ymin      xmax      ymax
#>  79.03223  31.71838 120.85553  36.36456

bind_cols(
  prov_centroid %>% st_drop_geometry(),
  prov_centroid %>% st_coordinates() %>%
    as_tibble() %>%
    set_names(c("经度", "纬度"))
) %>%
  dplyr::filter(between(经度, 79.0322, 120.85553))

#> # A tibble: 29 x 5
#>    省             省代码 省类型      经度  纬度
#>    <chr>           <dbl> <chr>      <dbl> <dbl>
#>  1 安徽省         340000 省          117.  31.8
#>  2 澳门特别行政区 820000 特别行政区  114.  22.2
#>  3 北京市         110000 直辖市      116.  40.2
#>  4 福建省         350000 省          118.  26.1
#>  5 甘肃省         620000 省          101.  37.9
#>  6 广东省         440000 省          113.  23.3
#>  7 广西壮族自治区 450000 自治区      109.  23.8
#>  8 贵州省         520000 省          107.  26.8
#>  9 海南省         460000 省          110.  17.9
#> 10 河北省         130000 省          116.  39.5
#> # … with 19 more rows

我们先以安徽省为例计算：

安徽省的质心坐标为：

prov_centroid %>%
  dplyr::filter(省 == "安徽省") %>%
  st_coordinates() -> ah

ah

#>          X       Y
#> 1 117.2261 31.8259

这个点往北往南的坐标分别是：

ah2 <- c(ah[1], 0)
ah3 <- c(ah[1], 90)

这两个点可以连成一条线：

c(ah2, ah3) %>%
  matrix(nrow = 2, byrow = T) %>%
  st_linestring() %>%
  st_sfc(crs = 4326) -> ahline

这条线和秦岭-淮河线的交点是：

# 计算交点
st_intersection(ahline, qh) -> ahpoint

我们可以看看这个计算的过程：

library(ggspatial)
# 绘图展示
mycrs <- "+proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
prov %>%
  st_transform(mycrs) %>%
  st_bbox() -> provbbox
ggplot(qh) +
  geom_sf() +
  geom_sf(data = ahline) +
  geom_sf(data = ahpoint, color = "green") +
  geom_sf(data = prov, size = 0.1,
          color = "black", fill = NA) +
  geom_sf(data = subset(prov_centroid, 省 == "安徽省"), color = "red") +
  geom_sf(data = jdx, size = 0.5, color = "black") +
  coord_sf(crs = mycrs,
           xlim = c(provbbox[1], provbbox[3]),
           ylim = c(provbbox[2], provbbox[4])) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_x = unit(0.75, "cm"),
    pad_y = unit(0.5, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
  ))

地图示意

这个交点距离安徽省的质心是：

# 距离
prov_centroid %>%
  dplyr::filter(省 == "安徽省") %>%
  st_distance(ahpoint) %>%
  .[,1] %>%
  as.numeric()

#> [1] 138964.9

那么我们就可以编写一个函数，这个函数的输入参数是省份的名字，输出的结果是纬度差距离：

# 编写函数
get_vertical_distance <- function(x) {
  print(x)
  prov_centroid %>%
    dplyr::filter(省 == x) %>%
    st_coordinates() -> ah

  ah2 <- c(ah[1], 0)
  ah3 <- c(ah[1], 90)

  c(ah2, ah3) %>%
    matrix(nrow = 2, byrow = T) %>%
    st_linestring() %>%
    st_sfc(crs = 4326) -> ahline

  # 计算交点
  st_intersection(ahline, qh) -> ahpoint

  # 距离
  prov_centroid %>%
    dplyr::filter(省 == x) %>%
    st_distance(ahpoint) %>%
    .[,1] %>%
    as.numeric()
}

# 例如计算北京的：
get_vertical_distance("北京市")
#> [1] 795812.4

然后我们就可以计算所的省份的了（去掉不能计算的）：

bind_cols(
  prov_centroid %>% st_drop_geometry(),
  prov_centroid %>% st_coordinates() %>%
    as_tibble() %>%
    set_names(c("经度", "纬度"))
) %>%
  dplyr::filter(between(经度, 79.0322, 120.85553)) %>%
  mutate(dist = map_dbl(省, get_vertical_distance)) -> provdist

provdist
#> # A tibble: 29 x 6
#>    省             省代码 省类型      经度  纬度     dist
#>    <chr>           <dbl> <chr>      <dbl> <dbl>    <dbl>
#>  1 安徽省         340000 省          117.  31.8  138965.
#>  2 澳门特别行政区 820000 特别行政区  114.  22.2 1177319.
#>  3 北京市         110000 直辖市      116.  40.2  795812.
#>  4 福建省         350000 省          118.  26.1  785467.
#>  5 甘肃省         620000 省          101.  37.9  585196.
#>  6 广东省         440000 省          113.  23.3 1049223.
#>  7 广西壮族自治区 450000 自治区      109.  23.8  951532.
#>  8 贵州省         520000 省          107.  26.8  645296.
#>  9 海南省         460000 省          110.  17.9 1624295.
#> 10 河北省         130000 省          116.  39.5  727145.
#> # … with 19 more rows

这里距离的单位是米。

所以计算思路就是，沿着该点所在的经度线选择两点，确保这两点的连线与秦岭淮河线相交，然后计算交点坐标。最后就可以计算该点距离交点的距离了！

点击这里跳转到 RStata 短书平台获取附件：如何使用 R 语言计算各省份质心到秦岭淮河线的纬度差距离