Proyección conforme cónica de Lambert

La proyección conforme cónica de Lambert es una de las muchas formas que existen de representar cartográficamente el mapa de la tierra sobre un plano.
Fue presentada por el matemático, físico, filósofo y astrónomo Johann Heinrich Lambert en 1772.
Esta proyección superpone un cono sobre la esfera de la Tierra, con dos paralelos de referencia secantes al globo.

¿Que tiene de especial esta proyección?
Pues que la distorsión es nula a lo largo de los paralelos de referencia y se incrementa fuera de ellos.
    El cuadriculado Lambert es nacional; y el punto cero de convergencia o meridiano cero, pasa por el Observatorio Astronómico de Madrid, en este punto la meridiana geográfica y la Lambert es la misma, pero a medida que nos alejamos al E o al W de Madrid varía el ángulo de convergencia.
Además es el tipo de proyección que se utiliza al realizar mapas en la mayoría de los países de Europa.


La formulación matemática para transformar un punto según un sistema de referencia geodésico a coordenadas según la proyección conforme cónica de Lambert es la siguiente:


donde:

donde  es la longitud,  la longitud de referencia,  la latitud,  la latitud de referencia y  y  los paralelos estándar.
Hay que tener en cuenta que la operación inversa es iterativa, es decir que si nos dan inicialmente una posición X,Y basada en PCCL y nos piden obtener la coordenada geográfica correspondiente tendremos que hacer lo siguiente (Con 5 iteraciones es suficiente para que converja)
....
  // ITERACION
  i1 := pi / 2 - 2 * ArcTAN(c);
  FOR i := 1 TO 5 DO
    BEGIN
      i10 := pi / 2 - 2 * ArcTAN(c * power(((1 - e * sin(i1)) / (1 + e * sin(i1))), (e / 2)));
      i1 := i10;
    END;
...
A continuación tenéis los dos procedimientos que he realizado para pasar coordenadas X,Y a Coordenadas Geográficas y viceversa.

{
PARÁMETROS INICIALES CON LOS QUE HE TRABAJADO
 (datos de la situación de España):
CONST
  FalsoEste = 500000;
  FalsoNorte = 0;
  LatitudReferencia = 40; 
  LongitudReferencia = -3; 
  RadioTierra = 6378137; 
  inversof = 298.257223563; // SI WGS84
  // Inversof = 298.257222101; // si ETRS89
  ParaleloNorte = 44; 
  ParaleloSur = 36; 
  pi = 3.141592653589793238;

}




PROCEDURE TfMapOffLine.PCCL_GeoAXY(long, lat: Double; VAR X, Y: Double);
VAR
  rLongRef, rLatRef: Double;
  n, F, Fi: Double;
  e2, m, m1, m2, t, t0, t1, t2, rho0, rho, k, gamma: Double;
  rLong, rLat, ParaleloStd1, ParaleloStd2: Double;
BEGIN
  rLatRef := DegToRad(LatitudReferencia); 
  rLongRef := DegToRad(LongitudReferencia); 
  ParaleloStd1 := DegToRad(ParaleloNorte); 
  ParaleloStd2 := DegToRad(ParaleloSur); 

  rLong := DegToRad(long); 
  rLat := DegToRad(lat); 

  e2 := 2 / inversof - 1 / power((inversof), 2); 
  m := cos(rLat) / SQRT(1 - e2 * power((sin(rLat)), 2)); 
  m1 := cos(ParaleloStd1) / SQRT(1 - e2 * power((sin(ParaleloStd1)), 2)); 
  m2 := cos(ParaleloStd2) / SQRT(1 - e2 * power((sin(ParaleloStd2)), 2)); 
  t := TAN(pi / 4 - (rLat / 2)) / power(((1 - SQRT(e2) * sin(rLat)) / (1 + SQRT(e2) * sin(rLat))), (SQRT(e2) / 2));
  
  t0 := TAN(pi / 4 - (rLatRef / 2)) / power(((1 - SQRT(e2) * sin(rLatRef)) / (1 + SQRT(e2) * sin(rLatRef))),
    (SQRT(e2) / 2));

  t1 := TAN(pi / 4 - (ParaleloStd1 / 2)) /
    power(((1 - SQRT(e2) * sin(ParaleloStd1)) / (1 + SQRT(e2) * sin(ParaleloStd1))), (SQRT(e2) / 2)); 
  t2 := TAN(pi / 4 - (ParaleloStd2 / 2)) /
    power(((1 - SQRT(e2) * sin(ParaleloStd2)) / (1 + SQRT(e2) * sin(ParaleloStd2))), (SQRT(e2) / 2));
  n := (LN(m1) - LN(m2)) / (LN(t1) - LN(t2)); 
  F := m1 / (n * power(t1, n)); 
  rho0 := RadioTierra * F * power(t0, n); 
  rho := RadioTierra * F * power(t, n); 
  Fi := n * (rLong - rLongRef); 
  X := rho * sin(Fi);
  Y := rho0 - rho * cos(Fi);
  k := (rho * n) / (RadioTierra * m);
  gamma := Fi;
  X := X + FalsoEste;
  Y := Y + FalsoNorte;

