使用Shapefile測量國家與赤道的距離

我想測量國家相對于赤道的距離，我想到了使用帶有 R 的 shape-file 來做到這一點。我不知道這有多可行，或者是否有可能，但我知道知道可以測量距離。

接下來測量國家和赤道線之間的距離。理想情況下，計算將產生一個包含國家和距離的面板資料。

這可能嗎？

uj5u.com熱心網友回復：

讓我們使用包獲取世界地圖 shapefile rnaturalearth：

library(sf)
library(rnaturalearth)
library(tidyverse)

map_world <- ne_countries(returnclass = 'sf')

為了簡化示例，讓我們選擇幾個不同的國家：

countries <- c("United States of America",
               "United Kingdom",
               "Australia",
               "Argentina",
               "Chad")

現在我們過濾掉這些國家并找到它們的質心。然后我們生成一個額外的 sfc 列，它是質心在赤道上的投影。最后，我們測量質心和赤道點之間的距離：

country_points <- map_world %>% 
  filter(admin %in% countries) %>% 
  st_centroid() %>%
  bind_cols(st_coordinates(.)) %>%
  select(admin, X, Y) %>%
  rowwise() %>%
  mutate(equator = st_sfc(lapply(X, function(x) st_point(c(x, 0))),
                          crs = st_crs(map_world))) %>%
  mutate(distance_equator = st_distance(geometry, equator)[,1]) %>%
  select(-equator)

這導致所有國家的質心到赤??道的距離測量：

country_points
#> Simple feature collection with 5 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -103.5729 ymin: -34.73683 xmax: 134.3143 ymax: 53.8083
#> CRS:            proj=longlat  datum=WGS84  no_defs  ellps=WGS84  towgs84=0,0,0
#> # A tibble: 5 x 5
#> # Rowwise: 
#>   admin                          X     Y              geometry distance_equator
#>   <chr>                      <dbl> <dbl>           <POINT [°]>              [m]
#> 1 Argentina                 -64.8  -34.7 (-64.75442 -34.73683)         3862566.
#> 2 Australia                 134.   -25.8   (134.3143 -25.7631)         2864730.
#> 3 United Kingdom             -2.76  53.8   (-2.759269 53.8083)         5983219.
#> 4 Chad                       18.6   15.3   (18.57333 15.27628)         1698647.
#> 5 United States of America -104.    44.8  (-103.5729 44.75598)         4976646.

此外，我們可以像這樣繪制結果：

ggplot(data = map_world)  
  geom_sf()  
  geom_hline(yintercept = 0, color = 'red', size = 1)  
  geom_segment(data = country_points, size = 1,
               aes(X, Y, xend = X, yend = 0), color = 'red')

^{使用reprex v2.0.2創建于 2022-10-27}

uj5u.com熱心網友回復：

找到從點到線的最短距離似乎不像我期望的那樣作業st_nearest_point，但從點到多邊形卻可以。艾倫卡梅隆的回答很好地解決了這個問題，因為到赤道的最短距離將具有與緯度設定為零的質心相同的經度坐標。

或者，赤道可以用多邊形表示，st_nearest_point可以直接使用。

library(tidyverse)
library(sf)
library(rnaturalearth)

# Use rnaturalearth for shapefile, keep only the name & geometry
world <- ne_countries(returnclass = 'sf') %>%
  select(name)
# Get country centroids. This will need cleaning!! Several countries aren't 
#  regular & continuous.
world_centroids <- st_centroid(st_make_valid(world))

# df for equator, polygons require at least 4 points
df <- data.frame(
  id = c(1,1,1,1)
  , x = c(179,0,-179,179)
  , y = c(0,0,0,0)
)
# sfheaders to make a polygon easily
equator <- sfheaders::sf_polygon( df, x = "x", y = "y", linestring_id = "id", keep = TRUE ) %>%
  st_set_crs(4326)

lines_to_equator <- st_nearest_points(world_centroids, equator)

ggplot()  
  geom_sf(data = world)  
  geom_sf(data = equator, color = 'blue')  
  geom_sf(data = world_centroids, color = 'green')  
  geom_sf(data = lines_to_equator, color = 'red')

使用 dplyr 將距離列添加到 world_centroids tibble，但st_distance 不適用于 equator 以南的質心。一種解決方法是st_length在st_nearest_points：

world_centroids %>% 
   mutate(dist_to_eq = st_distance(., equator, by_element = F))  %>%
   mutate(better_dist_to_eq = st_length(st_nearest_points(geometry, equator))) %>%
   head(9)
Simple feature collection with 9 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -64.75442 ymin: -84.96858 xmax: 134.3143 ymax: 41.1342
CRS:            proj=longlat  datum=WGS84  no_defs  ellps=WGS84  towgs84=0,0,0
                    name                    geometry  dist_to_eq better_dist_to_eq
            Afghanistan   POINT (66.00365 33.84035) 3762881 [m]       3762881 [m]
                 Angola   POINT (17.46526 -12.2259)       0 [m]       1359460 [m]
                Albania     POINT (20.0341 41.1342) 4573922 [m]       4573922 [m]
   United Arab Emirates    POINT (54.20068 23.8715) 2654394 [m]       2654394 [m]
              Argentina POINT (-64.75442 -34.73683)       0 [m]       3862566 [m]
                Armenia   POINT (45.00983 40.21435) 4471639 [m]       4471639 [m]
             Antarctica   POINT (82.5071 -84.96858)       0 [m]       9448089 [m]
 Fr. S. Antarctic Lands   POINT (69.53103 -49.3073)       0 [m]       5482730 [m]
              Australia   POINT (134.3143 -25.7631)       0 [m]       2864730 [m]
 

^{由reprex 包（v2.0.1）于 2022 年 10 月 27 日創建}

標籤：r科幻

上一篇：堆疊直方圖的X軸編輯

下一篇：獲取條件下列中每個向量的連續TRUE（或FALSE）值的數量