Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
executable file 505 lines (420 sloc) 16.9 KB
#!/usr/bin/env python3
import ephem
from datetime import datetime
import xml.dom.minidom
import json
import math
import sys
from pprint import pprint
import argparse
def nearest_time(targettime, t1, t2):
'''Given a target ephem.time and two other datetimes,
return the time closer to the target.
'''
d1 = abs(targettime - t1)
d2 = abs(targettime - t2)
if d1 <= d2:
return t1
return t2
def find_rise_set(observer, obj, d=None):
'''Given an object (like Sun or Moon), find its rising and setting time
closest to the given date d, either preceding or following it,
for the observer's location.
If date isn't specified, use the observer's date.
'''
if d:
observer.date = d
prevrise = observer.previous_rising(obj)
nextrise = observer.next_rising(obj)
prevset = observer.previous_setting(obj)
nextset = observer.next_setting(obj)
riseset_ret = {}
# Rise time
observer.date = nearest_time(observer.date, prevrise, nextrise)
obj.compute(observer)
riseset_ret['rise'] = { 'az': obj.az / ephem.degree,
'time': observer.date
}
# Set time
observer.date = nearest_time(observer.date, prevset, nextset)
obj.compute(observer)
riseset_ret['set'] = { 'az': obj.az / ephem.degree,
'time': observer.date
}
return riseset_ret
def find_azimuths(observer):
riseset = {}
# Find sunrise and sunset:
riseset['sun'] = find_rise_set(observer, ephem.Sun())
# Now find the full moon closest to the date,
# which may be the next full moon or the previous one.
lastfull = ephem.previous_full_moon(observer.date)
nextfull = ephem.next_full_moon(observer.date)
now = ephem.now()
if now - lastfull > nextfull - now:
observer.date = nextfull
else:
observer.date = lastfull
riseset['full moon'] = find_rise_set(observer, ephem.Moon())
return riseset
def bearing_to(wp1, wp2):
# https://www.movable-type.co.uk/scripts/latlong.html
# Don't trust any code you find for this: test it extensively;
# most posted bearing finding code is bogus.
# print("bearing from waypoints:", wp1, wp2)
lat1, lon1 = math.radians(wp1[1]), math.radians(wp1[2])
lat2, lon2 = math.radians(wp2[1]), math.radians(wp2[2])
y = math.sin(lon2 - lon1) * math.cos(lat2)
x = math.cos(lat1) * math.sin(lat2) - \
math.sin(lat1) * math.cos(lat2) * math.cos(lon2-lon1)
return math.degrees(math.atan2(y, x)) % 360
def find_alignments(observer, waypoints, year, allpoints=False):
'''Find all the alignments with solstice/equinox sun/moon rise/set.
Returns a dict: { 'vernal equinox': { 'moon': { 'rise': 94.17... } } }
of azimuth angles in decimal degrees
'''
azimuths = {}
# start_date = ephem.Date('%d/1/1' % year)
start_date = ephem.Date((year, 1, 1))
# date= ephem.date((-59000,1,1))
observer.date = ephem.next_equinox(start_date)
azimuths['vernal equinox'] = find_azimuths(observer)
observer.date = ephem.next_solstice(observer.date)
azimuths['summer solstice'] = find_azimuths(observer)
observer.date = ephem.next_equinox(observer.date)
azimuths['autumnal equinox'] = find_azimuths(observer)
observer.date = ephem.next_solstice(observer.date)
azimuths['winter solstice'] = find_azimuths(observer)
# pprint(azimuths)
# How many degrees is close enough?
DEGREESLOP = 2.
# If allpoints is set, check angles among all pairs of points.
# Otherwise, only check angles from observer to other points.
if allpoints:
print("Looking for alignments among all points")
observer_points = waypoints
else:
observer_points = [ [ observer.name,
observer.lat / ephem.degree,
observer.lon / ephem.degree,
observer.elevation ] ]
# Now go through all the angles between waypoints and see if
# any of them correspond to any of the astronomical angles.
matches = []
for wp1 in observer_points:
print("\nChecking observer", wp1)
for wp2 in waypoints:
if wp1 == wp2:
continue
angle = bearing_to(wp1, wp2)
print(" ... vs", wp2, angle)
# Does that angle match any of our astronomical ones?
for season in azimuths: # vernal equinox, etc.
for body in azimuths[season]: # sun, full moon
for event in azimuths[season][body]: # rise, set
event_az = azimuths[season][body][event]['az']
if abs(event_az - angle) < DEGREESLOP:
matches.append({
'observer': wp1[0],
'target': wp2[0],
'event': '%s %s%s' % (season, body, event),
'azimuth': event_az,
'slop': event_az - angle,
'time': azimuths[season][body][event]['time'],
'latitude': wp2[1],
'longitude': wp2[2],
})
return matches
def read_waypoint_file_CSV(filename):
"""Read a CSV waypoint file. Ignore tracks.
Return a list of [name, lat, lon, ele] floats for waypoints.
"""
import csv
points = []
observer = None
with open(filename) as csvfp:
reader = csv.DictReader(csvfp)
for row in reader:
# Each row is an OrderedDict
try:
if 'elevation' in row:
ele = float(row['elevation'])
else:
ele = 0
point = [ row['name'],
float(row['latitude']), float(row['latitude']),
ele ]
except ValueError:
print("ValueError on", row)
if row['name'].lower() == "observer":
observer = point
else:
points.append(point)
return observer, points
def read_waypoint_file_GPX(filename):
"""Read a GPX waypoint file. Ignore tracks.
Return a list of [name, lat, lon, ele] floats for waypoints.
"""
dom = xml.dom.minidom.parse(filename)
first_segment_name = None
observer = None
# Handle waypoints
waypts = dom.getElementsByTagName("wpt")
if not waypts:
return []
waypoints = []
pointno = 0
for pt in waypts:
lat = float(pt.getAttribute("lat"))
lon = float(pt.getAttribute("lon"))
try:
ele = float(get_DOM_text(pt, "ele"))
except:
ele = 500 # meters
# Some of these names have newlines in them.
# Just keep the part before the newline.
name = get_DOM_text(pt, "name").split('\n')[0]
# print(" { 'name': '%s', 'lat': %f, 'lon': %f }," % (name, lat, lon))
if not name:
pointno += 1
name = "Point %d" % pointno
if name.lower() == "observer":
observer = [ 'Observer', lat, lon, ele ]
else:
waypoints.append([ name, lat, lon, ele ])
# pprint(waypoints)
return observer, waypoints
def get_DOM_text(node, childname=None):
'''Get the text out of a DOM node.
Or, if childname is specified, get the text out of a child
node with node name childname.
'''
if childname:
nodes = node.getElementsByTagName(childname)
# print "node has", len(nodes), childname, "children"
if not nodes:
return None
node = nodes[0]
if not node:
return None
n = node.childNodes
if len(n) >= 1 and n[0].nodeType == n[0].TEXT_NODE:
return n[0].data
return None
def find_observer_point(obsname, waypoints):
'''Find the waypoint with the name obsname.'''
for wp in waypoints:
if wp[0] == obsname:
return wp
return None
def save_alignments_as_JSON(observer, alignments, waypoints, filename):
'''Given a list of alignments,
save a JSON file that can be plotted in various ways.
'''
print("Saving as JSON")
out = []
for a in alignments:
if observer:
obsname = observer.name
obslat = math.degrees(observer.lat)
obslon = math.degrees(observer.lon)
else:
obs = find_observer_point(a['observer'], waypoints)
if not obs:
print("Can't find observer for", a)
continue
obsname = obs[0]
obslat = obs[1]
obslon = obs[2]
if (a['latitude'] == obslat and
a['longitude'] == obslon):
continue
out.append({ 'observer_name': obsname,
'event': a['event'],
'observer_lat': math.degrees(obslat),
'observer_lon': math.degrees(obslon),
'target_name': a['target'],
'target_lat': a['latitude'],
'target_lon': a['longitude'],
})
with open(filename, 'w') as outfp:
outfp.write(json.dumps(out))
def save_alignments_as_GPX(observer, alignments, waypoints, filename):
'''Given a list of alignments,
save them as a GPX file with tracks between observer and each target,
and a waypoint for each target.
'''
with open(filename, 'w') as outfp:
print('''<?xml version='1.0' encoding='UTF-8' standalone='yes' ?>
<gpx version="1.1" creator="SkyAlignments~" xmlns="http://www.topografix.com/GPX/1/1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
<trk>''', file=outfp)
for a in alignments:
if observer:
obsname = observer.name
obslat = math.degrees(observer.lat)
obslon = math.degrees(observer.lon)
else:
obs = find_observer_point(a['observer'], waypoints)
if not obs:
print("Can't find observer for", a)
continue
obsname = obs[0]
obslat = obs[1]
obslon = obs[2]
if (a['latitude'] == obslat and
a['longitude'] == obslon):
continue
print(''' <trkseg>
<name>%s, %s</name>
<trkpt lat="%f" lon="%f">
<time>%s</time>
<name>%s</name>
</trkpt>
<trkpt lat="%f" lon="%f">
<time>%s</time>
<name>%s, %s</name>
</trkpt>
</trkseg>''' % (a['target'], a['event'],
obslat, obslon,
str(a['time']), obsname,
a['latitude'], a['longitude'], str(a['time']),
a['target'], a['event']),
file=outfp)
print("</trk>", file=outfp)
# Apparently naming trksegs doesn't do anything, at least in pytopo.
# Set waypoints too.
for a in alignments:
if (a['latitude'] == obslat and
a['longitude'] == obslon):
continue
print(''' <wpt lat="%f" lon="%f">
<time>%s</time>
<name>%s, %s</name>
</wpt>''' % (a['latitude'], a['longitude'], str(a['time']),
a['target'], a['event']),
file=outfp)
print("</gpx>", file=outfp)
print("Saved alignments to", filename)
if __name__ == '__main__':
parser = argparse.ArgumentParser(description=""
"""Find alignments between latitude/longitude coordinate pairs
and the sun. moon, and other objects on special dates such as
solstices and equinoxes.
Observer location may be specified either with -o lat,lon,ele or by
naming one of the GPX waypoints 'Observer'; otherwise the first
waypoint in the first file will be taken as the observer location.
When specifying location on the command line, latitude and longitude
are in decimal degrees. Elevation is optional; it will be assumed to be
meters unless followed by the letter f,
e.g. -o 34.8086585,-103.2011914,1650f""",
formatter_class=argparse.RawTextHelpFormatter)
# Specify an observer:
parser.add_argument('-o', '--observer', action="store", dest="observer",
help='Observer location (lat,lon[,ele])')
# Specify an observer name:
parser.add_argument('-n', '--observername', action="store",
dest="observername", help='Observer name')
# Different year from now:
parser.add_argument('-y', '--year', action="store", type=int,
dest="year", help='Year')
# Output JSON file:
parser.add_argument('-j', '--writejson', action="store",
dest="writejson", help='Save JSON to output file')
# Output GPX file:
parser.add_argument('-g', '--writegpx', action="store",
dest="writegpx", help='Save GPX to output file')
# Don't use an observer, check angles between all pairs of points:
parser.add_argument('-a', "--all", dest="allpoints", default=False,
action="store_true",
help="Don't use an observer, check angles "
"between all pairs of points")
parser.add_argument('waypointfiles', nargs='+',
help='GPX or CSV files containing waypoints')
args = parser.parse_args(sys.argv[1:])
if args.observer:
floats = args.observer.split(',')
lat = float(floats[0].strip())
lon = float(floats[1].strip())
if len(floats) > 2:
if floats[2].endswith('f'): # ends with f, convert feet to meters
ele = float(floats[2][:-1].strip()) * 0.3048
elif floats[2].endswith('m'): # ends with m, already meters
ele = float(floats[2][:-1].strip())
else: # assume meters
ele = float(floats[2].strip())
else:
ele = 0.
observer_point = [ 'Observer', lat, lon, ele ]
else:
observer_point = None
if args.year:
year = args.year
else:
year = datetime.now().year
waypoints = []
for filename in args.waypointfiles:
if filename.lower().endswith('gpx'):
obs, wp = read_waypoint_file_GPX(filename)
elif filename.lower().endswith('csv'):
obs, wp = read_waypoint_file_CSV(filename)
else:
print("Unknown file type %s: skipping" % filename)
continue
if wp:
waypoints += wp
else:
print("No waypoints in", filename)
if obs:
observer_point = obs
if not waypoints:
parser.print_help()
sys.exit(1)
if not observer_point:
print("Using first waypoint for observer:", waypoints[0])
observer_point = waypoints[0]
waypoints = waypoints[1:]
# find_alignments will need an initial observer to calculate
# things like the dates of equinoxes and solstices, and the
# azimuths. So set it even if allpoints is set.
observer = ephem.Observer()
# Observer will take degrees as a string, but if you pass it floats
# it expects radians, though that's undocumented.
observer.lat = observer_point[1] * ephem.degree
observer.lon = observer_point[2] * ephem.degree
if len(observer_point) > 3:
observer.elevation = observer_point[3]
else:
observer.elevation = 500.0 # meters
if args.observername:
observer.name = args.observername
else:
observer.name = "%s %f, %f, %dm" % (observer_point[0],
observer.lat / ephem.degree,
observer.lon / ephem.degree,
observer.elevation)
# print(observer)
alignments = find_alignments(observer, waypoints,
year=year, allpoints=args.allpoints)
# Now we no longer need the ephem.observer, and if allpoints is set,
# don't want to pass observer to save_alignments*.
if args.allpoints:
observer = None
if alignments:
# pprint(alignments)
cur_observer = None
for a in alignments:
if a['observer'] != cur_observer:
cur_observer = a['observer']
print("\nFrom %s:" % cur_observer)
print("%s, %s at %s (az %d +/- %.2f)"
% ( a['event'], a['target'],
a['time'],
a['azimuth'], a['slop']))
if args.writejson:
save_alignments_as_JSON(observer, alignments, waypoints,
args.writejson)
if args.writegpx:
save_alignments_as_GPX(observer, alignments, waypoints,
args.writegpx)
else:
print("Couldn't find any alignments with %s" % observer.name)
You can’t perform that action at this time.