Puntos intermedios en un arco de círculo máximo

Hace varias semanas recibimos en el foro de GeoTux una pregunta relacionada con el cálculo de distancias a partir de dos puntos definidos por coordenadas geográficas, latitud y longitud. Como el tema es interesante, decidí publicar la respuesta en el blog.

La pregunta parafraseada es: Dados dos puntos A y B definidos por sus coordenadas geográficas latitud y longitud (asumimos en WGS84), cómo obtener los puntos intermedios con una distancia de x metros entre sí.

Dependiendo de la distancia entre A y B, hay diferentes modelos de la Tierra para realizar los cálculos (sin tener en cuenta el relieve): 1) Plano, 2) Esfera o elipsoide.

1) Para distancias relativamente cortas, se puede considerar la Tierra como un plano. Como se parte de coordenadas geográficas se pueden usar dos enfoques, prescindir de la curvatura o utilizar proyecciones:

a) Usar trigonometría en el plano sobre un triángulo definido por los puntos A, B y C (C tendría la latitud de B y la longitud de A o viceversa) con lados m (diferencia de longitudes), p (diferencia de latitudes) y d (distancia del recorrido), y un ángulo CAB que corresponde al rumbo de salida (en el caso del plano, el rumbo siempre es el mismo). Las distancias vienen dadas por la conversión 1′ de arco = 1 milla naútica (nmi) = 1852 m. Según (Ayres, 1970), los triángulos con distancias menores que 200 nmi pueden considerarse planos.

Triángulo esférico y plano

Usando esta opción se debe calcular el rumbo (conociendo los tres lados, m, p y d) y luego se debe cambiar el problema: a partir del rumbo, la posición de A y la distancia x (desde A hasta el punto X sobre la línea AB), se calcula un nuevo m para conocer la longitud del punto X y un nuevo p para conocer su latitud. Este último paso se realizaría variando la distancia x tantas veces como quepa en d (distancia AB), esto es, primero se toma x, luego 2x, luego 3x, etc.

b) Convertir los datos a un Sistema de Coordenadas Proyectadas (en metros), como el UTM, y allí usar herramientas que funcionen con datos en el plano, por ejemplo, segmentación dinámica o funciones como ST_Segmentize de PostGIS [1]. Esta última solución se puede encontrar implementada en [2], en donde se proyectan los puntos A y B a un sistema de coordenadas planas, se crea una línea con base en los puntos proyectados y finalmente se divide la línea, todo esto en una consulta SQL. Se debe tener en cuenta que las proyecciones implican ciertas deformaciones geométricas (distancias, áreas, ángulos) de la superficie terrestre, por lo que al medir distancias sobre zonas extensas se puede estar incurriendo en errores. Una discusión más amplia de las consideraciones principales sobre el cálculo distancias usando proyecciones puede encontrarse en (Iliffe and Lott, 2008), en donde se indica que para lograr mayor precisión, debe tenerse en cuenta el factor de escala (k), uno de los parámetros que definen las proyecciones pero que a veces varía dependiendo de las direcciones. Para facilitar las cosas, puede usarse la proyección conforme, en la que se preservan formas (y por ende ángulos) y en las que el factor de escala a lo largo de los meridianos es igual al factor de escala a lo largo de los paralelos. Cabe notar que al utilizar proyecciones es posible retornar a las coordenadas geográficas originales, por lo que en un sentido estricto no se está prescindiendo de la curvatura terrestre, como en el caso a).

Por su sencillez computacional, la proyección cilíndrica equidistante es usada para aproximar el cálculo de distancias sobre la esfera [3] o sobre el elipsoide [4], esta última referencia menciona además el cálculo del rumbo usando trigonometría en el plano.

2) Cuando se considere la Tierra como una esfera o un elipsoide, se debe tener en cuenta que la primera es menos precisa (suficiente para muchas tareas prácticas) pero más sencilla en términos matemáticos, lo que repercute en la rapidez de una implementación informática.

La pregunta que motiva este post corresponde a un par de problemas de la navegación: Primero, dado un punto A (punto inicial), un rumbo inicial y una distancia, hallar el punto X (punto de llegada). Segundo, dados los puntos A y B encontrar la distancia entre ellos y el rumbo inicial. Estos problemas son conocidos como los problemas directo e inverso de la Geodesia, respectivamente.  

