Skip to content

Instantly share code, notes, and snippets.

@george-silva
Last active December 18, 2015 00:59
Show Gist options
  • Select an option

  • Save george-silva/5700878 to your computer and use it in GitHub Desktop.

Select an option

Save george-silva/5700878 to your computer and use it in GitHub Desktop.
def azimute_reto(self,p1, p2):
print "reto"
delta_x = round(p2.x - p1.x, 2)
delta_y = round(p2.y - p1.y, 2)
print "dX", delta_x
print "dY", delta_y
if delta_x == 0 and delta_y > 0:
return 0
if delta_x > 0 and delta_y == 0:
return 90
if delta_x == 0 and delta_y < 0:
return 180
if delta_x < 0 and delta_y == 0:
return 270
return None
def azimute_puissant(self, p1, p2):
"""
Calcula azimute puissant
"""
reto = self.azimute_reto(p1, p2)
if reto:
return reto
f = 0.00335281068118
inv_f = 298.257222101
e2 = 0.00669438002290
a = 6378137
b = 6356752.314140
seno_1segundo = 0.00000484813681108
lat_media = (radians(p1.y) + radians(p2.y)) / 2
seno_lat_media = math.sin(lat_media)
cos_lat_media = math.cos(lat_media)
pow_seno_20 = math.pow(seno_lat_media, 2)
nm = a / (math.pow(1 - (e2 * pow_seno_20), 0.5))
delta_lat = (p2.y - p1.y) * 3600
delta_lon = (p2.x - p1.x) * 3600
mm = (a * (1 - e2)) / math.pow(1 - (e2 * pow_seno_20), 1.5)
bm = 1 / (mm * seno_1segundo)
x = delta_lon * cos_lat_media * nm * seno_1segundo
y = delta_lat * cos(radians(delta_lon / 7200)) * mm * seno_1segundo + 0.0000000001
F = (1 / 12) * seno_lat_media * cos_lat_media * cos_lat_media * seno_1segundo * seno_1segundo
gamma = (delta_lon * seno_lat_media * (1 / cos(radians(delta_lat / 7200))) + (
F * delta_lon * delta_lon * delta_lon))
if x < 0:
sinal_x = -1
elif x == 0:
sinal_x = 0
else:
sinal_x = 1
if y < 0:
sinal_y = -1
elif y == 0:
sinal_y = 0
else:
sinal_y = 1
if sinal_y == 0:
if delta_lat == 0 and delta_lon > 0:
return 90
# oeste
if delta_lat == 0 and delta_lon < 0:
return 270
if delta_lon > 0 and delta_lat > 0:
return math.degrees(math.atan(delta_lon / delta_lat))
if (delta_lon > 0 > delta_lat) or (delta_lon < 0 and delta_lat < 0):
return math.degrees(math.atan(delta_lon / delta_lat)) + 180
if delta_lon < 0 < delta_lat:
return math.degrees(math.atan(delta_lon / delta_lat)) + 360
azimute = 180 * (1 - (0.5 * sinal_x) - (0.5 * sinal_x * sinal_y)) + (math.degrees(math.atan(x / y)) - (
gamma / 7200))
return azimute
def azimute_reto(self,p1, p2, arredondar=True):
print "reto"
if arredondar:
delta_x = round(p2.x - p1.x, 3)
delta_y = round(p2.y - p1.y, 3)
else:
delta_x = p2.x - p1.x
delta_y = p2.y - p1.y
print "dX", delta_x
print "dY", delta_y
if delta_x == 0 and delta_y == 0:
return self.azimute_reto(p1, p2, arredondar=False)
if delta_x == 0 and delta_y > 0:
return 0
if delta_x > 0 and delta_y == 0:
return 90
if delta_x == 0 and delta_y < 0:
return 180
if delta_x < 0 and delta_y == 0:
return 270
return None
def azimute_puissant(self, p1, p2):
"""
Calcula azimute puissant
"""
reto = self.azimute_reto(p1, p2)
if reto:
return reto
f = 0.00335281068118
inv_f = 298.257222101
e2 = 0.00669438002290
a = 6378137
b = 6356752.314140
seno_1segundo = 0.00000484813681108
lat_media = (radians(p1.y) + radians(p2.y)) / 2
seno_lat_media = math.sin(lat_media)
cos_lat_media = math.cos(lat_media)
pow_seno_20 = math.pow(seno_lat_media, 2)
nm = a / (math.pow(1 - (e2 * pow_seno_20), 0.5))
delta_lat = (p2.y - p1.y) * 3600
delta_lon = (p2.x - p1.x) * 3600
mm = (a * (1 - e2)) / math.pow(1 - (e2 * pow_seno_20), 1.5)
bm = 1 / (mm * seno_1segundo)
x = delta_lon * cos_lat_media * nm * seno_1segundo
y = delta_lat * cos(radians(delta_lon / 7200)) * mm * seno_1segundo + 0.0000000001
F = (1 / 12) * seno_lat_media * cos_lat_media * cos_lat_media * seno_1segundo * seno_1segundo
gamma = (delta_lon * seno_lat_media * (1 / cos(radians(delta_lat / 7200))) + (
F * delta_lon * delta_lon * delta_lon))
if x < 0:
sinal_x = -1
elif x == 0:
sinal_x = 0
else:
sinal_x = 1
if y < 0:
sinal_y = -1
elif y == 0:
sinal_y = 0
else:
sinal_y = 1
if sinal_y == 0:
if delta_lat == 0 and delta_lon > 0:
return 90
# oeste
if delta_lat == 0 and delta_lon < 0:
return 270
if delta_lon > 0 and delta_lat > 0:
return math.degrees(math.atan(delta_lon / delta_lat))
if (delta_lon > 0 > delta_lat) or (delta_lon < 0 and delta_lat < 0):
return math.degrees(math.atan(delta_lon / delta_lat)) + 180
if delta_lon < 0 < delta_lat:
return math.degrees(math.atan(delta_lon / delta_lat)) + 360
azimute = 180 * (1 - (0.5 * sinal_x) - (0.5 * sinal_x * sinal_y)) + (math.degrees(math.atan(x / y)) - (
gamma / 7200))
return azimute
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment