2013-12-11 5 views
7

У меня есть 2 векторов:нахождения точки пересечения в R

set.seed(1) 
x1 = rnorm(100,0,1) 
x2 = rnorm(100,1,1) 

Я хочу, чтобы построить их как линии, а затем найти точки пересечения линий, а также, если существует несколько точек пересечения, то я хочу найдите каждый из них.

enter image description here

я наткнулся на подобный вопрос, и пытались решить эту проблему с помощью spatstat, но я не был в состоянии преобразовать мой комбинированный кадр данных, содержащий оба значения вектора к psp object.

+0

Вы хотите сказать, вы хотите, чтобы найти все пересечения линий в 'участка (x1, x2, тип =«L»)'? –

+0

Или вы имеете в виду переходы 'plot (seq_along (x1), x1, type = 'l')' и 'lines (seq_along (x2), x2, type = 'l', col =" red ")' –

+0

Я хочу, чтобы координаты, где бы ни была интуиция, я привел приведенные выше векторы как примеры игрушек, но моя фактическая серия является нелинейной, уравнение которой не указано. – SBS

ответ

11

Если вы буквально просто имеете два случайных вектора чисел, вы можете использовать довольно простой метод, чтобы получить пересечение обоих. Просто найдите все точки, где x1 выше x2, а затем под ним на следующей точке или наоборот. Это точки пересечения. Затем просто используйте соответствующие наклоны, чтобы найти перехват для этого сегмента.

set.seed(1) 
x1=rnorm(100,0,1) 
x2=rnorm(100,1,1) 
# Find points where x1 is above x2. 
above<-x1>x2 
# Points always intersect when above=TRUE, then FALSE or reverse 
intersect.points<-which(diff(above)!=0) 
# Find the slopes for each line segment. 
x1.slopes<-x1[intersect.points+1]-x1[intersect.points] 
x2.slopes<-x2[intersect.points+1]-x2[intersect.points] 
# Find the intersection for each segment. 
x.points<-intersect.points + ((x2[intersect.points] - x1[intersect.points])/(x1.slopes-x2.slopes)) 
y.points<-x1[intersect.points] + (x1.slopes*(x.points-intersect.points)) 
# Plot. 
plot(x1,type='l') 
lines(x2,type='l',col='red') 
points(x.points,y.points,col='blue') 

enter image description here

+0

есть проблемы с вертикальными сегментами, не так ли? – baptiste

+0

@baptiste Я думаю, что вертикальные сегменты невозможны. Возможно, вы имели в виду горизонтальные сегменты, потому что, да, горизонтальные перекрывающиеся сегменты создавали бы проблему. На самом деле, это также было бы проблемой, если точки пересекались в точно оцененной точке. Удобно, вы могли бы проверить оба одновременно: проверьте, равны ли точки в вашем векторе, добавьте их в набор пересекающихся точек, затем проверьте, произошли ли эти точки один за другим. Может, это не то, что вы имели в виду? – nograpes

+0

нет, я просто думал об общей проблеме пересечения сегментов сегмента, но здесь вы правы, ни один сегмент не является вертикальным (что вызовет проблемы с бесконечным наклоном). – baptiste

5

Вот альтернативный сегмент-сегмент кода пересечения,

# segment-segment intersection code 
# http://paulbourke.net/geometry/pointlineplane/ 
ssi <- function(x1, x2, x3, x4, y1, y2, y3, y4){ 

    denom <- ((y4 - y3)*(x2 - x1) - (x4 - x3)*(y2 - y1)) 
    denom[abs(denom) < 1e-10] <- NA # parallel lines 

    ua <- ((x4 - x3)*(y1 - y3) - (y4 - y3)*(x1 - x3))/denom 
    ub <- ((x2 - x1)*(y1 - y3) - (y2 - y1)*(x1 - x3))/denom 

    x <- x1 + ua * (x2 - x1) 
    y <- y1 + ua * (y2 - y1) 
    inside <- (ua >= 0) & (ua <= 1) & (ub >= 0) & (ub <= 1) 
    data.frame(x = ifelse(inside, x, NA), 
      y = ifelse(inside, y, NA)) 

} 
# do it with two polylines (xy dataframes) 
ssi_polyline <- function(l1, l2){ 
    n1 <- nrow(l1) 
    n2 <- nrow(l2) 
    stopifnot(n1==n2) 
    x1 <- l1[-n1,1] ; y1 <- l1[-n1,2] 
    x2 <- l1[-1L,1] ; y2 <- l1[-1L,2] 
    x3 <- l2[-n2,1] ; y3 <- l2[-n2,2] 
    x4 <- l2[-1L,1] ; y4 <- l2[-1L,2] 
    ssi(x1, x2, x3, x4, y1, y2, y3, y4) 
} 
# do it with all columns of a matrix 
ssi_matrix <- function(x, m){ 
    # pairwise combinations 
    cn <- combn(ncol(m), 2) 
    test_pair <- function(i){ 
    l1 <- cbind(x, m[,cn[1,i]]) 
    l2 <- cbind(x, m[,cn[2,i]]) 
    pts <- ssi_polyline(l1, l2) 
    pts[complete.cases(pts),] 
    } 
    ints <- lapply(seq_len(ncol(cn)), test_pair) 
    do.call(rbind, ints) 

} 
# testing the above 
y1 = rnorm(100,0,1) 
y2 = rnorm(100,1,1) 
m = cbind(y1, y2) 
x = 1:100 
matplot(x, m, t="l", lty=1) 
points(ssi_matrix(x, m)) 

enter image description here

0

Поздний ответ, но вот "пространственный" метод с использованием пакета SP и RGEOS. Это требует, чтобы оба x и y были числовыми (или могут быть преобразованы в числовые). Проекция произвольно, но EPSG: 4269, казалось, хорошо работать:

library(sp) 
library(rgeos) 
# dummy x data 
x1 = rnorm(100,0,1) 
x2 = rnorm(100,1,1) 

#dummy y data 
y1 <- seq(1, 100, 1) 
y2 <- seq(1, 100, 1) 

# convert to a sp object (spatial lines) 
l1 <- Line(matrix(c(x1, y1), nc = 2, byrow = F)) 
l2 <- Line(matrix(c(x2, y2), nc = 2, byrow = F)) 
ll1 <- Lines(list(l1), ID = "1") 
ll2 <- Lines(list(l2), ID = "1") 
sl1 <- SpatialLines(list(ll1), proj4string = CRS("+init=epsg:4269")) 
sl2 <- SpatialLines(list(ll2), proj4string = CRS("+init=epsg:4269")) 

# Calculate locations where spatial lines intersect 
int.pts <- gIntersection(sl1, sl2, byid = TRUE) 
int.coords <- [email protected] 

# Plot line data and points of intersection 
plot(x1, y1, type = "l") 
lines(x2, y2, type = "l", col = "red") 
points(int.coords[,1], int.coords[,2], pch = 20, col = "blue") 

 Смежные вопросы

  • Нет связанных вопросов^_^