Bom, hoje vou voltar um pouco às origens do blog!!!
Um pouco de autolisp pra relembrar os velhos tempos de programação em POG!!!
A lisp abaixo na verdade são algumas subrotinas para conversão de coordenadas geográficas em UTM e viceversa.
É bem fácil de usar se você souber o que é UTM e coordenada geográfica e está familiarizado com Georeferenciamento.
Muitas pessoas me perguntam se eu tenho uma rotina pra converter e.... Bem, tenho!!
Está aí!!
;| Conversão de UTM para GEOGRAFICA
baseado em http://recursos.gabrielortiz.com/index.asp?Info=058a
metodo: Coticchia-Surace
ay: coorenada do semi eixo maior (m)
bx: coorenada do semi eixo menor (m)
pt: (cord_X coord_Y coord_Z)
fuso: fuso, inteiro
hmsf: hemisferio, "N" para norte e "S" para sul
se for conhecido f, temos: bx=(1-f)*ay
sad69 -> ay = 6.378.160,000m e f = 1/298,25
corrego alegre -> ay = 6.378.388,000m e f = 1/297,00
(utm2geo (GETPOINT) 6378160.0 298.25 22 "S") => 25º25'51"15216 49º17'02"51881
|;
(defun utm2geo (pt ay f fuso hmsf / tmp e el el² c alpha lo phil nu A1 A2 J2 J4 re
J6 beta gama Bo zeta xi eta sinhxi dl tau lat lon x y x1 y1 fe bx)
(setq ay (float ay)
bx (* ay (- 1 (/ 1.0 f)))
x (car pt)
y (cadr pt)
re 6366197.724 ;raio da terra
fe 0.9996 ;fator de escala
tmp (sqrt (- (expt ay 2) (expt bx 2)))
e (/ tmp ay) ;excentricidade
el (/ tmp bx) ;2ª excentricidade
el² (expt el 2)
c (/ (expt ay 2) bx);raio polar de curvatura
x1 (- x 500000.0)
y1 (if (or (= hmsf 'N) (= hmsf "N")) y (- y 10000000.0))
phil (/ y1 (* re fe))
lo (- (* 6 fuso) 183)
nu (/ (* c fe) (sqrt (1+ (* el² (expt (cos phil) 2)))))
a (/ x1 nu)
A1 (sin (* 2 phil))
A2 (* A1 (expt (cos phil) 2))
J2 (+ phil (/ A1 2.0))
J4 (/ (+ (* J2 3.0) A2) 4.0)
J6 (/ (+ (* 5.0 J4) (* A2 (expt (cos phil) 2))) 3.0)
alpha (/ (* 3.0 el²) 4.0)
beta (* (/ 5.0 3.0) (expt alpha 2))
gama (* (/ 35.0 27.0) (expt alpha 3))
Bo (* fe c (+ phil (* (- alpha) J2) (* beta J4) (* (- gama) J6)))
b (/ (- y1 Bo) nu)
zeta (* (/ (* el² (expt a 2)) 2.0) (expt (cos phil) 2))
xi (* a (- 1.0 (/ zeta 3.0)))
eta (+ (* b (- 1.0 zeta)) phil)
sinhxi (/ (- (exp xi) (exp (- xi))) 2.0)
dl (atan (/ sinhxi (cos eta)))
tau (atan (* (cos dl) (tan eta)))
lon (+ (* (/ 180.0 pi) dl) lo)
lat (* (/ 180.0 pi)
(+ phil (* (+ 1.0
(* el² (expt (cos phil) 2.0))
(* (/ -3.0 2.0) el² (sin phil) (cos phil) (- tau phil)))
(- tau phil)))))
(if (caddr pt)
(list lon lat (caddr pt))
(list lon lat 0.0)))
;pt -> long lat
;a -> semi eixo maior
;f -> achatamento
(defun geo2utm (pt a f / b e el el² c lamb fi fuso lo deltal Am eps n v
S A1 A2 J2 J4 J6 alfa beta gama bo)
(setq a (float a)
b (- a (/ a f))
el (/ (sqrt (- (expt a 2) (expt b 2))) b)
el² (expt el 2)
c (/ (expt a 2) b)
fuso (fix (+ (/ (car pt) 6.0) 31))
lamb (/ (* (car pt) pi) 180.0)
fi (/ (* (cadr pt) pi) 180.0)
lo (- (* fuso 6) 183) ;meridiano central
deltal (- lamb (/ (* lo pi) 180.0))
Am (* (cos fi) (sin deltal))
eps (* 0.5 (log (/ (+ 1 Am) (- 1 Am))))
n (- (atan (/ (tan fi) (cos deltal))) fi)
v (/ (* c 0.9996) (sqrt (+ 1 (* el² (expt (cos fi) 2)))))
S (/ (expt (* el eps (cos fi)) 2) 2.0)
A1 (sin (* 2.0 fi))
A2 (* A1 (expt (cos fi) 2.0))
J2 (+ fi (/ A1 2.0))
J4 (/ (+ (* 3.0 J2) A2) 4.0)
J6 (/ (+ (* 5 J4) (* A2 (expt (cos fi) 2))) 3.0)
alfa (/ (* 3.0 el²) 4.0)
beta (* (/ 5.0 3.0) (expt alfa 2))
gama (* (/ 35.0 27.0) (expt alfa 3))
bo (* 0.9996 c (+ fi (* (- alfa) J2) (* beta J4) (* (- gama) J6))))
(list (+ 500000.0 (* eps v (1+ (/ S 3.0)))) ;x
(+ bo (* n v (1+ S)) (if (< lat 0.0) 10000000.0 0.0));y
(caddr pt)
))
(defun LLA_wgs84->sad69 (pt)
(geo2geo pt 6378137.0 298.257223563 6378160.0 298.25 66.87 -4.37 38.52))
(defun LLA_sad69->wgs84 (pt)
(geo2geo pt 6378160.0 298.25 6378137.0 298.257223563 -66.87 4.37 -38.52))
(defun geo2geo (pt ;long_from lat_from h_from ;coordenadas geodesicas de origem
a_from f_from ;parametros geodesicos de origem
a_to f_to ;parametros geodesicos de destino
dx dy dz ;translação origem->destino
/ lat1 long1 f1 f2 a1 a2 e²1 e²2 N1 N2 Xw Yw Zw b2 p teta fi lamb hb ep²)
(setq lat1 (/ (* (cadr pt) pi) 180.0)
long1 (/ (* (car pt) pi) 180.0)
f1 (/ 1.0 f_from) ;wgs
a1 a_from ;wgs
e²1 (* f1 (- 2.0 f1))
N1 (/ a1 (sqrt (- 1.0 (* e²1 (expt (sin lat1) 2.0)))))
;coord carteziana no sistema de origem:
Xw (* (+ N1 (caddr pt)) (cos lat1) (cos long1))
Yw (* (+ N1 (caddr pt)) (cos lat1) (sin long1))
Zw (* (+ (* N1 (- 1.0 e²1)) (caddr pt)) (sin lat1))
;coord carteziana do ponto no novo sistema:
Xb (+ Xw dx)
Yb (+ Yw dy)
Zb (+ Zw dz)
;converter carteziana para geodesica:
f2 (/ 1.0 f_to)
a2 a_to
e²2 (* f2 (- 2.0 f2))
b2 (* a2 (- 1.0 f2))
p (sqrt (+ (expt Xb 2) (expt Yb 2)))
teta (atan (/ (* Zb a2) (* p b2)))
ep² (/ (- (expt a2 2.0) (expt b2 2.0)) (expt b2 2.0))
fi (atan (/ (+ Zb (* ep² b2 (expt (sin teta) 3.0))) (- p (* e²2 a2 (expt (cos teta) 3.0)))))
lamb (atan (/ Yb Xb))
N2 (/ a2 (sqrt (- 1.0 (* e²2 (expt (sin fi) 2.0)))))
hb (- (/ p (cos fi)) N2))
(list (/ (* 180.0 lamb) pi);long
(/ (* 180.0 fi) pi) ;lat
hb) ;altitude
)
Link(s) da(s) subrotina(s) usada(s): tan
A dificuldade nem está nos cálculos em si, mas na sintaxe das fórmulas, não acham??
;| Conversão de UTM para GEOGRAFICA
baseado em http://recursos.gabrielortiz.com/index.asp?Info=058a
metodo: Coticchia-Surace
ay: coorenada do semi eixo maior (m)
bx: coorenada do semi eixo menor (m)
pt: (cord_X coord_Y coord_Z)
fuso: fuso, inteiro
hmsf: hemisferio, "N" para norte e "S" para sul
se for conhecido f, temos: bx=(1-f)*ay
sad69 -> ay = 6.378.160,000m e f = 1/298,25
corrego alegre -> ay = 6.378.388,000m e f = 1/297,00
(utm2geo (GETPOINT) 6378160.0 298.25 22 "S") => 25º25'51"15216 49º17'02"51881
|;
(defun utm2geo (pt ay f fuso hmsf / tmp e el el² c alpha lo phil nu A1 A2 J2 J4 re
J6 beta gama Bo zeta xi eta sinhxi dl tau lat lon x y x1 y1 fe bx)
(setq ay (float ay)
bx (* ay (- 1 (/ 1.0 f)))
x (car pt)
y (cadr pt)
re 6366197.724 ;raio da terra
fe 0.9996 ;fator de escala
tmp (sqrt (- (expt ay 2) (expt bx 2)))
e (/ tmp ay) ;excentricidade
el (/ tmp bx) ;2ª excentricidade
el² (expt el 2)
c (/ (expt ay 2) bx);raio polar de curvatura
x1 (- x 500000.0)
y1 (if (or (= hmsf 'N) (= hmsf "N")) y (- y 10000000.0))
phil (/ y1 (* re fe))
lo (- (* 6 fuso) 183)
nu (/ (* c fe) (sqrt (1+ (* el² (expt (cos phil) 2)))))
a (/ x1 nu)
A1 (sin (* 2 phil))
A2 (* A1 (expt (cos phil) 2))
J2 (+ phil (/ A1 2.0))
J4 (/ (+ (* J2 3.0) A2) 4.0)
J6 (/ (+ (* 5.0 J4) (* A2 (expt (cos phil) 2))) 3.0)
alpha (/ (* 3.0 el²) 4.0)
beta (* (/ 5.0 3.0) (expt alpha 2))
gama (* (/ 35.0 27.0) (expt alpha 3))
Bo (* fe c (+ phil (* (- alpha) J2) (* beta J4) (* (- gama) J6)))
b (/ (- y1 Bo) nu)
zeta (* (/ (* el² (expt a 2)) 2.0) (expt (cos phil) 2))
xi (* a (- 1.0 (/ zeta 3.0)))
eta (+ (* b (- 1.0 zeta)) phil)
sinhxi (/ (- (exp xi) (exp (- xi))) 2.0)
dl (atan (/ sinhxi (cos eta)))
tau (atan (* (cos dl) (tan eta)))
lon (+ (* (/ 180.0 pi) dl) lo)
lat (* (/ 180.0 pi)
(+ phil (* (+ 1.0
(* el² (expt (cos phil) 2.0))
(* (/ -3.0 2.0) el² (sin phil) (cos phil) (- tau phil)))
(- tau phil)))))
(if (caddr pt)
(list lon lat (caddr pt))
(list lon lat 0.0)))
;pt -> long lat
;a -> semi eixo maior
;f -> achatamento
(defun geo2utm (pt a f / b e el el² c lamb fi fuso lo deltal Am eps n v
S A1 A2 J2 J4 J6 alfa beta gama bo)
(setq a (float a)
b (- a (/ a f))
el (/ (sqrt (- (expt a 2) (expt b 2))) b)
el² (expt el 2)
c (/ (expt a 2) b)
fuso (fix (+ (/ (car pt) 6.0) 31))
lamb (/ (* (car pt) pi) 180.0)
fi (/ (* (cadr pt) pi) 180.0)
lo (- (* fuso 6) 183) ;meridiano central
deltal (- lamb (/ (* lo pi) 180.0))
Am (* (cos fi) (sin deltal))
eps (* 0.5 (log (/ (+ 1 Am) (- 1 Am))))
n (- (atan (/ (tan fi) (cos deltal))) fi)
v (/ (* c 0.9996) (sqrt (+ 1 (* el² (expt (cos fi) 2)))))
S (/ (expt (* el eps (cos fi)) 2) 2.0)
A1 (sin (* 2.0 fi))
A2 (* A1 (expt (cos fi) 2.0))
J2 (+ fi (/ A1 2.0))
J4 (/ (+ (* 3.0 J2) A2) 4.0)
J6 (/ (+ (* 5 J4) (* A2 (expt (cos fi) 2))) 3.0)
alfa (/ (* 3.0 el²) 4.0)
beta (* (/ 5.0 3.0) (expt alfa 2))
gama (* (/ 35.0 27.0) (expt alfa 3))
bo (* 0.9996 c (+ fi (* (- alfa) J2) (* beta J4) (* (- gama) J6))))
(list (+ 500000.0 (* eps v (1+ (/ S 3.0)))) ;x
(+ bo (* n v (1+ S)) (if (< lat 0.0) 10000000.0 0.0));y
(caddr pt)
))
(defun LLA_wgs84->sad69 (pt)
(geo2geo pt 6378137.0 298.257223563 6378160.0 298.25 66.87 -4.37 38.52))
(defun LLA_sad69->wgs84 (pt)
(geo2geo pt 6378160.0 298.25 6378137.0 298.257223563 -66.87 4.37 -38.52))
(defun geo2geo (pt ;long_from lat_from h_from ;coordenadas geodesicas de origem
a_from f_from ;parametros geodesicos de origem
a_to f_to ;parametros geodesicos de destino
dx dy dz ;translação origem->destino
/ lat1 long1 f1 f2 a1 a2 e²1 e²2 N1 N2 Xw Yw Zw b2 p teta fi lamb hb ep²)
(setq lat1 (/ (* (cadr pt) pi) 180.0)
long1 (/ (* (car pt) pi) 180.0)
f1 (/ 1.0 f_from) ;wgs
a1 a_from ;wgs
e²1 (* f1 (- 2.0 f1))
N1 (/ a1 (sqrt (- 1.0 (* e²1 (expt (sin lat1) 2.0)))))
;coord carteziana no sistema de origem:
Xw (* (+ N1 (caddr pt)) (cos lat1) (cos long1))
Yw (* (+ N1 (caddr pt)) (cos lat1) (sin long1))
Zw (* (+ (* N1 (- 1.0 e²1)) (caddr pt)) (sin lat1))
;coord carteziana do ponto no novo sistema:
Xb (+ Xw dx)
Yb (+ Yw dy)
Zb (+ Zw dz)
;converter carteziana para geodesica:
f2 (/ 1.0 f_to)
a2 a_to
e²2 (* f2 (- 2.0 f2))
b2 (* a2 (- 1.0 f2))
p (sqrt (+ (expt Xb 2) (expt Yb 2)))
teta (atan (/ (* Zb a2) (* p b2)))
ep² (/ (- (expt a2 2.0) (expt b2 2.0)) (expt b2 2.0))
fi (atan (/ (+ Zb (* ep² b2 (expt (sin teta) 3.0))) (- p (* e²2 a2 (expt (cos teta) 3.0)))))
lamb (atan (/ Yb Xb))
N2 (/ a2 (sqrt (- 1.0 (* e²2 (expt (sin fi) 2.0)))))
hb (- (/ p (cos fi)) N2))
(list (/ (* 180.0 lamb) pi);long
(/ (* 180.0 fi) pi) ;lat
hb) ;altitude
)
Link(s) da(s) subrotina(s) usada(s): tan
A dificuldade nem está nos cálculos em si, mas na sintaxe das fórmulas, não acham??
