在R中使用SF-向大点数据集添加几何的最佳方法是什么?

问题描述 投票:0回答:1

免责声明:我只是开始使用sf,所以(希望!)这里可能缺少明显的东西。

我有AusGeoid2020数据,该数据由15,454,800点和一些属性组成,这些属性可在椭球高度(即GPS高度)和AHD之间进行转换。

尽管文件很大(914Mb),它很容易读入:

library(plyr)
library(magrittr)
library(dplyr)
library(readr)
library(sf)

AusGeoid2020 <- read_fwf(
  file = "AUSGeoid2020_20170908_win.dat",
  col_positions = fwf_widths(
    widths = c(3L,9L,2L,2L,3L,7L,2L,3L,3L,7L,10L,10L),
    col_names = c(
      "ID",
      "ellipsoid to AHD separation (m)",
      "Latitude (hem)",
      "Latitude (deg)",
      "Latitude (min)",
      "Latitude (sec)",
      "Longitude (hem)",
      "Longitude (deg)",
      "Longitude (min)",
      "Longitude (sec)",
      "deflection of the vertical (seconds, xi)",
      "deflection of the vertical (seconds, eta)"
    )
  ),
  col_types = cols(
    ID = col_character(),
    `ellipsoid to AHD separation (m)` = col_double(),
    `Latitude (hem)` = col_character(),
    `Latitude (deg)` = col_double(),
    `Latitude (min)` = col_double(),
    `Latitude (sec)` = col_double(),
    `Longitude (hem)` = col_character(),
    `Longitude (deg)` = col_double(),
    `Longitude (min)` = col_double(),
    `Longitude (sec)` = col_double(),
    `deflection of the vertical (seconds, xi)` = col_double(),
    `deflection of the vertical (seconds, eta)` = col_double()
  ),
  skip = 1L
)

AusGeoid2020 <- AusGeoid2020 %>% 
  mutate(
    Latitude = `Latitude (deg)` + (`Latitude (min)`/60) + (`Latitude (sec)`/3600),
    Latitude = case_when(
      `Latitude (hem)` == "S" ~ -1 * Latitude,
      TRUE ~ Latitude
    ),
    Longitude = `Longitude (deg)` + (`Longitude (min)`/60) + (`Longitude (sec)`/3600),
    Longitude = case_when(
      `Longitude (hem)` == "W" ~ -1 * Longitude,
      TRUE ~ Longitude
    )
  ) %>% 
  select(
    ID,
    `ellipsoid to AHD separation (m)`,
    Latitude,
    Longitude,
    `deflection of the vertical (seconds, xi)`,
    `deflection of the vertical (seconds, eta)`
  )

我的问题是:向此大型数据框添加几何的最佳方法是什么?我相信我想要的函数是st_point()而不是矢量化的,所以我不得不使用{plyr}中的alply()创建几何列,但这是[[very资源密集型的,这使我认为必须更好。

st_geometry(AusGeoid2020) <- st_sfc( alply(AusGeoid2020, 1, function(row) { st_point(x = c(row$Longitude, row$Latitude), dim = "XY") }), crs = 7844L )
这需要很长时间。任何建议表示赞赏!
r geospatial spatial sf sp
1个回答
0
投票
我们可以如下使用st_as_sf。默认设置将删除带有坐标信息的列(在这种情况下为LongitudeLatitude)。如果要保留这些列,请设置remove = FALSE

AusGeoid2020_sf <- AusGeoid2020 %>% st_as_sf(coords = c("Longitude", "Latitude"), crs = 7844L, remove = FALSE)

© www.soinside.com 2019 - 2024. All rights reserved.