END;

PROCEDURE TfMapOffLine.PCCL_XYAGeo(XInic, yInic: Double; VAR Latitud, Longitud: Double);
VAR
  i: Integer;
  rLongRef, rLatRef, ParaleloStd1, ParaleloStd2: Double;
  n, F: Double;
  e, e2, m, m1, m2, t, t0, t1, t2, rho0, rho, X, Y: Double;
  rLong, rLat, long, lat: Double;
  a, b, c: Double;
  i1, i2, i3, i4, i5, i6, i7, i8, i9, i10: Double; // iteracion
BEGIN
  rLatRef := DegToRad(LatitudReferencia); 
  rLongRef := DegToRad(LongitudReferencia); 

  ParaleloStd1 := DegToRad(ParaleloNorte); 
  ParaleloStd2 := DegToRad(ParaleloSur);

  rLong := DegToRad(LongitudReferencia); 
  rLat := DegToRad(LatitudReferencia); 

  e2 := 2 / inversof - 1 / power((inversof), 2); 
  e := SQRT(e2);
  X := XInic - FalsoEste; 
  Y := yInic - FalsoNorte; 

  m := cos(rLat) / SQRT(1 - e2 * power((sin(rLat)), 2)); 
  m1 := cos(ParaleloStd1) / SQRT(1 - e2 * power((sin(ParaleloStd1)), 2)); 
  m2 := cos(ParaleloStd2) / SQRT(1 - e2 * power((sin(ParaleloStd2)), 2)); 
  t := TAN(pi / 4 - (rLat / 2)) / power(((1 - SQRT(e2) * sin(rLat)) / (1 + SQRT(e2) * sin(rLat))), (SQRT(e2) / 2));
  
  t0 := TAN(pi / 4 - (rLatRef / 2)) / power(((1 - SQRT(e2) * sin(rLatRef)) / (1 + SQRT(e2) * sin(rLatRef))),
    (SQRT(e2) / 2));
  
  t1 := TAN(pi / 4 - (ParaleloStd1 / 2)) /
    power(((1 - SQRT(e2) * sin(ParaleloStd1)) / (1 + SQRT(e2) * sin(ParaleloStd1))), (SQRT(e2) / 2)); 
  t2 := TAN(pi / 4 - (ParaleloStd2 / 2)) /
    power(((1 - SQRT(e2) * sin(ParaleloStd2)) / (1 + SQRT(e2) * sin(ParaleloStd2))), (SQRT(e2) / 2)); 
  n := (LN(m1) - LN(m2)) / (LN(t1) - LN(t2)); 
  F := m1 / (n * power(t1, n)); 
  rho0 := RadioTierra * F * power(t0, n); 
  rho := SQRT(power(X, 2) + power((rho0 - Y), 2)); 
  c := power((rho / (RadioTierra * F)), (1 / n)); 
  a := pi / 2 - 2 * ArcTAN(c); 
  b := ArcTAN(X / (rho0 - Y)); 
  Longitud := b / n + rLongRef;
  
  // ITERACION
  i1 := pi / 2 - 2 * ArcTAN(c);
  //con 5 iteraciones es suficiente 
  FOR i := 1 TO 5 DO
    BEGIN
      i10 := pi / 2 - 2 * ArcTAN(c * power(((1 - e * sin(i1)) / (1 + e * sin(i1))), (e / 2)));
      i1 := i10;
    END;

  Latitud := i10 * 180 / pi; 
  Longitud := Longitud * 180 / pi;
END;


Si os interesa ver la aplicación de esta proyección y más temas relacionados con la medida de longitudes, posición de satélites, medida de áreas, cálculo de rutas, mapas offline, etc... os recomiendo descargaros la app GPS TOTAL RUN (en un principio la hice pensando en los runners que les encanta el deporte, con módulos para mostrar la ruta realizada, hacer estadísticas, etc..., pero luego ya ví que podía mejorarla con mas utilidades)