Delphi Magic: Proyección conforme cónica de Lambert

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:

$x=\rho \sin[n(\lambda -\lambda _{0})]$

$y=\rho _{0}-\rho \cos[n(\lambda -\lambda _{0})]$

donde:

$n={\frac {\ln(\cos \phi _{1}\sec \phi _{2})}{\ln[\tan({\frac {1}{4}}\pi +{\frac {1}{2}}\phi _{2})\cot({\frac {1}{4}}\pi +{\frac {1}{2}}\phi _{1})]}}$

$\rho =F\cot ^{n}({\frac {1}{4}}\pi +{\frac {1}{2}}\phi )$

$\rho _{0}=F\cot ^{n}({\frac {1}{4}}\pi +{\frac {1}{2}}\phi _{0})$

$F={\frac {\cos \phi _{1}\tan ^{n}({\frac {1}{4}}\pi +{\frac {1}{2}}\phi _{1})}{n}}$

donde

\lambda

es la longitud,

\lambda _{0}

la longitud de referencia,

\phi

la latitud,

\phi _{0}

la latitud de referencia y

\phi _{1}

\phi _{2}

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;

Tener en cuenta que si se van a combinar en un SIG datos procedentes de mapas topográficos (datum europeo) con posiciones
tomadas con GPS (datum WGS-84) es necesario establecer la correspondencia entre ambos. Las posiciones
tomadas con GPS deberán ser desplazadas 0.07 minutos al Norte y 0.09 minutos al Este.

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)

Delphi Magic

Proyección conforme cónica de Lambert

No hay comentarios:

Publicar un comentario

Simulación del movimiento de los electrones en un campo electrico