Skip to content

Instantly share code, notes, and snippets.

@bbinet
Forked from pklaus/sunplot.py
Created August 7, 2014 16:48
Show Gist options
  • Save bbinet/f1348c994a92446fbe2e to your computer and use it in GitHub Desktop.
Save bbinet/f1348c994a92446fbe2e to your computer and use it in GitHub Desktop.
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment