Convertir coordenadas Geograficas a UTM con Tcl

Buenas y santas,

siguiendo los pasos descritos en el excelente trabajo publicado en la página GabrielOrtiz.com, he realizado un script TCL que realiza la conversión de Coordenadas Geométricas a UTM. Para la transformación se ha utilizado el elipsoide internacional o de Hayford, que es el que toma como referencia el DATUM ED50, utilizado en la mayor parte de la cartografía española.


proc GEOtoUTM {lat lon} {

# Elipsoide de Hayford
set Semieje_a 6378388.0
set Semieje_b 6356911.946130

# Constantes
set Pi 3.1415926535897932384626433832795

##########################
# CALCULOS

# Calculo de la excentricidad

set e [expr (sqrt(pow($Semieje_a,2)-pow($Semieje_b,2)))/$Semieje_a]

# Calculo de la segunda excentricidad y su cuadrado
set e_segunda [expr (sqrt(pow($Semieje_a,2)-pow($Semieje_b,2)))/$Semieje_b]
set e_segunda_cuadrado [expr pow($e_segunda,2)]

# Calculo del radio polar de curvatura
set r_pol_curv [expr pow($Semieje_a,2)/$Semieje_b]

# Calculo del aplanamiento
set aplanamiento [expr ($Semieje_a-$Semieje_b)/$Semieje_a]

set lat_rad [expr $lat*$Pi/180]
set lon_rad [expr $lon*$Pi/180]

set huso [expr entier(($lon_rad/6)+31)]

set merid_central [expr ($huso * 6) – 183]
set dist_ang [expr $lon_rad – (-3 * $Pi / 180)]

##########################################

set A [expr cos($lat_rad) * sin($dist_ang)]

set epsilon [expr (1.0/2.0)*(log((1.0+$A)/(1.0-$A)))]
set nu [expr atan(tan($lat_rad)/cos($dist_ang))-$lat_rad]
set v [expr ($r_pol_curv/sqrt(1+$e_segunda_cuadrado*pow(cos($lat_rad),2))*0.9996)]
set tao [expr ($e_segunda_cuadrado/2)*pow($epsilon,2)*pow(cos($lat_rad),2)]

set A1 [expr sin(2*$lat_rad)]
set A2 [expr $A1*pow(cos($lat_rad),2)]

set J2 [expr $lat_rad+($A1/2.0)]
set J4 [expr (3.0*$J2 + $A2)/4.0]
set J6 [expr (5.0*$J4 + $A2 *pow(cos($lat_rad),2))/3.0]

set alfa [expr (3.0/4.0)*$e_segunda_cuadrado]
set beta [expr (5.0/3.0)*pow($alfa,2)]
set gamma [expr (35.0/27.0)*pow($alfa,3)]
set B_phi [expr 0.9996*$r_pol_curv*($lat_rad-$alfa*$J2+$beta*$J4-$gamma*$J6)]

set X [expr $epsilon*$v*(1.0+($tao/3.0))+500000]
set Y [expr $nu*$v*(1.0+$tao)+$B_phi]

#############################################
# Imprimir resultados intermediarios
# para comprobar

# Seleccionar imprimir o no
set print_res FALSE

if {$print_res ==”TRUE”} {

puts $e
puts $e_segunda
puts $e_segunda_cuadrado
puts “$r_pol_curv $r_pol_curv”
puts $aplanamiento
puts “Lat_Rad $lat_rad”
puts $huso
puts $merid_central
puts $dist_ang
puts $A
puts $epsilon
puts $nu
puts $v
puts $tao
puts “$A1 $A2”
puts “$J2 $J4 $J6”
puts “$alfa $beta $gamma”
puts “$B_phi = 0.9996 . $r_pol_curv . ($lat_rad – $alfa . $J2 + $beta . $J4 – $gamma . $J6)”

}

# Imprimir resultado en la consola de salida
puts “X = [format %.3f $X] || Y = [format %.3f $Y]”

}

GEOtoUTM 43.4884075 -3.801873306

Para probarlo podéis consultar la guía para utilizar Textpad para programar en Tcl.

Un saludo

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

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