Javier Abadia Miranda (abadia++at++sgonyx.ita.es)
Thu, 30 Oct 1997 09:46:48 +0100
I give it for free, but not for nothing.
PLEASE, if you use this code, let me know.
This is the only thing I ask for.
Any comments welcome.
Aaadios. (Bye)
/*
* Javier Abadia Miranda - Oct 1997
* Ingeniero en Informatica (casi)
*
* Send comments to abadia++at++pandora.ita.es
*
* Centro Politecnico Superior - Instituto Tecnologico de Aragon
* Zaragoza (Spain)
*
* utm <-> geodetic coordinates conversion
*
* Using ecuations published in Boletino di Geodesia #1 (Italy)
* by Alberto Coticchia and Luciano Surace
*
* This is a good aproximation to exact expressions.
*
*/
/*
* -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
* Please, I *want* to know that you are using this code
* An e-mail stating "I'm X and I am using your utm2geo
* routines to ..." is ok. Thank you.
* Remember mailto://abadia++at++pandora.ita.es
* -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
*/
#include <stdio.h>
#include <math.h>
#ifdef MEGAPRECISION
#define REAL long double
#define cos(a) cosl(a)
#define sin(a) sinl(a)
#define tan(a) tanl(a)
#define atan(a) atanl(a)
#define atan2(a,b) atan2l(a,b)
#define log(a) logl(a)
#define exp(a) expl(a)
#define pow(a,b) powl(a,b)
#else
#define REAL double
#endif
/*
* For more info about UTM, zones or coordinates follow this link
*
* http://www.utexas.edu/ftp/depts/grg/gcraft/notes/mapproj/mapproj.html
*/
typedef struct {
REAL e2; /* excentricity squared */
REAL c; /* meters */
REAL a; /* meters */
char *datum; /* name of ellipsoid's datum */
} Elipsoide;
Elipsoide hayford = {
0.0067681703, /* e2 */
6399936.609, /* c */
6378388.0, /* a */
"Postdam"
};
/*
* latitude: measured in degrees, with origin at ecuator
* latitude > 0, North
* latitude < 0, South
*/
/*
* longitude: measured in degrees, with origin at Greenwich
* meridian
* longitude > 0, East
* longitude < 0, West
*/
/*
* X,Y: measured in metres; huso = UTM zone
*/
void utm2geo(int huso, REAL X, REAL Y, REAL *lat, REAL *lon);
void geo2utm(REAL lat, REAL lon, int *huso, REAL *X, REAL *Y);
/*
* returns a pointer to a STATIC buffer where there
* is a representation of g in [gg§mm'ss"mmm N] printable
* form. NEVER use this function twice in the same printf.
* Second parameter indicates if g is a latitude (true)
* or a longitude (false)
*/
char *grados_a_cadena(REAL g, int latitud);
/*
* main() only to test
*/
void main()
{
int huso;
REAL X,Y;
REAL latitud, longitud;
printf("Conversion UTM a Geograficas ");
printf("usando %d bits de precision\n", sizeof(REAL) * 8);
/*
* entrada de datos
*/
printf("huso.....: "); scanf("%d", &huso);
printf("x........: "); scanf("%f", &X);
printf("y........: "); scanf("%f", &Y);
printf("huso: %d, ",huso);
printf("X: %.3f, ",X);
printf("Y: %.3f\n",Y);
utm2geo(huso, X,Y, &latitud, &longitud);
printf("lat(phi) %s, ", grados_a_cadena(latitud, 1));
printf("lon(lambda) %s\n", grados_a_cadena(longitud, 0));
geo2utm(latitud, longitud, &huso, &X, &Y);
printf("huso: %d, ",huso);
printf("X: %.3f, ",X);
printf("Y: %.3f\n",Y);
}
/*
* STATIC functions
*/
/*
* given a UTM zone, calculates longitude (Easting)
* of its central meridian
*/
static REAL meridiano_central_del_huso(int huso);
/*
* calculates which zone is the place in,
* given its Easting
*/
static int determinar_huso(REAL lambda);
/*
* intermediate calculation
*/
static REAL calcular_B_sub_phi(REAL phi)
{
REAL B;
REAL A1,A2,J2,J4,J6,alpha,beta,gamma;
REAL cos2_phi = pow(cos(phi), 2.0);
A1 = sin( 2.0 * phi);
A2 = A1 * cos2_phi;
J2 = phi + A1 / 2.0;
J4 = (3.0 * J2 + A2) / 4.0;
J6 = (5.0 * J4 + A2 * cos2_phi) / 3.0;
alpha = 3.0/4.0 * hayford.e2;
beta = 5.0/3.0 * pow(alpha, 2.0);
gamma = 35.0/27.0 * pow(alpha, 3.0);
B = (0.9996 * hayford.c * (phi - alpha * J2 + beta * J4 - gamma * J6));
return B;
}
/*
* Caution: Never call this function twice inside a
* printf, unless you want to check for your compiler's
* parameter order of evaluation.
* Second call will destroy result of previous one
*/
char *grados_a_cadena(REAL g, int latitud)
{
static char buffer[100];
int grados, minutos, segundos, milesimas;
char c;
if(latitud)
c = (g >= 0) ? 'N' : 'S';
else
c = (g >= 0) ? 'E' : 'W';
g = fabs(g);
grados = (int) g;
g = (g - grados) * 60.0;
minutos = (int) g;
g = (g - minutos) * 60.0;
segundos = (int) g;
g = (g - segundos) * 1000.0;
milesimas = (int) g;
sprintf(buffer, "%2d§ %2d' %2d\"%03d %c", grados, minutos, segundos, milesimas, c);
return buffer;
}
void utm2geo(int huso, REAL X, REAL Y, REAL *lat, REAL *lon)
{
REAL y,x,phi,nu,a,b,tseta,xi,eta,senh_xi,delta_lambda,t;
REAL cos2_phi;
REAL B_sub_phi;
REAL sumando_1, sumando_2, corchete;
REAL latitud, longitud;
REAL lambda_0;
y = Y;
x = X - 500000.0 ;
phi = y / 6366197.724 / 0.9996 ;
cos2_phi = pow(cos(phi), 2.0);
nu = hayford.c * 0.9996 / pow((1.0 + hayford.e2 * cos2_phi), 0.5);
a = x / nu;
B_sub_phi = calcular_B_sub_phi(phi);
b = (y - B_sub_phi) / nu;
tseta = hayford.e2 * a * a / 2.0 * cos2_phi;
xi = a * (1.0 - tseta / 3.0);
eta = b * (1.0 - tseta) + phi;
senh_xi = (exp(xi) - exp(-xi)) / 2.0;
delta_lambda = atan2( senh_xi, cos(eta));
t = atan(cos(delta_lambda) * tan(eta));
sumando_1 = hayford.e2 * cos2_phi;
sumando_2 = 3.0 / 2.0 * hayford.e2 * sin(phi) * cos(phi) * (t - phi);
corchete = 1.0 + sumando_1 - sumando_2;
lambda_0 = meridiano_central_del_huso(huso) * M_PI / 180.0;
latitud = phi + corchete * (t-phi); /* phi */
longitud = delta_lambda + lambda_0; /* lambda */
*lat = latitud * 180.0 / M_PI;
*lon = longitud * 180.0 / M_PI;
}
void geo2utm(REAL lat, REAL lon, int *huso, REAL *X, REAL *Y)
{
REAL xi, A, n, nu, tseta, B_sub_phi, delta_lambda;
REAL lambda_0, lambda, phi, cos2_phi;
phi = lat * M_PI / 180.0;
lambda = lon * M_PI / 180.0;
*huso = determinar_huso(lon);
lambda_0 = meridiano_central_del_huso(*huso) * M_PI / 180.0;
delta_lambda = lambda - lambda_0;
A = cos(phi) * sin(delta_lambda);
xi = 1.0 / 2.0 * log( (1.0+A) / (1.0-A) );
n = atan2( tan(phi), cos(delta_lambda) ) - phi;
cos2_phi = pow(cos(phi), 2.0);
nu = hayford.c * 0.9996 / pow( (1.0 + hayford.e2 * cos2_phi), 0.5);
tseta = hayford.e2 / 2.0 * xi * xi * cos2_phi;
B_sub_phi = calcular_B_sub_phi(phi);
*Y = n * nu * (1.0 + tseta) + B_sub_phi;
*X = xi * nu * (1.0 + tseta/3.0) + 500000.0;
}
/*
* UTM zones calculations (huso = UTM zone)
*/
static REAL meridiano_central_del_huso(int huso)
{
REAL lambda;
lambda = 6.0 * huso - 183.0;
return lambda;
}
static int determinar_huso(REAL lambda)
{
int huso;
huso = (lambda / 6.0) + 31;
while(huso < 1)
huso += 60;
return huso;
}
===================================List Archives, FAQ, FTP: http://www.sgi.com/Technology/Performer/
Submissions: info-performer++at++sgi.com
Admin. requests: info-performer-request++at++sgi.com
This archive was generated by hypermail 2.0b2 on Mon Aug 10 1998 - 17:56:08 PDT