извлекать координаты точек в простом фрейме данных объекта - PullRequest
Новая
       48

извлекать координаты точек в простом фрейме данных объекта

0 голосов
/ 23 октября 2018

У меня есть 2 линии, и мне нужно извлечь координаты, где эти 2 линии пересекаются.Я пытаюсь сделать это с пакетом sf в R. Я могу сделать пересечение, и я могу построить их, но я не могу извлечь координаты.Это сработало, когда на линии имеется только 1 точка пересечения, но когда их несколько (а в реальных данных их много), координаты точек пересечения сохраняются в одном поле на линию.

inters_pt выглядит следующим образом: геометрия класса A c (3, 3, 3.2, 4.2) B c (1.5, 1.5) A c (2, 4)

Я получаю сообщение:

Ошибка в st_coordinates.sfc (st_geometry (x)): не реализовано для объектов класса sfc_GEOMETRY

Я хотел бы извлечь координаты каждого пересечения и сохранить атрибутв столбце «класс».

Вывод должен выглядеть следующим образом: геометрия класса A (3, 3,2) A (3, 4,2) B (1,5, 1,5) A (2, 4)

Хотя изначально я думал, что это будет выглядеть так: класс геометрии A (3, 3) A (3.2, 4.2)

Я был не прав.

Поиск меня здесь получил: https://github.com/r-spatial/sf/issues/114

Но решения там либо отбрасывают точки (st_cast (mp, "POINT"), либо теряют атрибут 'class', либо выдают следующую ошибку:

Ошибка в vapply(lst, class, rep (NA_character_, 3)): значения должны иметь длину 3, но результат FUN (X [ 1 ]) - длина 2 

#First create my data, I create points, and convert them into SpatialLines
lines1.df <- data.frame(
  x = c(1,2,3,3,3,3,1,2,3),
  y = c(1,2,2,3,4,5,4,4,4),
  id = c(rep("A",6), rep("B",3))
) 

#with Kyle Walker's functions I convert the points to lines
#https://rpubs.com/walkerke/points_to_line 

library(sp)
library(maptools)

points_to_line <- function(data, long, lat, id_field = NULL, sort_field = NULL) {

  # Convert to SpatialPointsDataFrame
  coordinates(data) <- c(long, lat)

  # If there is a sort field...
  if (!is.null(sort_field)) {
    if (!is.null(id_field)) {
      data <- data[order(data[[id_field]], data[[sort_field]]), ]
    } else {
      data <- data[order(data[[sort_field]]), ]
    }
  }

  # If there is only one path...
  if (is.null(id_field)) {

    lines <- SpatialLines(list(Lines(list(Line(data)), "id")))

    return(lines)

    # Now, if we have multiple lines...
  } else if (!is.null(id_field)) {  

    # Split into a list by ID field
    paths <- sp::split(data, data[[id_field]])

    sp_lines <- SpatialLines(list(Lines(list(Line(paths[[1]])), "line1")))

    # I like for loops, what can I say...
    for (p in 2:length(paths)) {
      id <- paste0("line", as.character(p))
      l <- SpatialLines(list(Lines(list(Line(paths[[p]])), id)))
      sp_lines <- spRbind(sp_lines, l)
    }

    return(sp_lines)
  }
}


lines1 <- points_to_line(data = lines1.df, 
                        long = "x", 
                        lat = "y", 
                        id_field = "id")

proj4string(lines1) <- CRS("+proj=utm +zone=32 +datum=WGS84")

library(sf)

lines1_sf <- st_as_sf(lines1,
                         crs = "+proj=utm +zone=32 +datum=WGS84")

lines2.df <- data.frame(
  x = c(2,2,2.5,3.5,3.5,2.5,1,2),
  y = c(3,4.2,4.2,4.2,3.2,3.2,2,1),
  id = c(rep("A",6), rep("B",2))
) 

lines2 <- points_to_line(data = lines2.df, 
                            long = "x", 
                            lat = "y", 
                            id_field = "id")

proj4string(lines2) <- CRS("+proj=utm +zone=32 +datum=WGS84")

lines2_sf <- st_as_sf(lines2,
                       crs = "+proj=utm +zone=32 +datum=WGS84")
#add attributes
lines2_sf$class <- c("A", "B")

#plot both lines
library(ggplot2)
ggplot(lines1_sf) +
  geom_sf(aes(color="red")) +
  geom_sf(data=lines2_sf) +
  coord_sf(crs = "+proj=utm +zone=32 +datum=WGS84")  +
  theme_bw()

#find intersecting points
#intersect 
inters_pt <- st_intersection(lines2_sf, lines1_sf)

#plot shows the intersecting points
ggplot(lines1_sf) +
  geom_sf(aes(color="red")) +
  geom_sf(data=lines2_sf) +
  geom_sf(data=inters_pt) +
  coord_sf(crs = "+proj=utm +zone=32 +datum=WGS84")  +
  theme_bw()

#extract the coordinates from the geometry column
st_coordinates(inters_pt) #<-- gives error

st_cast(inters_pt, "POINT") #removes points

y = as(inters_pt, "Spatial") #loses the 'class' attribute
SpatialPoints(y)

st_cast(st_sfc(inters_pt), "POINT", group_or_split = FALSE) #error

1 Ответ

0 голосов
/ 23 октября 2018

Во-первых, стоит использовать sf для создания линий из точек, что уменьшает количество требуемого кода: 

library(sf)
library(dplyr)

#First create my data, I create points, and convert them into SpatialLines
lines1.df <- data.frame(
  x = c(1,2,3,3,3,3,1,2,3),
  y = c(1,2,2,3,4,5,4,4,4),
  id = c(rep("A",6), rep("B",3))
) 

lines2.df <- data.frame(
  x = c(2,2,2.5,3.5,3.5,2.5,1,2),
  y = c(3,4.2,4.2,4.2,3.2,3.2,2,1),
  id = c(rep("A",6), rep("B",2))
) 

points_to_line <- function(data, group = "id"){
  data <- data %>% 
    group_by_at(group) %>%
    summarise(do_union = FALSE) %>%
    st_cast("LINESTRING") %>%
    ungroup %>%
    select(-do_union)
}

crs <- "+proj=utm +zone=32 +datum=WGS84"
coords <- c("x", "y")
lines1_sf <- st_as_sf(lines1.df, coords = coords, crs = crs) %>%
  points_to_line()
lines2_sf <- st_as_sf(lines2.df, coords = coords, crs = crs) %>%
  points_to_line()

#add attributes
lines2_sf$class <- c("A", "B")

#find intersecting points
#intersect 
inters_pt <- st_intersection(lines1_sf, lines2_sf)

Здесь вам нужно сначала привести к MULTIPOINT, а затем к POINT 

#extract the coordinates from the geometry column
inters_pt <- inters_pt %>%
  st_cast("MULTIPOINT") %>%
  st_cast("POINT")

> st_coordinates(inters_pt) 
#    X   Y
#1 3.0 3.2
#2 3.0 4.2
#3 2.0 4.0
#4 1.5 1.5
...