Compruebe si la lat-long caer en mi ruta

Question

Estoy usando ggmap paquete para generar rutas. Los datos resultantes marco contendrá las rutas de la información desde la fuente hasta el destino.

   library(ggmap)
    route_df <- route(from = "Cambridge, 02138, USA",
              to = "Massachusetts Ave, 02139, USA",
              mode = "driving", structure = "route", alternatives = TRUE)

Quiero verificar si un determinado coordinar (-71.13857, 42.38039) caen en esta ruta o no. También quiero comprobar si una coordenada (-71.13857, 42.38039) es parte (uno de los waypoints de la ruta "A" resultado de la "route_df" de la variable. ¿Hay algún paquete que puedo utilizar para hacerlo?

Answer 1

Pruebe la reproducibles ejemplo a continuación. Usted necesidad de utilizar objetos Espaciales.

# Load libraries
library('ggmap')
library('sp')
library('rgeos')
library('mapview') # Interactive maps in R

# Grab a route from Google
route_df <- route(from = "Cambridge, 02138, USA",
                  to = "Massachusetts Ave, 02139, USA",
                  mode = "driving", structure = "route", alternatives = TRUE)


# Find UTM projection
# https://stackoverflow.com/questions/9186496/determining-utm-zone-to-convert-from-longitude-latitude
long2UTM <- function(long) {
  (floor((long + 180)/6) %% 60) + 1
}

long2UTM(route_df$lon)

# Convert route points to SpatialPointsDataFrame
routeSPDF <- route_df
coordinates(routeSPDF) <- c("lon", "lat")
proj4string(routeSPDF) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") # original CRS
routeSPDF <- spTransform(x = routeSPDF, CRSobj = CRS("+proj=utm +zone=19 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) # Transform (project)

# Convert routes to SpatialLinesDataFrame
possibleRoutes <- unique(routeSPDF$route)

# List of lines (one line is one possible route)
linesList <- list()

for(i in 1:length(possibleRoutes)) {
  linesList[[i]] <- Lines(slinelist = list(Line(coordinates(routeSPDF[which(routeSPDF$route == possibleRoutes[i]),]))), ID = possibleRoutes[i])
}

# SpatialLines
routesSL <- SpatialLines(LinesList = linesList, proj4string = CRS("+proj=utm +zone=19 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

# Create particular coordinate as SpatialPoints object
point1 <- SpatialPoints(coords = cbind(-71.13857, 42.38039), proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
point1 <- spTransform(x = point1, CRSobj = CRS("+proj=utm +zone=19 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) # Transform (project)

# Measure distance from points to lines (in meters)
gDistance(spgeom1 = point1, spgeom2 = routesSL, byid = TRUE)

# Verify if point1 fall in this route or not
gIntersects(spgeom1 = point1, spgeom2 = routesSL, byid = TRUE)
gWithinDistance(spgeom1 = point1, spgeom2 = routesSL, byid = TRUE, dist = 13.7)

# Verify if point1 is part (one of the waypoints) of the route "A" resulted from the "route_df"
gWithinDistance(spgeom1 = point1, spgeom2 = routesSL["A"], byid = TRUE, dist = 13.7)

# Test with a point from the line
testPoint <- routeSPDF[5,]
gDistance(spgeom1 = testPoint, spgeom2 = routesSL, byid = TRUE)
gWithinDistance(spgeom1 = testPoint, spgeom2 = routesSL, byid = TRUE, dist = 8.703543e-10)

# Plot
mapview(routeSPDF, zcol = 'route') + mapview(routesSL, zcol = 'ID') + mapview(point1, col = "black") + mapview(testPoint)