Descuento del 40% en productos de Embarcadero

Limited Time Offer

Sólo hasta el 31 de Diciembre de 2017 tienes un 40% de descuento en las versiones de DELPHI, C++ BUILDER y RAD STUDIO Professional, Entreprise y Architect 



Aprovecha que quedan pocas horas para hacerte con esta oferta.
Aquí tienes el link donde, si te interesa, los puedes comprar on-line. 


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;


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)