在R中找到直线和三角形之间的交点

Finding intersection points between a line and a triangle in R

我试图找到路径和三角形的交点。输出可以是具有一组与三角形相交的路径的经度或纬度点的简单矩阵。这些点可以在图中突出显示。谁能帮我解决这个问题?

# path 01 
longitude <- round(c(48.7188021250007,
                 48.7188133749999,
                 48.7188291249998,
                 48.7188336250004), 8);

Latitude <- round(c (2.39514523661229,
                 2.39512477447308, 
                 2.39472235230961,
                 2.39467460730113), 8);

# path 02
longitude_2 <- round(c(48.71881,
                   48.71882,
                   48.71883,
                   48.71884), 8);

Latitude_2 <- round(c (2.39479,
                   2.3947,
                   2.39469396980779,
                   2.3945), 8);

plot(longitude,Latitude,ylim=range(c(Latitude,Latitude_2)),
xlim=range(c(longitude,longitude_2)), type="l",col="blue")
lines(longitude_2,Latitude_2,col="darkgreen")
title(main = "Two different paths");
legend("right", legend=c("path A", "path B"),
   col=c("blue", "darkgreen"), lty=1:1, cex=0.8, 
   box.lty=0)


# Drawing polygon 

lon <- c(48.71882,48.71883, 48.71884);
lat <- c (2.3945, 2.39478,  2.39479);
x <- cbind(lon, lat)
polygon(x, col='lightblue')

我认为您可以在创建 SpatialLines 和 SpatialPolygons 后使用 rgeos 包中的 gIntersection 函数。

require(sp)
require(rgeos)
longitude <- round(c(48.7188021250007,
                     48.7188133749999,
                     48.7188291249998,
                     48.7188336250004), 8);

Latitude <- round(c (2.39514523661229,
                     2.39512477447308, 
                     2.39472235230961,
                     2.39467460730113), 8);

longitude_2 <- round(c(48.71881,
                       48.71882,
                       48.71883,
                       48.71884), 8);

Latitude_2 <- round(c (2.39479,
                       2.3947,
                       2.39469396980779,
                       2.3945), 8);

lon <- c(48.71882,48.71883, 48.71884);
lat <- c (2.3945, 2.39478,  2.39479);


lines <- SpatialLines(list(Lines(list(Line(cbind(longitude, Latitude)),
                                      Line(cbind(longitude_2, Latitude_2))), 1)))

pol_1 <- SpatialPolygons(list(Polygons(list(Polygon(coords = cbind(lon, lat))),1)))

result <- gIntersection(lines, pol_1)

plot(lines, col = "red")
plot(pol_1, col = "blue", add = T)
plot(result, col = "green", add = T)

或者,如果您想要具有正确缩放比例的相同绘图,您可以使用以下方法重新设计它:

plot(x = longitude,
     y = Latitude,
     ylim = range(c(Latitude,Latitude_2)),
     xlim = range(c(longitude,longitude_2)), 
     type="l",
     col="blue")
lines(longitude_2,Latitude_2,
      col = "darkgreen")
polygon(pol_1@polygons[[1]]@Polygons[[1]]@coords, 
        col = 'lightblue')
title(main = "Two different paths");
legend("right", 
       legend=c("path A", "path B", "Intersect"),
       col = c("blue", "darkgreen", "red"), 
       lty = 1:1, 
       cex = 0.8, 
       box.lty = 0)

lines(result@lines[[1]]@Lines[[1]]@coords, col = "red")
lines(result@lines[[1]]@Lines[[2]]@coords, col = "red")
points(result@lines[[1]]@Lines[[1]]@coords, col = "red")
points(result@lines[[1]]@Lines[[2]]@coords, col = "red")
text(result@lines[[1]]@Lines[[1]]@coords, labels = pasteCols(round(t(result@lines[[1]]@Lines[[1]]@coords),4), sep = " "), cex= 0.7)
text(result@lines[[1]]@Lines[[2]]@coords, labels = pasteCols(round(t(result@lines[[1]]@Lines[[1]]@coords),4), sep = " "), cex= 0.7)