Python Online: Projeção TRANSVERSA de Mercator
- Adauto Costa
- 26 de abr. de 2023
- 7 min de leitura
Acesse o https://colab.research.google.com/, crie um novo "notebook" e teste o código abaixo:
import math
from math import sin
from math import cos
from math import sqrt
from math import tan
from math import atan
class Elipsoide:
def __init__(self, nome, a, f):
fi = 1 / f
self.a = a
self.b = a - a / f
self.e = sqrt(2 * fi - fi ** 2)
self.e2 = sqrt((a ** 2 - self.b ** 2) / self.b ** 2)
self.Nome = nome
def GrandeNormal(self, lat):
e = self.e
N = self.a / sqrt(1 - e**2 * sin(lat)**2)
return N
def RaioCurvaturaSecaoMeridiana(self, lat):
return self.a * (1 - self.e ** 2) / (1 - (self.e * sin(lat)) ** 2) ** (3/2)
def ArcoDeElipseMeridiana(self, lat):
a = self.a
e = self.e
A = 1 + (3 / 4) * e ** 2 + (45 / 64) * e ** 4 + (175 / 256) * e ** 6 + (11025 / 16384) * e ** 8 + (43659 / 65536) * e ** 10
B = (3 / 4) * e ** 2 + (15 / 16) * e ** 4 + (525 / 512) * e ** 6 + (2205 / 2048) * e ** 8 + (72765 / 65536) * e ** 10
C = (15 / 64) * e ** 4 + (105 / 256) * e ** 6 + (2205 / 4096) * e ** 8 + (10395 / 16384) * e ** 10
D = (35 / 512) * e ** 6 + (315 / 2048) * e ** 8 + (31185 / 131072) * e ** 10
E = (315 / 16384) * e ** 8 + (3465 / 65536) * e ** 10
F = (693 / 131072) * e ** 10
o2 = lat
o1 = 0
s0 = o2 - o1
s1 = sin(2 * o2) - sin(2 * o1)
s2 = sin(4 * o2) - sin(4 * o1)
s3 = sin(6 * o2) - sin(6 * o1)
s4 = sin(8 * o2) - sin(8 * o1)
s5 = sin(10 * o2) - sin(10 * o1)
S = a * (1 - e ** 2) * (A * s0 - 1 / 2 * B * s1 + 1 / 4 * C * s2 - 1 / 6 * D * s3 + 1 / 8 * E * s4 - 1 / 10 * F * s5)
return S
class TransversaMercator:
def __init__(self, el, lat, lon, m=0):
lat *= math.pi / 180 #latitude convertida para radianos
lon *= math.pi / 180 #longitude convertida para radianos
m *= math.pi / 180 #meridiano central convertido para radianos
self.Elip = el #objeto Elipsoide
self.m = m
self.dl = lon - m
self.S = el.ArcoDeElipseMeridiana(lat)
self.N = el.GrandeNormal(lat) #cálculo da Grande Normal
e2 = el.e2 #segunda excentricidade do elipsóide
t = tan(lat)
na = e2 * cos(lat)
self.x = self.N * (self.dl * cos(lat) +
1 / 6 * self.dl ** 3 * cos(lat) ** 3 * (1 - t ** 2 + na ** 2) +
1 / 120 * self.dl ** 5 * cos(lat) ** 5 * (5 - 18 * t ** 2 + t ** 4 + 14 * na ** 2 - 58 * t ** 2 * na ** 2 + 13 * na ** 4 + 4 * na ** 6 - 64 * na ** 4 * t ** 2 - 24 * na ** 6 * t ** 2) +
1 / 5040 * self.dl ** 7 * cos(lat) ** 7 * (61 - 479 * t ** 2 + 179 * t ** 4 - t ** 6))
self.y = self.N * (self.S / self.N +
1 / 2 * self.dl ** 2 * sin(lat) * cos(lat) +
1 / 24 * self.dl ** 4 * sin(lat) * cos(lat) ** 3 * (5 - t ** 2 + 9 * na ** 2 + 4 * na ** 4) +
1 / 720 * self.dl ** 6 * sin(lat) * cos(lat) ** 5 * (61 - 58 * t ** 2 + t ** 4 + 270 * na ** 2 - 330 * t ** 2 * na ** 2 + 445 * na ** 4 + 324 * na ** 6 - 680 * na ** 4 * t ** 2 + 88 * na ** 8 - 600 * na ** 6 * t ** 2 - 192 * na ** 8 * t ** 2) +
1 / 40320 * self.dl ** 8 * sin(lat) * cos(lat) ** 7 * (1385 - 311 * t ** 2 + 543 * t ** 4 - t ** 6))
c2 = (1 / math.sin(1 / 3600) * math.pi/180) ** 2
self.c = self.dl * sin(lat) * (1 + (1 / (3 * c2) * (self.dl * cos(lat)) ** 2) * (1 + 3 * na ** 2 + 2 * na ** 4) + (1 / (15 * c2 ** 2) * (self.dl * cos(lat)) ** 4) * (2 - t ** 2))
self.c *= 180/math.pi #convergência meridiana convertida para graus decimais
self.c = round(self.c, 4)
self.k = 1 + 1 / 2 * self.dl ** 2 * cos(lat) ** 2 * (1 + na ** 2) + 1 / 24 * self.dl ** 4 * cos(lat) ** 4 * (5 - 4 * t ** 2 + 14 * na ** 2 + 13 * na ** 4 - 28 * t ** 2 * na ** 2 + 4 * na ** 6 - 48 * t ** 2 * na ** 4 - 24 * t ** 2 * na ** 6) + 1 / 720 * self.dl ** 6 * cos(lat) ** 6 * (61 - 148 * t ** 2 + 16 * t ** 4)
self.k = round(self.k, 4)
def ProblemaInverso(el, x, y, m):
a = el.a #semieixo maior do elipsóide
e = el.e #primeira excentricidade
e2 = el.e2 #segunda excentricidade
lat = [0] * 50
lat[0] = y / a
A0 = 1 - 1/4 * e ** 2 - 3/64 * e ** 4 - 5/256 * e ** 6 - 175/16384 * e ** 8
A2 = 3/8 * (e ** 2 + 1/4 * e ** 4 + 15/128 * e ** 6 - 455/4096 * e ** 8)
A4 = 15/256 * (e ** 4 + 3/4 * e ** 6 - 77/128 * e ** 8)
A6 = 35/3072 * (e ** 6 - 41/32 * e ** 8)
A8 = -315/131072 * e ** 8
i = 0
for i in range(1, 51):
f = a * (A0 * lat[i - 1] - A2 * sin(2 * lat[i - 1]) + A4 * sin(4 * lat[i - 1]) - A6 * sin(6 * lat[i - 1]) + A8 * sin(8 * lat[i - 1])) - y
fl = a * (A0 - 2 * A2 * cos(2 * lat[i - 1]) + 4 * A4 * cos(4 * lat[i - 1]) - 6 * A6 * cos(6 * lat[i - 1]) + 8 * A8 * cos(8 * lat[i - 1]))
lat[i] = lat[i - 1] - f / fl
if abs(lat[i] - lat[i - 1]) < 0.000000000001:
break
lat1 = lat[i]
t1 = tan(lat1)
na1 = e2 * cos(lat1)
N1 = el.GrandeNormal(lat1)
M1 = el.RaioCurvaturaSecaoMeridiana(lat1)
lat2 = lat1 - (t1 * x ** 2 / (2 * M1 * N1)) + (t1 * x ** 4 / (24 * M1 * N1 ** 3)) * (5 + 3 * t1 ** 2 + na1 ** 2 - 4 * na1 ** 4 - 9 * na1 ** 2 * t1 ** 2) - (t1 * x ** 6 / (720 * M1 * N1 ** 5)) * (61 - 90 * t1 ** 2 + 46 * na1 ** 2 + 45 * t1 ** 4 - 252 * t1 ** 2 * na1 ** 2 - 3 * na1 ** 4 + 100 * na1 ** 6 - 66 * t1 ** 2 * na1 ** 4 - 90 * t1 ** 4 * na1 ** 2 + 88 * na1 ** 8 + 225 * t1 ** 4 * na1 ** 4 + 84 * t1 ** 2 * na1 ** 6 - 192 * t1 ** 2 * na1 ** 8) + (t1 * x ** 8 / (40320 * M1 * N1 ** 7)) * (1385 + 3633 * t1 ** 2 + 4095 * t1 ** 4 + 1575 * t1 ** 6)
lon2 = 1 / cos(lat1) * (x / N1 -
1 / 6 * (x / N1) ** 3 * (1 + 2 * t1 ** 2 + na1 ** 2) +
1 / 120 * (x / N1) ** 5 * (5 + 6 * na1 ** 2 + 28 * t1 ** 2 - 3 * na1 ** 4 + 8 * t1 ** 2 * na1 ** 2 + 24 * t1 ** 4 - 4 * na1 ** 6 + 4 * t1 ** 2 * na1 ** 4 + 24 * t1 ** 2 * na1 ** 6) -
1 / 5040 * (x / N1) ** 7 * (61 + 662 * t1 ** 2 + 1320 * t1 ** 4 + 720 * t1 ** 6))
c = atan(t1 / N1 * x -
t1 / 3 * (x / N1) ** 3 * (1 - na1 ** 2 - 2 * na1 ** 4) +
t1 / 15 * (x / N1) ** 5 * (2 + 2 * na1 ** 2 + 9 * na1 ** 4 + 6 * t1 ** 2 * na1 ** 2 + 20 * na1 ** 6 + 3 * t1 ** 2 * na1 ** 4 - 27 * t1 ** 2 * na1 ** 6 + 11 * na1 ** 8 - 24 * t1 ** 2 * na1 ** 8) -
17 * t1 / 315 * (x / N1) ** 7)
k = 1 + (1 + na1 ** 2) / 2 * (x / N1) ** 2 + 1 / 24 * (1 + 6 * na1 ** 2 + 9 * na1 ** 4 + 4 * na1 ** 6 - 24 * t1 ** 2 * na1 ** 4 - 24 * t1 ** 2 * na1 ** 6) * (x / N1) ** 4 + 1 / 720 * (x / N1) ** 6
return [round(lat2 * 180 / math.pi,8), round(lon2 * 180 / math.pi + m,8), round(c * 180 / math.pi,4), round(k, 4)] #ângulos convertidos para graus decimais
lat = -7.5 #Latitude em graus decimais
lon = -34.2 #Longitude em graus decimais
m = -33.0 #Meridiano Central em graus decimais
el = Elipsoide("GRS80", 6378137, 298.257222101)
tmer = TransversaMercator(el, lat, lon, m)
tmer_inv = TransversaMercator.ProblemaInverso(el, tmer.x, tmer.y, m)
str = f"Coord. de Entrada (lat, lon, merid central):\n({lat},{lon},{m})"
print(str)
print("Problema Direto:")
str = f"(x, y, conv. meridiana, fator escala):\n({round(tmer.x,3)},{round(tmer.y,3)},{tmer.c},{tmer.k})"
print(str)
print("Problema Inverso:")
str = f"(lat, lon, conv. meridiana, fator escala):\n({tmer_inv[0]},{tmer_inv[1]},{tmer_inv[2]},{tmer_inv[3]})"
print(str)Para acessar o formulário matemático da Projeção Transversa de Mercator clique aqui.
No código acima, altere as coordenadas de entrada e insira as suas. Neste exemplo as coordenadas de entrada foram latitude -7.5°, longitude -34.2° e Meridiano Central -33º. Considerada latitude de origem igual a zero.
Comments