Método para convertir coordenadas geográficas en UTM: Clase Coordenadas (C/C++)

En el artículo anterior, presenté un método para convertir coordenadas geocéntricas en geográficas (cambio de datum 7 parámetros). Este método es complementario a la transformación que se puede realizar del cambio de coordenadas UTM a UTM, con cambio de datum, que pasaría por los siguientes pasos: de UTM a geográficas, de geográficas a geocéntricas, de geocéntricas a geográficas (cambio de datum 7 parámetros o de Bursa-Wolf) y de geográficas a UTM; todos métodos de la Clase Coordenadas. Ahora me toca implementar el último de ellos, es decir, de geográficas a UTM.

El código es el siguiente:

void Coordenadas::

	CambioGeoUtm (int selector, char long_, char lat_, Coordenadas v, Coordenadas &t){

	const double pi = 3.14159265358979323846;

	double a, b;

	switch (selector) {

	case 1:

		a = 6378388.0, b = 6356911.946130;    //Hayford
	
                break;

	case 2:

		a = 6378137.0;  b = 6356752.3142;     // WGS 84

                break;
	
	}

// Sobre la geometria del elipsoide
// excentricidad = e1; segunda excentricidad = e2

	double e1, e2;

	e1 = sqrt (pow(a,2) - pow(b, 2)) / a;
	e2 = sqrt (pow(a,2) - pow(b, 2)) / b;

// Radio polar de curvatura y aplanamiento
// radio polar de curvatura = c; aplanamiento = alpha

	double c, alpha;

	c = (a*a)/b;

	alpha = (a -b) / a;

	double long_gd_rad, lat_gd_rad;

	if (long_ == 'E' || long_ == 'e') {	

		long_gd_rad = -(v.a[0] * pi)/180;

	}

		else {

		long_gd_rad = (v.a[0] * pi)/180;
		
		}

	lat_gd_rad = (v.a[1] * pi) / 180;

// Determinacion del Huso

	double huso_dec;
	int huso;

	huso_dec = v.a[0]/6 + 31;

	huso = int (huso_dec);

// Obtencion del meridiano central del uso = lambda0

	double lambda0, delta_lambda;

	lambda0 = (huso * 6.0 - 183.0)*(pi/180.0);

// Determinacion de la distancia angular que existe entre la longitud del punto (long_gd_rad) y
// el meridiano central del huso (lamda0)

	delta_lambda = long_gd_rad - lambda0;

// Ecuaciones de Coticchia-Surace para el Problema Directo (Paso de Geograficas a UTM)
// Calculo de Parametros

	double A, xi, eta, nu, zeta, A1, A2, J2, J4, J6, alpha2, beta, gamma, B_phi;

	A = cos(lat_gd_rad) * sin(delta_lambda);

	xi = 0.5 * log((1+A)/(1-A));

	eta = atan(tan(lat_gd_rad)/cos(delta_lambda)) - lat_gd_rad;

	nu = (c*0.9996)/sqrt((1 + e2*e2*cos(lat_gd_rad)*cos(lat_gd_rad)));

	zeta = (e2*e2/2)*(xi*xi)*(cos(lat_gd_rad)*cos(lat_gd_rad));

	A1 = sin(2.0*lat_gd_rad);

	A2 = A1 * (cos(lat_gd_rad)*cos(lat_gd_rad));

	J2 = lat_gd_rad + A1/2.0;

	J4 = (3*J2 + A2)/4;

	J6 = (5*J4 + A2 * (cos(lat_gd_rad)*cos(lat_gd_rad)))/3;

	alpha2 = (3.0/4.0)*(e2*e2);

	beta = (5.0/3.0)*(alpha2*alpha2);

	gamma = (35.0/27.0)*(pow(alpha2,3));

	B_phi = 0.9996 * c * (lat_gd_rad - alpha2 * J2 + beta * J4 - gamma * J6);

	t.a[0] = xi*nu*(1+zeta/3.0)+500000.00;

	t.a[1] = eta*nu*(1+zeta)+B_phi;

	if (lat_ == 'S' || lat_ == 's')

		t.a[1] += 10000000.0;
	
	t.a[2] = v.a[2];

}

Implementado en el programa principal para el vértice de Carbonera (España); cuyas coordenadas geográficas decimales son (-3.598075790068, 39.547561311850, 697.897194861660) y el datum es ED50, es decir, opción 1 (elipsoide de salida Hayford).

#include "Coordenadas.h"
#include <gdal/ogr_geometry.h>

int main(){

	system("clear");

	Coordenadas a, b; // Declara objetos de la Clase Coordenadas 

	double u[3];
	
	u[0]=-3.598075790068; u[1]=39.547561311850; u[2]=697.897194861660;

	a.Iniciar(u); // Inicia el objeto a (vertice de Carbonera en coordenadas geograficas)

	cout << "vertice de Carbonera en Coordenadas geograficas:\n" << a;

	b.CambioGeoUtm(1,'o', 'n', a, b);        //De geograficas a UTM (cambio de datum) 

	cout << "vertice de Carbonera en Coordenadas UTM (datum ED50):\n" << b;

	return 0;

}

La salida es:

vertice de Carbonera en Coordenadas geograficas:
(-3.598075790068, 39.547561311850, 697.897194861660)
vertice de Carbonera en Coordenadas UTM (datum ED50):
(448610.601283641299, 4377788.156761185266, 697.897194861660)

que concuerda con lo esperado.

About these ads
Esta entrada fue publicada en Código C++, SIG, Transformar Coordenadas. Guarda el enlace permanente.

Deja un comentario

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s