Instantly share code, notes, and snippets.

# pklaus/sunplot.py

Last active February 8, 2024 19:24
Show Gist options
• Save pklaus/52ed125b781d635975c7 to your computer and use it in GitHub Desktop.
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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 A nice online tool: http://bl.ocks.org/mbostock/7784f4b2c7838b893e9b """ 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() # 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}} """ def date_parser(string): return datetime.strptime(string, '%Y-%m-%d').date() 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 today (overrides --include-date).') parser.add_argument('--include-date', type=date_parser, default=None, help='Include a date to choose freely.') 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)}) elif args.include_date: lines.append({'varname': 'TODAY', 'date': args.include_date}) plot(lines, args.lat, args.lon, args.tzoffset) if __name__ == "__main__": main()

### RaoulAbrutat commented Oct 30, 2016

How does a Python 3 write command look like to write the above sun path data (JUNE, DECEMBER and TODAY) to a file?
I understand that write() takes exactly 1 argument.

I would highly appreciate your feedback (this is my 1. exposure to Python). Thank you.

Observation:
By-the-by, for gnuplot 5.0.4 to plot the diagram "at" in set key, and "e" after line {{TODAY}} was missing.