Generate iCalendar .ics files with events for astrological aspects

Question

I'm relatively new to Python, coming from a deep C++ background. I'm mostly looking for feedback on how to make my code more idiomatic/pythonic, but I would welcome and appreciate any and all other feedback as well.

This code uses the Python library PyEphem, which builds upon XEphem, an X (UNIX GUI) program that's an astrological ephemeris, which is (was historically) a table of celestial bodies and their positions in the sky. PyEphem is used to track the altitudes of astronomical bodies, which are used as an indication of when astrological aspects have or are going to occur.

I'm deploying this code as an AWS Lambda microservice, thus the lambda_handler function. Consider also reviewing the test and deployment Makefile that corresponds to the module below.

Here is the entirety of my module:

#!/usr/bin/env python2.7

# aspectus: aspect+prospectus. iCalendars with astrological aspect events.
# Copyright (C) 2017 Frederick Eugene Aumson
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""aspectus: aspect+prospectus. Generates iCalendar .ics files to populate your
calendar application with events for astrological aspects (angles) between two
celestial bodies.  Currently supports only the trine and sextant aspects, and
only between the sun and the earth. Could easily be extended to support other
bodies, and perhaps with not much more work extended to support aspects between
two bodies outside of Earth.
"""

from dateutil.tz import gettz
from geopy.geocoders import Nominatim
from ephem import Observer, Sun  # pylint:disable=no-name-in-module
from datetime import timedelta
from math import degrees, radians
from icalendar import Calendar, Event

def find_altitude(targets, body, observer, step, latest_date):
    """returns a dict with 'target' and 'datetime' properties, where 'datetime'
    is the next point in time (UTC) at which the ephem Body :param body will be
    found at one of the altitudes in the list :param targets, and 'target' is
    the altitude of the body at that datetime.  begins the search at the (UTC!)
    datetime contained by ephem Observer :param observer and proceeds searching
    by stepping the observer's datetime forward by datetime.timedelta :param
    step in each iteration. :param targets should be a list of floats
    representing altitudes in degrees.  :param observer should be already
    initialized with a datetime and location. body.compute(observer) should
    already have been called.  Returns None if :param latest_date is reached
    without finding a target.
    """

    #print targets, body, observer, step

    altitude = degrees(float(repr(body.alt)))

    while True:
        # capture each target's relation to altitude before taking step forward
        before = [target < altitude for target in targets]

        # step observer date forward
        observer.date = observer.date.datetime() + step
        if observer.date.datetime() > latest_date:
            #print 'reached latest date'
            return None
        body.compute(observer)
        altitude = degrees(float(repr(body.alt)))
        #print altitude, observer

        # if any target is "close enough", return it
        for target in targets:
            if abs(target - altitude) < abs(target * 0.00001):
                # step forward to prevent double-catching this target
                observer.date = observer.date.datetime() + \
                        timedelta(seconds=step.total_seconds()*10)
                body.compute(observer)
                #print 'found', altitude, 'at', observer
                return {'target': target, 'datetime': observer.date.datetime()}

        # capture each target's relation to altitude after step forward
        after = [target < altitude for target in targets]

        # a target whose relation changed was passed by the step forward.
        # if there is such a target, step the observer date backwards, to
        # "un-pass" the target, and then recurse, with a finer-grained step.
        changes = [delta[0] != delta[1] for delta in zip(before, after)]
        for index, changed in enumerate(changes):
            if changed:
                observer.date = observer.date.datetime() - step
                body.compute(observer)
                return find_altitude([targets[index]], body, observer,
                        timedelta(seconds=step.total_seconds()/10), latest_date)

def generate_icalendar(place, lookaheaddays, start):
    """generate icalendar containing events for all sun alignments at the place
    described by string :param place, looking by int()-compatible string :param
    lookaheaddays number of days, optionally starting from the optional
    dateutil.parser-compatible string :param start, which is interpreted as the
    datetime with local time (NOT UTC) in the timezone at :param place, as
    determined by tzwhere.
    """

    location = Nominatim().geocode(place)

    observer = Observer()
    observer.lat = radians(location.latitude)
    observer.lon = radians(location.longitude)
    if start is not None:
        from dateutil.parser import parse
        observer.date = parse(start).astimezone(gettz('UTC'))

    startdate = observer.date.datetime()

    sun = Sun()
    sun.compute(observer)

    cal = Calendar()
    cal.add('prodid',
            '-//[email protected]//'
            'https://github.com/feuGeneA/aspectus//EN')
            # embed revision number into URL above
    cal.add('version', '2.0') # iCalendar spec version, not prodid version
    cal_has_events = False

    while observer.date.datetime() < \
            startdate + timedelta(days=int(lookaheaddays)):
        find = find_altitude([30, -30], sun, observer, timedelta(hours=1),
                startdate + timedelta(days=int(lookaheaddays)))
        if find == None:
            break
        event = Event()
        if find['target'] == 30:
            event.add('summary', 'Sun-earth Trine, at '+place)
        elif find['target'] == -30:
            event.add('summary', 'Sun-earth Sextant, at '+place)
        event.add('dtstart', find['datetime'].replace(tzinfo=gettz('UTC')))
        event.add('dtend', \
                find['datetime'].replace(tzinfo=gettz('UTC'))\
                +timedelta(seconds=1))

        cal.add_component(event)
        cal_has_events = True

    if cal_has_events:
        return cal
    else:
        return None

def lambda_handler(event, context): # pylint:disable=unused-argument
    """handle event from AWS lambda"""

    print 'event: ', event
    print 'context: ', context

    cal = generate_icalendar( \
            event['queryStringParameters']['place'],
            event['queryStringParameters']['lookaheaddays'],
            event['queryStringParameters']['startdatetime'] \
                if 'startdatetime' in event['queryStringParameters'] else None)

    if cal is None:
        return {'statusCode': '204'}
    else:
        return {
            'statusCode': '200',
            'body': cal.to_ical(),
            'headers': {
                'Content-Type': 'text/calendar',
                'Content-Disposition': \
                        'attachment; filename="Sun-earth trines at '+\
                                event['queryStringParameters']['place']+'.ics"'
            },
        }

def main():
    from sys import argv
    cal = generate_icalendar(argv[1], argv[2], None)
    print cal.to_ical()

if __name__ == "__main__":
    main()

"""Backlog of improvements to this file:
As a user, I want to pass calendar-generation options via the command line, so
that I can generate a calendar by running the python script directly, rather
than relying on an online instance of the service.

