sunplot.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
This tool can draw the path of the sun for a given
position on earth. It will draw two lines:
One for June 21 and another for December 21.

You need to have PyEphem and gnuplot installed
for this script to work properly.

Interesting circles of latitude:
+- 66:33:44  Polar Circle
+- 23:26:16  Tropic of Cancer / Tropic of Capricorn 
    0:00:00  Equator


How to plot it with matplotlib:
http://stackoverflow.com/questions/12858806/stereographic-sun-diagram-matplotlib-polar-plot-python

Convert the radial presentation to make it stereographic?
http://www.l-e-s-s.co.uk/Guides/Physics/SolarGeometry.htm#28
"""

import ephem
import math
from datetime import datetime, date, time, timedelta

sun = ephem.Sun()

def values(date, lat, lon, tz_offset):
    ob = ephem.Observer()
    plotstep = timedelta(minutes=5)
    ob.lat, ob.lon = lat, lon
    ts = datetime.combine(date, time(12)) - timedelta(hours=tz_offset)
    ob.date = ts
    try:
        sunrise = ob.previous_rising(sun).datetime()
        noon = ob.next_transit(sun, start=sunrise).datetime()
        sunset = ob.next_setting(sun).datetime()
        #print(sunrise, noon, sunset)
    except (ephem.AlwaysUpError, ephem.NeverUpError):
        sunrise = ts - timedelta(hours=12)
        sunset  = ts + timedelta(hours=12)
    ret = []
    ts = sunrise
    while ts < sunset:
        ret.append(alt_az(ts, ob))
        ts += plotstep
    ret.append(alt_az(sunset, ob))
    return ret

def alt_az(ts, ob):
    ob.date = ts
    sun.compute(ob)
    return (sun.alt*180/math.pi, sun.az*180/math.pi)

def transform(data):
    ret = []
    for coords in data:
        if coords[0] < -0.5: continue
        ret.append((90.0-coords[1], 90-coords[0]))
    return ret

def plot(lines, lat, lon, tz_offset):
    from subprocess import Popen,PIPE
    gnuplot = 'gnuplot'
    plot = Popen([gnuplot, '-persist'], stdin=PIPE, stdout=PIPE, stderr=PIPE)
    command = template

    command = command.replace("{{TITLE}}", "Sun Path for (Lat: %s, Lon: %s)" % (lat,lon) )

    for line in lines:
        data = transform(values(line['date'], lat, lon, tz_offset))
        command = command.replace("{{" + line['varname'] + "}}", "\n".join("%f %f"%d for d in data))
        command = command.replace("{{LEGEND_" + line['varname'] + "}}", line['date'].strftime('%B %d %Y'))

    plot.stdin.write(command.encode())

    plot.stdin.flush()

def main():
    import argparse
    parser = argparse.ArgumentParser(description='Plot the course of the sun. The location defaults to the RGO.')
    parser.add_argument('--lat', type=ephem.degrees, default=ephem.degrees('51:28:40'), help='Latitude of the place to plot. North is +.')
    parser.add_argument('--lon', type=ephem.degrees, default=ephem.degrees('-0:0:5'), help='Longitude of the place to plot. East is +.')
    parser.add_argument('--tzoffset', type=int, help='Offset of the local time compared to UTC in full hours. + means ahead of UTC.')
    parser.add_argument('--year', type=int, default=date.today().year, help='The year you want to take a look at')
    parser.add_argument('--include-today', action='store_true', help='Include a line for the the day/month today')
    args = parser.parse_args()

    args.lat = args.lat.znorm # should even restrict to +- 90.0 deg
    args.lon = args.lon.znorm

    if not args.tzoffset:
        args.tzoffset = round((args.lon / math.pi) * 12)

    june = date(args.year, 1, 1).replace(month=6,day=21)
    december = date(args.year, 1, 1).replace(month=12,day=21)

    lines = [{'varname': 'JUNE', 'date': june}, {'varname': 'DECEMBER', 'date': december}]
    if args.include_today: lines.append({'varname': 'TODAY', 'date': date.today().replace(year=args.year)})

    plot(lines, args.lat, args.lon, args.tzoffset)

# The gnuplot command template (data will be inserted by Python):
template = """
set encoding utf8
#set term svg enhanced mouse size 600,400
#set output "sunpath.svg"

set angles degrees
set polar
set size ratio 1
set tmargin 3
set bmargin 3

set style line 11 lc rgb 'gray80' lt -1
set grid polar ls 11
set grid
set key 140,100

unset border
unset xtics
unset ytics

r=90
set rrange [0:r]
set rtics 10 format '' scale 0

set label '0° - North' center at first 0, first r*1.05
set label '180° - South' center at first 0, first -r*1.05
set label '270° - West' right at first -r*1.05, 0
set label '90° - East' left at first r*1.05, 0

set label '{{TITLE}}' left at first -r*1.4, first r*1.15

set for [i=1:8] label at first r*0.02, first r*((i/9.0) + 0.03) sprintf("%d°", 90-i*10)
unset raxis

#plot 50*(1+sin(t)) linewidth 2 t ''

#plot "-" u 1:2 smooth bezier title '{{LEGEND_JUNE}}', "-" u 1:2 smooth bezier title '{{LEGEND_DECEMBER}}', "-" u 1:2 smooth bezier title '{{LEGEND_TODAY}}'
plot "-" u 1:2 w l title '{{LEGEND_JUNE}}', "-" u 1:2 w l title '{{LEGEND_DECEMBER}}', "-" u 1:2 w l title '{{LEGEND_TODAY}}'
{{JUNE}}
e
{{DECEMBER}}
e
{{TODAY}}
"""

if __name__ == "__main__":
    main()