2011-09-26

Algoritmo para calcular la distancia entre dos puntos en la tierra

En algún lugar de la red me encontré con este algoritmo para calcular la distancia entre dos puntos del planeta, lo dejo para su conocimiento y utilidad:

Calcular la distancia entre dos puntos sobre un plano podria llegar a ser relativamente sencillo. Sin embargo, cuando estos dos puntos los ubicamos sobre la esfera terrestre, es decir lo que pretendemos es calcular la distancia lineal entre dos posiciones dadas (latitud + longitud) la cosa se complica.
Básicamente se complica por que en el calculo de la distancia entre ambas posiciones debemos contemplar la curvatura terrestre. Es aqui dónde entra en escena la Formula del Haversine.
Sin entrar en demasiados detalles en términos matemáticos, la Formula del Haversine es:
R = radio de la Tierra
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
Para utilizar la Formula del Haversine necesitamos, además de las dos posiciones (lat + lon) el radio de la Tierra. Este valor es relativo a la latitud pues al no ser la Tierra perfectamente redonda, el valor del radio ecuatorial es de 6378km mientras que el polar es de 6357km. El radio equivolumen es de 6371km. Para este post utilizaremos el valor del radio equatorial.
Por tanto lo primero que haremos es representar la clase Posicion que no es más que una clase con un par de propiedades relativas a la longitud y latitud.
public class Posicion{

    public Posicion(float latitud, float longitud){         
        Latitud = latitud;         
        Longitud = longitud;     
    }     

    public float Latitud {
        get;
        set; 
    }     

    public float Longitud { 
        get; 
        set; 
    }
}
Hablando en C#, lo que pretendo es que dados dos objetos del tipo Posicion pueda hacer algo tal que así para calcular la distancia entre ambos. Si tenemos el objeto Igualada representando la ubicación georeferenciada de la ciudad y queremos saber la distancia lineal con Granada haremos lo siguiente:
[Test]

public void Test(){  
   
    var Igualada = new Posicion(41.57879F,  1.617221F);     
    var Granada = new Posicion(37.176487F, -3.597929F);     
    float distancia = Igualada.DistanciaKm(Granada);     
    Assert.AreEqual(664.0D, System.Math.Round(distancia));
}
Para ello vamos a crear un método extensor llamado DistanciaKm sobre la clase Posicion el cual tendra la responsabilidad de calcular la distancia entre ambas posiciones; en definitiva contendrá la Fórmula del Haversine:
public static class Extensiones {
     public const float RadioTierraKm = 6378.0F;     

     public static float DistanciaKm(this Posicion posOrigen, Posicion posDestino){ 
         //TODO   
     }
}
  
Ahora pasamos ya a las operaciones matemáticas, el primer calculo que debemos hacer es obtener la diferencia entre las latitudes y longitudes de ambas posiciones. Si P1 con Lat1 y Lon2 y P2 con Lat2 y Lon2 entonces tendremos algo así:
DiferenciaLats = P2.Lat2 – P1.Lat1
DiferenciaLongs = P2.Lon2 – P2.Lon2
El primer “inconveniente” que encontramos es que las latitudes y longitudes se exponen en grados / minutos / segundos y por tanto debemos pasarlas a Radianes. Si pasamos todo esto a código dentro del método extensor DistanciaKm tenemos:
public static float DistanciaKm(this Posicion posOrigen, Posicion posDestino){
     
    var difLatitud = (posDestino.Latitud – posOrigen.Latitud).EnRadianes();        
    var difLongitud = (posDestino.Longitud -posOrigen.Longitud).EnRadianes();
}

Como veis hemos vuelto a hacer uso de otro método extensor, esta vez sobre System.Single o float y no es más que:

static float EnRadianes(this float valor){
    return Convert.ToSingle(Math.PI / 180) * valor; 
}
El segundo paso es calcular la mitad del cuadrado de la distancia en línea recta (acorde a la longitud) entre los dos puntos y que vamos a representar en la variable implicitamente tipada “a” y la distancia del angulo en radianes representada por “c”:
public static float DistanciaKm(this Posicion posOrigen, Posicion posDestino){

     var difLatitud = (posDestino.Latitud – posOrigen.Latitud).EnRadianes();     
          var difLongitud = (posDestino.Longitud -posOrigen.Longitud).EnRadianes();
          var a = Math.Sin(difLatitud/2).AlCuadrado() + Math.Cos(posOrigen.Latitud.EnRadianes())*     Math.Cos(posDestino.Latitud.EnRadianes())* Math.Sin(difLongitud/2).AlCuadrado(); var c= 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 – a)); 
}

En este último punto he vuelto a utilizar otro método extensor sobre el tipo double llamado AlCuadrado cuya firma es:

static double  AlCuadrado(this double valor){     
    return Math.Pow(valor, 2);
}
Finalmente, la distancia en Km será dada por el valor obtenido en “c” y multiplicado por el radio de la tierra, con lo que el método DistanciaKm quedará así:

public static float DistanciaKm(this Posicion posOrigen, Posicion posDestino){

     var difLatitud = (posDestino.Latitud – posOrigen.Latitud).EnRadianes();             
     var difLongitud = (posDestino.Longitud -posOrigen.Longitud).EnRadianes();
     var a = Math.Sin(difLatitud/2).AlCuadrado() + Math.Cos(posOrigen.Latitud.EnRadianes())* 
Math.Cos(posDestino.Latitud.EnRadianes())* Math.Sin(difLongitud/2).AlCuadrado(); var c= 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 – a)); 

         return RadioTierraKm*Convert.ToSingle©; 
}
A tener en cuenta
Pese a que la Formula del Haversine es de las mas utilizadas para el cáculo de distancias entre dos puntos (además de ésta también se utiliza la Ley Esférica del Coseno, entre otras) tened en cuenta que:
  • La formula asume que la Tierra es completamente redonda, con lo que cabe esperar una tasa de error que se podria llegar a asumir.
  • Durante la redacción de este post he encontrado información en la que afirma que es uso de la Formula del Haversine es menos precisa para calculo entre dos posiciones cuya distancia es inferior a 20 km (12millas).
  • Para el calculo de la distancia en millas seria necesario representar el radio de la Tierra en millas nauticas.