As a tester of the web service, I want any Python exceptions to be conveyed as
part of the body of an HTTP 500 response, so that if something goes wrong I
have some indication of what the problem was.

As a user, I want to specify a number of minutes ahead of the alignment at
which I will be notified by my calendar application, so that I don't get
stuck with the application's default (30 minutes for Google Calendar), or,
worse yet, so that I don't have no alarm at all. (Add VALARM components to
the events.)

As a user, rather than the alignment events being a single moment in time, I
want them to have a duration, using an astrological "orb" based on the size
of the body in angular distance, so that I understand how long the alignment
is having an effect.  (Use pyehem.body.size, which is diameter in
arcseconds)

As a maintainer, I want the commented-out print statements in
find_altitude() to replaced with debug/trace logger statements, so that I
can enable/disable them in one place, via logger configuration, rather than
hunting for and changing each statement.
""" # PEP258: Additional Docstrings # pylint:disable=pointless-string-statement

I have two unit test modules, which live in a tests subdirectory alongside the module above. Here is the first:

#!/usr/bin/env python2.7

import unittest

import os
import sys
sys.path.insert(0, os.path.abspath('..'))
import aspectus

from datetime import datetime, timedelta
import icalendar

def mock_find_altitude(targets, body, observer, step, latest_date):
    observer.date = observer.date.datetime()+timedelta(days=1)
    return {'target': targets[0], 'datetime': observer.date.datetime()}

aspectus.find_altitude = mock_find_altitude

class TestGenerateiCalendar(unittest.TestCase):
    def test_one_day_out(self):
        cal = aspectus.generate_icalendar('avon nc', '1', None)
        self.assertIsInstance(cal, icalendar.Calendar)
        self.assertTrue(len(cal.content_lines()) > 0)

if __name__ == "__main__":
    unittest.main()

And here is the second unit test module:

#!/usr/bin/env python2.7

import unittest

import os
import sys
sys.path.insert(0, os.path.abspath('..'))
import aspectus

import icalendar

def mock_gen_ical_returning_calendar(place, lookaheaddays, start):
    cal = icalendar.Calendar()
    return cal

def mock_gen_ical_returning_None(place, lookaheaddays, start):
    return None

class TestLambdaHndler(unittest.TestCase):
    def test_good(self):
        aspectus.generate_icalendar = mock_gen_ical_returning_calendar
        resp = aspectus.lambda_handler(
            {
                'httpMethod': 'GET',
                'queryStringParameters': {
                    'place': 'Avon, North Carolina',
                    'lookaheaddays': '1',
                    'startdatetime': '2017/01/01 00:00:00Z-0500'
                }
            },
            None)
        self.assertIn('statusCode', resp)
        self.assertEqual(resp['statusCode'], '200')

Answer 1

Just reviewing the function find_altitude.

1. Review

1. The docstring is comprehensive, which is good.

2. But the specification is very complicated, which makes the function hard to use. In particular, the caller has to remember to set observer.date to the start of the span to be searched, and to call body.compute(observer). It would be easy to forget to do one or both of these, and then the function would go wrong. It would be simpler if find_altitude did both of these tasks. This would require adding an extra parameter giving the start of the period to be searched.

3. In Python, it is conventional to use the names start, stop, and step (in that order) when describing an arithmetic progression: see for example the arguments to range and slice.

4. When a function needs to return multiple results (as here, where you have an altitude and a datetime) it is usually most convenient to return them as a tuple (rather than a dictionary). That's because the caller can conveniently use tuple unpacking to assign the results to local variables, like this:

   altitude, datetime = find_altitude(...)
   

   If you like having a result with named fields, then you can use a collections.namedtuple, which would give you the best of both worlds.

5. I think altitudes would be a better name than targets.

6. In the case where none of the altitudes is crossed, it would be a good idea to raise an exception rather than returning None. The reason is that it would be easy for the caller to forget to check for the exceptional case, possibly leading to incorrect results, whereas an exception is unambiguous.

7. Formatting the docstring into paragraphs would improve its readability. See below for how I would write it (though you might need to mark it up with :param: and whatnot).

8. Here you go through a rather winding path to convert body.alt to degrees, by converting it to a string representation, parsing that as a float, and then converting that to degrees:

   altitude = degrees(float(repr(body.alt)))
   

   It would be better to take a shorter path:

   altitude = degrees(body.alt)
   

   But an alternative approach would be to convert all the target altitudes to radians once on input, avoiding the need for further conversions during the search.

9. Here you compare altitude with each target altitude in turn, taking \$Θ(n)\$ if there are \$n\$ target altitudes:

   before = [target < altitude for target in targets]
   

   But if the target altitudes were sorted into numerical order, then you would only need to know the position of altitude in the sorted list, which can be found in \$O(\log n)\$ using bisect.bisect. If the position in the sorted list changes, then you know that body has crossed one (or more) of the altitudes.

10. Each time around the loop you convert observer.date to a datetime twice:

    while True:
        # ...
        observer.date = observer.date.datetime() + step
        if observer.date.datetime() > latest_date:
    

    but it would be simpler to convert the inputs into PyEphem data structures once, and then no conversions would be required in the loop:

    date = ephem.Date(start)
    stop = ephem.Date(stop)
    step = step.total_seconds() * ephem.second
    
    while True:
        # ...
        date += step
        if date > stop:
            raise NotFound()
        observer.date = date
    

11. Once you've found a timestep in which the body has crossed a target altitude, then you can use PyEphem's newton function to efficiently find the exact crossing point. This uses the Newton-Raphson method to get quadratic convergence (the number of correct digits roughly doubles in each step), whereas the the approach in the post only has linear convergence.

2. Revised code

import ephem
from bisect import bisect
from math import radians

class NotFound(Exception):
    pass

def find_altitude(altitudes, body, observer, start, stop, step):
    """Find the first altitude crossed by a body in a time period.

    Parameters
    altitudes: list(float) -- the target altitudes in degrees
    body: Body -- the body
    observer: Observer -- observer from whose location the
        altitude is calculated.
    start: datetime -- start of the period to be searched
    stop: datetime -- end of the period to be searched
    step: timedelta -- size of the search steps

    Returns a pair (altitude, datetime) whose first element is the first
    altitude crossed by the body in the time period, and whose second
    element is the datetime at which that happens.

    If none of the altitudes is crossed in the time period, raise
    NotFound.

    """
    orig_altitudes = sorted(altitudes)

    # Convert inputs to PyEphem data structures.
    altitudes = list(map(radians, orig_altitudes))
    date = ephem.Date(start)
    stop = ephem.Date(stop)
    step = step.total_seconds() * ephem.second

    def alt(d):
        # Return altitude of body as seen by observer at date d.
        observer.date = d
        body.compute(observer)
        return body.alt

    # Position of alt(date) in sorted list of altitudes.
    pos = bisect(altitudes, alt(date))
    while True:
        # Advance current date forward by one step.
        old_date = date
        date += step
        if date > stop:
            raise NotFound()

        old_pos = pos
        pos = bisect(altitudes, alt(date))
        if old_pos != pos:
            # At least one of the target altitudes was crossed.
            break

    # Find the index i of the first target altitude that was crossed.
    if old_pos < pos:       # Moving forward in list.
        i = old_pos
    else:                   # Moving backward in list.
        i = old_pos - 1

    # Solve alt(d) == altitudes[i] using Newton's method.
    def f(d):
        return alt(d) - altitudes[i]
    crossing_date = ephem.newton(f, old_date, date)

    # Convert outputs back again.
    return orig_altitudes[i], ephem.Date(crossing_date).datetime()

3. General point

By default, PyEphem uses the standard astronomical epoch J2000.0. But is that appropriate for this use case? For example, if you need results computed using epoch-of-date, then change:

body.compute(observer)

to:

body.compute(observer, epoch=d)