En esta oportunidad presento un algoritmo en R para ajustar puntos a líneas (en inglés conocido como “snapping”). La aplicación de ejemplo es la localización de puntos de trayectorias (obtenidos por ejemplo desde navegadores GPS) en una red vial.
For English click here.
Ajustar puntos a líneas (en inglés “snap”) es una tarea frecuente cuando se trabaja con puntos que por alguna razón deben estar localizados sobre líneas.
Una de las aplicaciones más comunes es ajustar datos de trayectorias a una red vial. Como tal, se trata de un problema de georreferenciación que busca mejorar las posiciones de las trayectorias al eliminar sus desviaciones utilizando una red vial como referencia.
Este problema, también conocido como “point-to-curve matching” (algo así como emparejamiento de punto a curva), es uno de los más simples enfoques para construir algoritmos para “map-matching” (emparejamiento de mapas), los cuales pueden incluir topología, probabilidad y lógica difusa, entre otros. Una revisión completa de algoritmos para “map-matching” se da en (Quddus et al., 2007).
En cuanto a software libre para geomática (FOSS4G), librerías como Java Topology Suite (JTS) y su traducción a C, GEOS, proveen clases y funciones que pueden ser usadas para ajustar puntos a líneas. Algunos ejemplos que muestran su uso son descritos por P. Ramsey (usando SQL y funciones PostGIS), Marcello Benigno (usando GRASS) y la guía de usuario de GeoTools.
Sin embargo, el paquete de R “rgeos” (interfaz de R para GEOS) aún no provee dichas funciones y uno debe implementarlas por si mismo. En lo que queda de este post les mostraré mi solución, que ajusta los puntos de las trayectorias a los puntos más cercanos que caen sobre la red vial, en términos geométricos, se hace una proyección de los puntos sobre las líneas.
Primero, se debe importar una red vial a R (pueden descargar un shapefile desde este enlace, fuente: OpenStreetMap), no olviden ajustar la ruta a los datos:
1 2 3 4 5 |
require("rgdal") crs = CRS("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs") l = readOGR(dsn="/docs/geodata/munster/osm/", layer="subset_street_graph") proj4string(l) = crs |
Importar los datos de la trajectoria desde http://ifgi.uni-muenster.de/~epebe_01/pm10/All4_1a.txt y asegurarse de que ambos (red vial y trayectoria) compartan el mismo sistema de referencia proyectado (para el cálculo de distancias):
6 7 8 9 10 11 12 |
x = read.csv("http://ifgi.uni-muenster.de/~epebe_01/pm10/All4_1a.txt") x$t = as.POSIXct(strptime(do.call(paste, x[c("Date", "Time")]), "%Y/%m/%d %H:%M:%S")) m = matrix(c(x[,2],x[,1]),ncol=2) df = data.frame(x$t, x[,5],x[,6]) p = SpatialPointsDataFrame(coords=m,data=df,proj4string=CRS("+proj=longlat +datum=WGS84")) p = spTransform(p,crs) |
Definir las funciones nearestPointOnSegment, nearestPointOnLine, de esta forma:
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
nearestPointOnSegment = function(s, p){ # Adapted from http://pastebin.com/n9rUuGRh ap = c(p[1] - s[1,1], p[2] - s[1,2]) ab = c(s[2,1] - s[1,1], s[2,2] - s[1,2]) t = sum(ap*ab) / sum(ab*ab) t = ifelse(t<0,0,ifelse(t>1,1,t)) x = s[1,1] + ab[1] * t y = s[1,2] + ab[2] * t c(x, y, (x-p[1])^2 + (y-p[2])^2) # Return nearest point and distance } nearestPointOnLine = function(coordsLine, coordsPoints){ nearest_points = vapply(2:nrow(coordsLine), function(x) nearestPointOnSegment(coordsLine[(x-1):x,], coordsPoints), FUN.VALUE=c(0,0,0)) # Return coordinates of the nearest point in this line nearest_points[1:2, which.min(nearest_points[3,])] } |
Finalmente, definir la función principal, llamada snapPointsToLines:
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 |
snapPointsToLines <- function( points, lines, maxDist=NA, withAttrs=TRUE) { require("rgeos") if (!is.na(maxDist)) { w = gWithinDistance(points, lines, dist=maxDist, byid=TRUE) validPoints = apply(w,2,any) validLines = apply(w,1,any) points = points[validPoints,] lines = lines[validLines,] } d = gDistance(points, lines, byid=TRUE) nearest_line_index = apply(d, 2, which.min) # Position of each nearest line in lines object coordsLines = coordinates(lines) coordsPoints = coordinates(points) # Get coordinates of nearest points lying on nearest lines mNewCoords = vapply(1:length(points), function(x) nearestPointOnLine(coordsLines[[nearest_line_index[x]]][[1]], coordsPoints[x,]), FUN.VALUE=c(0,0)) # Recover lines' Ids (Ids and index differ if maxDist is given) if (!is.na(maxDist)) nearest_line_id = as.numeric(rownames(d)[nearest_line_index])+1 else nearest_line_id = nearest_line_index # Create data frame and sp points if (withAttrs) df = cbind(points@data, nearest_line_id) else df = data.frame(nearest_line_id, row.names=names(nearest_line_index)) SpatialPointsDataFrame(coords=t(mNewCoords), data=df, proj4string=CRS(proj4string(points))) } |
Ahora pueden usar la función snapPointsToLines con los datos completos:
68 |
snappedPoints = snapPointsToLines(p, l) |
O pueden establecer una distancia máxima (en este caso en metros) para evitar ajustar puntos que se encuentran más lejos:
69 |
snappedPoints = snapPointsToLines(p, l, 100) |
También pueden prescindir de los atributos originales y obtener únicamente las coordenadas de los puntos ajustados y el identificador de su línea más cercana:
70 |
snappedPoints = snapPointsToLines(p, l, 100, FALSE) |
Esto es lo que obtengo ejecutando la línea 69 y llamando la función plot de esta manera (las cruces rojas son los puntos de la trayectoria original y las verdes son los puntos ajustados):
71 72 73 |
plot(p,col="red") plot(snappedPoints,add=TRUE,col="green") plot(l,add=TRUE) |
Una vista más detallada (las cruces rojas son los puntos de la trayectoria original y las azules son los puntos ajustados):
Una manera sencilla de mejorar el resultado es pasar como parámetro una red vial filtrada, por ejemplo, solo con vías principales. Por supuesto que esto dependerá de la trayectoria. Este es el resultado de ajustar los puntos a la red vial filtrada (usando 100 metros como distancia máxima):
La solución presentada puede ser mejorada en términos de rendimiento. Con los datos de este post, la versión actual es casi dos veces más lenta que la de PostGIS.
REFERENCIAS
- M. Quddus, W. Ochieng and R. Noland. Current map-matching algorithms for transport applications : State-of-the art and future research directions. In Elsevier, Kidlington, ROYAUME-UNI, editor, Transportation research. Part C, Emerging technologies, vol. 15, no5, pages 312-328, 2007.
ACTUALIZACIÓN (2012.06.30)
Nos complace anunciar que las funciones snapPointsToLines, nearestPointOnSegment y nearestPointOnLine han sido incluidas en el reconocido paquete ‘maptools‘ (v. 0.8-16) de R.
Las funciones pueden ser llamadas directamente después de importar maptools a la sesión de R, de esta forma:
# install.packages("maptools") # Para instalar el paquete require("maptools") # Crear las líneas y los puntos (ver líneas 1-5 y 6-12 en el post) # Finalmente, llamar la función principal snappedPoints = snapPointsToLines(p, l, 100)
Comentarios recientes