You are here:GeoTux»Geo-Blogs»Librerías para Geomática»Ajustando puntos a líneas en R

Estadísticas

Invitados: 41
Usuarios registrados: 3131
Usuarios en línea:
-
Registrados hoy:
-

Registro

Redifusión (RSS)

Blogs y Noticias:
Recibe las actualizaciones en Geo-Noticias y Geo-Blogs

Recibir por e-mail
Recibir Geo-Noticias y Geo-Blogs por e-mail

¿Qué es esto?

En Twitter

Lunes 20 de Febrero de 2012 19:44

Ajustando puntos a líneas en R

Written by  German Carrillo
Rate this item
(0 votes)

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)



Last modified on Sábado 30 de Junio de 2012 19:56

Comentarios  

 
0 # ??Roger Bivand 18-04-2012 18:36
If you could see your way to trying to add these operations to rgeos, I'm sure that such a contribution would be widely appreciated. We've only ported some of what is possible through the C API, based on known use cases for us. If you'd like to be added as a developer on R-Forge, please let me know.
Responder | Responder con una citación | Citar
 
 
0 # Re:tuxman 18-04-2012 20:17
Well, if this code is valuable enough to be added to rgeos, I'd definitely do it. But I'm afraid porting the GEOS functions from C to R would take quite a long time, since I'm not familiar with it.
Responder | Responder con una citación | Citar
 
 
0 # codeRoger Bivand 19-04-2012 06:55
OK, so then it would make more sense to add your functions to the maptools package. Interfacing the GEOS library in C is what the rgeos package does. If you could send me source code and a help page, I can try them out and add them to maptools.
Responder | Responder con una citación | Citar
 
 
0 # Re:tuxman 19-04-2012 10:17
Alright, I'll do it as soon as I can.
Responder | Responder con una citación | Citar
 
 
+2 # Updatetuxman 30-06-2012 19:50
The functions published in this post have been included into the R package maptools! See the update [1]. (Thanks to Roger Bivand)

Tuxman
----
[1] http://geotux.tuxfamily.org/index.php/en/geo-blogs/item/296-snapping-points-to-lines-in-r#update
Responder | Responder con una citación | Citar
 
 
0 # Update #2tuxman 22-01-2014 19:57
ACTUALIZACIÓN #2 (2014.01.22):

A partir de ahora, la función snapPointsToLin es del paquete maptools acepta como parámetro el campo que contiene los IDs de las líneas.

Por ejemplo, si la capa de líneas tiene sus IDs almacenados en el campo "codigos", es posible indicarlo así:

Código:snappedPoints = snapPointsToLin es(p, l, idField="codigos")

De esta manera, la capa de puntos ajustados contendrá el ID de la línea a la que se ajustó cada punto, tomándolo del campo "codigos".

De no indicarse el parámetro idField, se tomará un ID interno del paquete sp, el cual es un consecutivo.

============================================
UPDATE #2 (2014.01.22):

From now on, it is possible to specify the ID field of the line data set by providing the parameter idField, this way (assuming the IDs are stored in the field "codes"):

Código:snappedPoints = snapPointsToLin es(p, l, idField="codes")

The snapped points data set will contain the ID of the line which each point was snapped to and will take it from the field "codes" of the line data set.

If no idField is specified, the line ID will be obtained from the sp package, which corresponds to a sequential number.
Responder | Responder con una citación | Citar
 
 
0 # errorzara 28-05-2016 14:45
Thanks for the helpful information.
However I couldn't get the result.
It shows the following error :


> proj4string(l) = crs
Warning message:
In `proj4string
Responder | Responder con una citación | Citar
 
 
0 # errorzara 30-05-2016 19:28
Seems my post is not complete. I meant This error :

proj4string(l) = crs
Warning message:
In `proj4string
Responder | Responder con una citación | Citar
 

Escribir un comentario


Código de seguridad
Refescar

 

¿Dónde nos leen?

Últimos comentarios