En este caso, no se usa la línea recta como la mínima distancia entre dos puntos, sino que se usa el geodésico, que en la esfera corresponde a un arco de círculo máximo (ver ejemplo gráfico en [5]). La trigonometría esférica da herramientas para realizar este tipo de cálculos sobre la esfera, que corresponden a resolver un triángulo esférico con uno de sus vértices coincidiendo con uno de los polos terrestres (Ayres, 1970). Por otro lado, la geodesia ofrece soluciones para cálculos sobre los elipsoides por medio de integrales de ecuaciones diferenciales para las que Torge (1991, pp.219-221) menciona tres grupos de soluciones. Según Torge (1991, p.218), para distancias menores de 100 km. se puede  utilizar la esfera en lugar del elipsoide.

En cuanto a una implementación, pueden emplearse las fórmulas y aplicaciones en Javascript para la esfera dadas en [6], [7] y [8], y para el elipsoide en [9], [10] (problema directo) y [11] (problema inverso). El paquete Geosphere [12] del software R también las implementa.

En Geosphere hay 4 funciones para cálculos de distancias sobre círculos máximos que, en orden de complejidad son: ’ley de los cosenos esféricos’ [13], ’Haversine’, ’Vincenty Sphere’ y ’Vincenty Ellipsoid’. Las primeras tres usan la Tierra como esfera mientras que la última usa el elipsoide.

Ejemplo usando el paquete Geosphere de R:
Como conocemos los puntos A y B se necesita calcular el rumbo de partida (en el caso de los círculos máximos varía a medida que se navega) para posteriormente calcular la distancia AB. Luego se automatiza la tarea de obtener las coordenadas del punto X situado a x distancia del punto de partida A sobre el círculo máximo descrito, dado el rumbo calculado. Dicha automatización se puede realizar en R, mediante el paquete Geosphere de esta forma (coordenadas geográficas en grados decimales, distancia de separación en metros):

 require(geosphere)
 puntos_intermedios<-function(lon1, lat1, lon2, lat2, dist) {
     p1<-c(lon1, lat1)
     p2<-c(lon2, lat2)
     d<-dist
     dp1p2<-distVincentyEllipsoid(p1, p2)
     b<-bearing(p1,p2)
 
     count<-1
     x<-c()
     y<-c() 
     m<-c()
     while (d<dp1p2){             
         p<-destPoint(p1, b, d=d, r=6378137)        
         x<-c(x, p[1])
         y<-c(y, p[2])
         m<-c(m, d)
         count<-count+1
         d<-dist*count
     }
     points <- matrix(c(x,y,m),ncol=3)
     return(points)
 }
 p<-puntos_intermedios(-74.138332, 4.703056, 7.684831, 52.134642, 300000)

Obteniendo finalmente las coordenadas deseadas (longitud, latitud, distancia):

Coordenadas y distancias

Que en un mapa lucen así:

 require(maps)
 data(world.cities)
 map("world", fill=TRUE, col="#a2d75e")
 map.axes()
 points(p[,1],p[,2],col="red",pch=20,cex=0.5)
 points(c(-74.138332, 7.684831),c(4.703056, 52.134642),col="yellow",pch=20)
 title("Puntos intermedios en un arco de círculo máximo","P1: 74°08'18\\"W,4°42'11\\"N  P2: 7°41'05\\"E, 52°08'05\\"N", cex.sub=0.8)
 text(c(-74.138332, 7.684831),c(4.703056, 52.134642),labels=c("P1","P2"),col="yellow",pos=4,offset=0.6)
 text(c(-74.138332, 7.684831),c(4.703056, 52.134642),labels=c("P1","P2"),col="magenta",pos=4)

Puntos intermedios en un arco de círculo máximo

REFERENCIAS

  • Frank, Ayres. Teoría y problemas de Trigonometría plana y esférica. Traducción y adaptación de Antonio Linares. Serie de compendios Shaum. Ed. MacGraw-Hill. 1970.
  • Iliffe, J. and Lott, R. Datums and map projections: For remote sensing, GIS and surveying. Section 3.8. pp.86-89. 2008
  • Torge, Wofgang. Geodesy. 2nd edition. Transl. from the German by Christopher Jekeli. Berlin, New York: de Gruyter. 1991.

[1] http://www.postgis.org/docs/ST_Segmentize.html
[2] http://underdark.wordpress.com/2011/08/20/visualizing-global-connections/
[3] http://www.movable-type.co.uk/scripts/latlong.html#equirectangular
[4] http://williams.best.vwh.net/avform.htm#flat
[5] http://www.gcmap.com/mapui?P=BOG-FMO&MS=bm
[6] http://williams.best.vwh.net/avform.htm

[7] http://www.movable-type.co.uk/scripts/latlong.html
[8] http://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
[9] http://williams.best.vwh.net/gccalc.htm
[10] http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html
[11] http://www.movable-type.co.uk/scripts/latlong-vincenty.html
[12] http://cran.r-project.org/web/packages/geosphere/index.html
[13] http://williams.best.vwh.net/avform.htm#Triangle