<p>In the desktop version of Stellarium, you can do this on the Ephemeris tab of the Astronomical calculations window.
Set the beginning of the date range with the desired time of day,
and make the time step a whole number of solar days.
Check the box H.C. (horizontal coordinates), otherwise you get a ring around the ecliptic.
You can write a CSV file using the Export ephemeris... dialog.</p>
<p><img src="https://i.sstatic.net/pyFjh.png" alt="Analemma generated by Stellarium" /></p>

CC BY-SA 4.0
2023-06-27T21:38:07.830Z
2023-06-27T21:38:07.830Z
Mike G
<p>Here's a Sage / Python program which creates analemma diagrams using data from JPL <a href="https://ssd.jpl.nasa.gov/horizons/" rel="nofollow noreferrer">Horizons</a>. The azimuth + elevation data is displayed in 3D on a globe, which can be rotated &amp; zoomed. (Sage uses three.js for interactive 3D graphics).</p>
<p>Please note that this globe displays the azimuth and elevation, so its poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe).</p>
<p>You can supply the location using coordinates in longitude, latitude, altitude form, or you can supply the site ID using the <a href="https://en.wikipedia.org/wiki/List_of_observatory_codes" rel="nofollow noreferrer">MPC code</a> of an observatory. The altitude, given in kilometres, is relative to the WGS-84 reference ellipsoid (as used by GPS), so it's very close to altitude relative to mean sea level.</p>
<p>If you supply a site ID the data in the <code>lon, lat, alt</code> field is ignored. To specify <code>lon, lat, alt</code>, put <code>coord</code> in the <code>site</code> field, or leave it blank. If you put <code>*</code> in the site field, Horizons will list all the MPC sites. You can also fetch the list with <a href="https://ssd.jpl.nasa.gov/api/horizons.api?format=text&amp;COMMAND=0&amp;CENTER=%27*@399%27&amp;OBJ_DATA=NO" rel="nofollow noreferrer">this URL</a>.</p>
<p>To create a correct analemma you need to specify the time zone as part of the start time field. Eg, to make a noon analemma for longitude 30°E, the time should be 12:00 UT+2 or 12:00 UT+2:00.</p>
<p>Horizons can estimate the effects of atmospheric refraction, if desired. (The true refraction depends on the weather). The effects of refraction are fairly small around noon, but are significant near dawn and dusk.</p>
<p>Here's the relevant entry from the Horizons manual:</p>
<blockquote>
<h3>4. Apparent azimuth &amp; elevation (AZ-EL) :</h3>
<p>Apparent azimuth and elevation of target. Adjusted for light-time, the gravitational deflection of light, stellar aberration, precession and nutation. There is an optional (approximate) adjustment for atmospheric refraction (Earth only).</p>
<p>Azimuth is measured clockwise from north: North(0) -&gt; East(90) -&gt; South(180) -&gt; West(270) Elevation angle is with respect to plane perpendicular to local zenith direction.<br />
TOPOCENTRIC ONLY. Units: DEGREES</p>
</blockquote>
<p>Here's a noon analemma for Greenwich observatory, using a 7 day timestep, and a datestep of 4 (so every 4th dot is larger). This plot is for 2023. It doesn't use refraction.</p>
<p><a href="https://i.sstatic.net/JKkf1.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JKkf1.png" alt="Analemma globe" /></a></p>
<p>The program can display numeric labels near the datestep points, but they tend to clutter the diagram, so I toggled the labels off for that example.</p>
<p>Here's the code.</p>
<pre><code>&quot;&quot;&quot; Retrieve solar azimuth &amp; elevation from Horizons for a specified Earth location
    and plot it in 3D on a globe.

    Written by PM 2Ring 2024.01.13
&quot;&quot;&quot;

import re, requests
from functools import lru_cache

url = &quot;https://ssd.jpl.nasa.gov/api/horizons_file.api&quot;
api_version = &quot;1.0&quot;

# Greenwich Lon, Lat, Alt (km)
greenwich = '0.0, 51.4773207, 0.06707'

base_cmd = &quot;&quot;&quot;
MAKE_EPHEM=YES
EPHEM_TYPE=OBSERVER
QUANTITIES=4
CAL_FORMAT=CAL
CAL_TYPE=GREGORIAN
CSV_FORMAT=YES
OBJ_DATA=NO
COMMAND=10
&quot;&quot;&quot;

@lru_cache(maxsize=int(3))
def fetch_data(site, coord, start, stop, step, refracted, verbose=True):       
    if site == &quot;&quot;:
        site = &quot;coord&quot;
    apparent = &quot;REFRACTED&quot; if refracted else &quot;AIRLESS&quot;
    
    cmd = f&quot;&quot;&quot;
!<span class="math-container">$$SOF
{base_cmd}
CENTER='{site}@399'
SITE_COORD='{coord}'
START_TIME='{start}'
STOP_TIME='{stop}'
STEP_SIZE='{step}'
APPARENT='{apparent}'
!$$</span>EOF
&quot;&quot;&quot;
    #print(cmd)
    req = requests.post(url, data={'format': 'text'}, files={'input': ('cmd', cmd)})
    version = re.search(r&quot;API VERSION:\s*(\S*)&quot;, req.text).group(1)
    if version != api_version:
        print(f&quot;Warning: API version is {version}, but this script is for {api_version}&quot;)

    m = re.search(r&quot;(?s)\\\$\\\$SOE(.*)\\\$\\\$EOE&quot;, req.text)
    if m is None:
        print(&quot;NO EPHEMERIS&quot;)
        print(req.text)
        return None

    if verbose:
        print(req.text)
    else:
        lines = req.text.splitlines()
        print(&quot;\n&quot;.join(lines[5:19]), &quot;\n&quot;)

    return m.group(1)[1:]

def extract_data(data, datestep):
    data = data.splitlines()
    data = [s.split(',')[:-1] for s in data]
    #for row in data: print(row)
    #return

    pos, dates = [], []
    for i, row in enumerate(data):
        # Extract azimuth &amp; elevation
        az, el = float(row[3]), float(row[4])
        pos.append((az, el))

        # Make timestamps
        if datestep and i % datestep == 0:
            # Extract calendar date
            ymd, hms = row[0].split()
            dates.append(ymd)

    return pos, dates

# Globe with lat-lon grid
var('u,v')
pp = parametric_plot3d
deg = pi / 180
dth = 15 * deg

def xyz(lat, lon):
    r = cos(lat)
    return vector([r*cos(lon), r*sin(lon), sin(lat)])

udom, vdom = (u, -pi/2, pi/2), (v, -pi, pi)
G = pp(xyz(u, v), udom, vdom, color='#007', plot_points=75)
G += sum(pp(xyz(u, lon), udom, color='#00c' if lon != pi else '#0ac') for lon in srange(0, pi*2, dth))
G += sum(pp(xyz(lat, v), vdom, color='#00c' if lat else '#0ac') for lat in srange(dth-pi/2, pi/2, dth))

@interact
def main(txt=HtmlBox(&quot;&lt;h3&gt;Analemma globe&lt;/h3&gt;&quot;),
  site='', coord=('lon, lat, alt (km)', greenwich),
  start=&quot;2023-Jan-01 12:00&quot;, stop=&quot;2024-Jan-01&quot;, step=&quot;7d&quot;,
  datestep=4, show_labels=True, refracted=False, color=&quot;red&quot;,
  verbose=False, auto_update=False):

    result = fetch_data(site.strip(), coord.strip(),
      start.strip(), stop.strip(), step.strip(), refracted, verbose=verbose)
    if result is None:
        print(&quot;No ephemeris data found!&quot;)
        return

    pos, dates = extract_data(result, datestep)
    #for row in pos: print(row)

    for i, s in enumerate(dates):
        print(i, s)

    # Make a blank globe
    P = G

    for i, (az, el) in enumerate(pos):
        do_label = datestep and i % datestep == 0
        v = xyz(el * deg, -az * deg)
        P += point3d(v, size=5 if do_label else 3, color=color)
        if show_labels and do_label:
            label = str(i // datestep)
            vv = v * 1.1
            P += text3d(label, vv, fontsize=8)

    P.show(frame=False, theme='dark', projection='orthographic', online=True)
</code></pre>
<p>The data-fetching part of the code is plain Python. The Sage features are only needed to do the 3D display.</p>
<p>Here's a <a href="https://sagecell.sagemath.org/?z=eJyNV21z2kYQ%5Fq5fcZHTSHJlGYxdN0zVhtiKQxsbCrSd1vYwZ3GAEkmn3p2wY4%5F%5Fe3fvJBAk7ZQBvezu7T2379i2TUZMiYStGJE8pYLQxyQr1ZK8IixlK6oSnpO54Bl5z0XyyHNJ5hykiCxYnMwTNiMRFSCf8lgLWwQ%2DNJ%2DRIuWKJPDNSeecgBZKFim%5FY4GlRf4QiVIsJ3efyfCSHI2SfEGOWkfHQasdtDuWbduWlWQFF4oI5sPv75JJJS2NZV7mseI8laQSSUU5jWm8ZJZVipSExF4qVcju4aGUs%2DBjkQY5lTRY8NUhLZLDZXWW6TxJWQAU24LLdMWExPPC8nbQAgB75EIwlt8n8ZJ84LlPPlDlk16qiPsp86zFmhkSpxW0fHLSDo5PTztHrVOfAOG709apg2r%2DFxoE8hNYN6MqVOxBvTobXF72rs7D1quz6GoSjcJvjk7333Rev4b7q8Hbn6fnvUkvvBpY1h2VbBpnM4QOlrvs%5FRJNo%2DH76DL8Mxpb%2Dmk6%2DXMYhYO342j0ezSyfv2tdzXpT%5FrRODy2znofpu8Go8veJIRH%5FaqlL0bRxWDU711ZZ%2DPfawnU2Ny8RtluGa%2D9WTvDzeiDTB5ZmOTK7XieNWNzMmcqXk5nVFFXJgpcG3MuZj6RCuIIb7zAKyvQ6XNBY8WAC66545KFE1Eyr0vMRwdSMieoh4R49K5VsSoasbV220RlUVDBcoXkUfRu1DubROc2KlhvBEEvGbF7%5FdGHaDw2y%5FTF2HaOB3zx8uV48M56qm3%2DbFXecZ5w02d0kGON%2D5NoejYYjM6BrkE8A3HSG02mk%5F5lhMJ4YE0cDDc0XmhSNJyO%2D38ZEkNSbzjsjWAfoNTnACpgiQALwkKUe4VAUwMoT79D1gDqOneCgkvlQn74BM0fPjkm2JwucTDenGefYEJI4CR5USLDdUCZ4%2DP5vWejdJMlggWSUREvXWH3hn0CcTXuD666N3LfvRnve7bO2wB1e8FC8LJw217ttFrNi5A0cm%5FjQHOUuf0HFTkUhy7BHepFiSRP1TOAvisVUUugyVgkhUIuFqmnht5n2zN1J9vB7f4kvZuXN%2DDSyA329WM0iJrIa8AZ6r3iOdvFaF8NiE6xaNQf294Od1uP8YoqRa5VWQ1rYHzvqt5ejMG5kUiTnEnjXi0TyCJNlKa6uyDsm9wOPvIkdzX%5F%2DqTbfn3r%2DQTplWEqVNnaUdft7q2lUxaUY3qYpMWLDiCGkekZPEgEKHj7EkbFvJaG5Tq%2D4113D9q32ksSOwSK3JoQRprg9zW1W1uC3xttewaoAQ0RXWHBDW59%2DGk6Kkn8Wg%5FLy4wJkNLgvY0J90hkjva1treWoo8%2DUDH9U041kusOGm%5FzenzbMDjHQl6wfOa6ZqVXWdjseEk%5FMaKSDMxHs0KuORAEtVF1%5F0zINxsCFLfWBvY29JimsBm0bpTekvmcQeVcZjpIAGXrtnKAtyWlN6khf8bS0YyHjYV1O8QOTu4TbPhUHaSQiguRzKwVFa5T%2DivHswpAS6BA0QwHi3iKc0BnBnG0QHpCDkn7%2D5Y1U9g02ydknwDHhNnD50c3xQYLaisnCRCKuUSy10S1YtD%5FhXst9jUX5MHZ%2DxLjWz%5FrJ1gDfrHKGc%2DggcAVlLmlTw6g5x75BK8g6q40Bd896wIxFi4iAcEVsDersVOlXITOXgu6uq8HnGkBOaVkeHqCa78NiSwzd6PAgCl3FscOOjs1tQ8sonsOMGjseDpykQVhKwXNF8xtIbR9AAw2877cRlsMke5irLah6isbUNXYAPQ2TFLvY72BkzGMMO2cjIJJ1YMK36ssfcsfXPuHZefHXg7Bl2XVaPfDIZBsz7dM%2Dw0dp%2DruoeukOD5prLQan4C5HqDMGmyIoQ1TYOfgZ5oftNqkfdRttWwzFWjOccWxzYwQ2qczG9fWmRIeA2PJ76cpvYNj63mhMUiE7yjYoraSLZhZXU8XFZeWik%5FLAnUaEkRjFX2yTHF%2D2BliAgmhXrheddz1a5Vn%2DmAbGTxM84013r4y8lT3dR%2DqQPxrM%2DKEFUsGBQ8kdOmd8zKfvbB3G9BXKuhWoTf7NEr9F%5FUZ1m6V52bplV8UXia9XbAoVy2rCiMldynNP5l40owh4LrYUl2X1e0dAExD%5F4ybCDBN6T%2Dq6nrFCkQxn2CNrklQFeijedxYbojJp3O%2DM8PCoafbE1276w11snXqENNXr1nkG9GpEdULt8t7DR4Cw4WaebjjhTVoRL0ClPCnaYuhgeJUADi1LognwDvnUK0Q8%5FeV2YcB4nHnWK7r8FcYPaEDLeUTljnBP0K1hYYYOvBXa8kXghbLJAYWz7HNm5nc%2DgenGUH2&amp;lang=sage" rel="nofollow noreferrer">live version of the program</a>, running on the SageMathCell server.</p>
<p>The program maintains a cache  of its last 3 queries of Horizons data, so if you only change &quot;cosmetic&quot; options the plot is generated without re-fetching the data from Horizons.</p>
<p>Here are the 3D interface controls.</p>
<ul>
<li>Orbit - right mouse, or left mouse + ctrl/meta/shiftKey
<ul>
<li>touch: one-finger rotate</li>
</ul>
</li>
<li>Zoom - middle mouse, or mousewheel
<ul>
<li>touch: two-finger spread or squish</li>
</ul>
</li>
<li>Pan - left mouse, or arrow keys
<ul>
<li>touch: two-finger move</li>
</ul>
</li>
</ul>

CC BY-SA 4.0
2024-01-13T16:41:46.497Z
Awesome. It can be called, the AIO Earth Analemma generator  :). But can you Label the coordinate system and directions in the diagram? Or only numbers allowed? Because it would be easy for layman to understand the diagram.
Jacob Miller
CC BY-SA 4.0
2024-01-13T16:53:51.070Z
@Jacob I suppose I *could* add coordinates, but they would get in the way, especially for noon analemmas. The grid steps are 15°. The "equator" is 0° elevation. The prime "meridian" is the azimuth=180° (south) line. They're in a lighter, greener shade than the other azimuth & elevation circles, but I guess that may be not so easy to see.
PM 2Ring
CC BY-SA 4.0
2024-01-13T17:11:51.060Z
So this view is outside the celestial sphere. Then if you made the globe transparent and looked from the opposite side (in the direction of the north at noon), view correspond to how observers on Earth see?
Jacob Miller
CC BY-SA 4.0
2024-01-13T17:40:29.080Z
@Jacob Right, but you can't actually zoom the camera into the sphere. I guess it's possible by playing with the viewpoint parameter... But it's easy to make the globe transparent: just add something like `opacity=0.5` to the function call that creates the globe: `G = pp(xyz(u, v), ...`
PM 2Ring
CC BY-SA 4.0
2024-01-13T18:00:10.277Z
I think if there is a way to apply some solid colour below the horizon but still with lines visible that would give a better perspective for a novice to analemma. By the way, Transparency works fine, lines are not too much of a bother.
Jacob Miller
CC BY-SA 4.0
2024-01-16T10:44:28.707Z
Re: "The "equator" is 0° elevation.":  0° elevation depicted in the diagram (horizontal greenish line in the default view) is the observer's local horizon which will be different from the equator unless the observer is in the North or South poles. If it had been the celestial equator then whichever the time the analemma shall be drawn with its axis orthogonal to the equator.
Jacob Miller
CC BY-SA 4.0
2024-01-16T11:15:49.260Z
@Jacob Right, it's not the celestial equator, which is why I wrote "equator" in quotes. My globe displays azimuth & elevation, so the poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe). I suppose I should clarify that in my answer...
PM 2Ring
CC BY-SA 4.0
2024-01-13T12:38:59.220Z
2024-01-16T11:20:29.930Z
2024-01-16T11:20:29.930Z
PM 2Ring
<p>If you're open to using a bit of Python, you can use the following script to generate an analemma for arbitrary times and locations.</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
import matplotlib.projections
import astropy.units as u
import astropy.time
import astropy.coordinates
import astropy.visualization
import cartopy

time = astropy.time.Time(&quot;2023-06-26T13:30&quot;) + np.arange(365) * u.day
location = astropy.coordinates.EarthLocation.of_address(&quot;Boston&quot;)
frame = astropy.coordinates.AltAz(obstime=time, location=location)

sun = astropy.coordinates.get_body(&quot;sun&quot;, time=time).transform_to(frame)

alt = sun.alt.deg
az = sun.az.deg

projection = cartopy.crs.Orthographic(az.mean(), alt.mean())
fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
ax.plot(az, alt, transform=cartopy.crs.PlateCarree());
ax.gridlines(draw_labels=True)
plt.show()
</code></pre>
<p>which outputs:</p>
<p><a href="https://i.sstatic.net/EB3q7.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/EB3q7.png" alt="analemma for Boston" /></a></p>

CC BY-SA 4.0
2023-06-27T13:54:30.970Z
Nice. Thanks. Can you provide a tweak to your code for a projection where the meridian lines converge to the pole, and the shape of the analemma doesn’t “droop” as much relative to a more symmetrical figure-8? This might make the graph look more like the di Cicco photograph. Perhaps polar coordinates, or maybe there’s a projection that models the field of view of a camera? (It’s been quite a while since I coded graphs in python.)
Bruce Simonson
CC BY-SA 4.0
2023-06-27T15:58:51.920Z
@BruceSimonson, I've updated the answer with an improved projection. Let me know if that's closer to what you're looking for.
Roy Smart
CC BY-SA 4.0
2023-06-27T18:10:34.337Z
Do you know if astropy takes the ellipsoid into account?   The angles might be slightly different based on being closer/further from the Equator (thus Earth's center).
Spencer
CC BY-SA 4.0
2023-06-27T18:51:12.633Z
Astropy does take the ellipsoid into account, but the `of_address()` method does not compute the height above the ellipsoid unless you provide a Google API key. More information here https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html
Roy Smart
CC BY-SA 4.0
2023-06-27T20:39:26.210Z
@Roy Smart - excellent projection choice. Thank you.
Bruce Simonson
CC BY-SA 4.0
2024-01-14T04:06:18.877Z
What are those horizontal coordinates (110deg E, 120deg E) on the output diagram? They seem wrong for the time.
Jacob Miller
CC BY-SA 4.0
2024-01-14T18:08:16.623Z
@JacobMiller, the horizontal coordinate is the local azimuthal coordinate in Boston, MA. Why do they seem wrong? Stellarium produces roughly the same result for 13:30 UTC using MikeG's answer.
Roy Smart
CC BY-SA 4.0
2024-01-14T22:49:30.047Z
I asked this earlier because the code seemingly produced some abnormal results for different attitudes and times. Later it became clear that it was a combination of using Boston as a location and UT in the code. So the code is correct for Boston.  I just applied it incorrectly for other locations :) Also why not add plt show line at the end for outputting the graph?
Jacob Miller
CC BY-SA 4.0
2024-01-14T23:13:49.263Z
Ah yes, it is a continual annoyance to me that Astropy does not support timezones, but oh well. I meant to add the `plt.show()`, but I must've forgotten. In Jupyter notebooks the `plt.show()` is implicit which is why I forgot.
Roy Smart
CC BY-SA 4.0
2023-06-27T03:41:03.560Z
2024-01-14T23:14:06.653Z
2024-01-14T23:14:06.653Z
Roy Smart
<p>I have in mind to do a yearlong photo project of the solar analemma.</p>
<p>Can someone suggest a software program that will provide (&amp; chart) the altitude and azimuth of the sun, for a given location, at a given time of day throughout the year?</p>
<p>I’d like a graph, and perhaps a table of values.</p>
<p>For example, for Dennis di Cicco’s famous analemma photograph, taken from Boston, at 8:30 AM (throughout the course of a year), the parameters would be</p>
<p>Boston (42.3601° N, 71.0589° W), and
8:30 AM (13:30 UT).</p>
<p>I would like to be able to use a different location and time of day than the di Cicco photograph.</p>
<p>(I have Stellarium on my iphone, if there’s a way to use it).</p>
<p><a href="https://i.sstatic.net/jBfwZ.jpg" rel="noreferrer"><img src="https://i.sstatic.net/jBfwZ.jpg" alt="Analemma, Dennis di Cicco" /></a></p>

CC BY-SA 4.0
2023-06-27T15:40:39.750Z
Another project that's similar that might interest you is a solarigraph.  It's a continuous 6 month exposure.  It has a higher success rate as one cloudy day won't ruin a year's worth of work.
Greg Miller
<modeling><analemma>
analemma for a specified lat/long at a specific time of day?
Bruce Simonson
# analemma for a specified lat/long at a specific time of day?

(Asked by: username_0 on 2023-06-27)

I have in mind to do a yearlong photo project of the solar analemma.

Can someone suggest a software program that will provide (& chart) the altitude and azimuth of the sun, for a given location, at a given time of day throughout the year?

I’d like a graph, and perhaps a table of values.

For example, for Dennis di Cicco’s famous analemma photograph, taken from Boston, at 8:30 AM (throughout the course of a year), the parameters would be

Boston (42.3601° N, 71.0589° W), and
8:30 AM (13:30 UT).

I would like to be able to use a different location and time of day than the di Cicco photograph.

(I have Stellarium on my iphone, if there’s a way to use it).

[![Analemma, Dennis di Cicco](https://i.sstatic.net/jBfwZ.jpg)](https://i.sstatic.net/jBfwZ.jpg)

**Tags:** modeling, analemma
> **Comments:**
> **username_1**: Another project that's similar that might interest you is a solarigraph.  It's a continuous 6 month exposure.  It has a higher success rate as one cloudy day won't ruin a year's worth of work.  

---

## 3 answers

### Answer by username_5 on 2023-06-27. Score: 12 | Accepted answer

If you're open to using a bit of Python, you can use the following script to generate an analemma for arbitrary times and locations.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.projections
import astropy.units as u
import astropy.time
import astropy.coordinates
import astropy.visualization
import cartopy

time = astropy.time.Time("2023-06-26T13:30") + np.arange(365) * u.day
location = astropy.coordinates.EarthLocation.of_address("Boston")
frame = astropy.coordinates.AltAz(obstime=time, location=location)

sun = astropy.coordinates.get_body("sun", time=time).transform_to(frame)

alt = sun.alt.deg
az = sun.az.deg

projection = cartopy.crs.Orthographic(az.mean(), alt.mean())
fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
ax.plot(az, alt, transform=cartopy.crs.PlateCarree());
ax.gridlines(draw_labels=True)
plt.show()

which outputs:

> **Comments:**
> **username_0**: Nice. Thanks. Can you provide a tweak to your code for a projection where the meridian lines converge to the pole, and the shape of the analemma doesn’t “droop” as much relative to a more symmetrical figure-8? This might make the graph look more like the di Cicco photograph. Perhaps polar coordinates, or maybe there’s a projection that models the field of view of a camera? (It’s been quite a while since I coded graphs in python.)  
> **username_5**: @BruceSimonson, I've updated the answer with an improved projection. Let me know if that's closer to what you're looking for.  
> **username_6**: Do you know if astropy takes the ellipsoid into account?   The angles might be slightly different based on being closer/further from the Equator (thus Earth's center).  
> **username_5**: Astropy does take the ellipsoid into account, but the `of_address()` method does not compute the height above the ellipsoid unless you provide a Google API key. More information here https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html  
> **username_0**: @Roy Smart - excellent projection choice. Thank you.  
> **username_4**: What are those horizontal coordinates (110deg E, 120deg E) on the output diagram? They seem wrong for the time.  
> **username_5**: @JacobMiller, the horizontal coordinate is the local azimuthal coordinate in Boston, MA. Why do they seem wrong? Stellarium produces roughly the same result for 13:30 UTC using MikeG's answer.  
> **username_4**: I asked this earlier because the code seemingly produced some abnormal results for different attitudes and times. Later it became clear that it was a combination of using Boston as a location and UT in the code. So the code is correct for Boston.  I just applied it incorrectly for other locations :) Also why not add plt show line at the end for outputting the graph?  
> **username_5**: Ah yes, it is a continual annoyance to me that Astropy does not support timezones, but oh well. I meant to add the `plt.show()`, but I must've forgotten. In Jupyter notebooks the `plt.show()` is implicit which is why I forgot.  

---

### Answer by username_2 on 2023-06-27. Score: 9

In the desktop version of Stellarium, you can do this on the Ephemeris tab of the Astronomical calculations window.
Set the beginning of the date range with the desired time of day,
and make the time step a whole number of solar days.
Check the box H.C. (horizontal coordinates), otherwise you get a ring around the ecliptic.
You can write a CSV file using the Export ephemeris... dialog.

![Analemma generated by Stellarium](https://i.sstatic.net/pyFjh.png)


---

### Answer by username_3 on 2024-01-13. Score: 3

Here's a Sage / Python program which creates analemma diagrams using data from JPL Horizons. The azimuth + elevation data is displayed in 3D on a globe, which can be rotated & zoomed. (Sage uses three.js for interactive 3D graphics).
Please note that this globe displays the azimuth and elevation, so its poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe).
You can supply the location using coordinates in longitude, latitude, altitude form, or you can supply the site ID using the MPC code of an observatory. The altitude, given in kilometres, is relative to the WGS-84 reference ellipsoid (as used by GPS), so it's very close to altitude relative to mean sea level.
If you supply a site ID the data in the lon, lat, alt field is ignored. To specify lon, lat, alt, put coord in the site field, or leave it blank. If you put * in the site field, Horizons will list all the MPC sites. You can also fetch the list with this URL.
To create a correct analemma you need to specify the time zone as part of the start time field. Eg, to make a noon analemma for longitude 30°E, the time should be 12:00 UT+2 or 12:00 UT+2:00.
Horizons can estimate the effects of atmospheric refraction, if desired. (The true refraction depends on the weather). The effects of refraction are fairly small around noon, but are significant near dawn and dusk.
Here's the relevant entry from the Horizons manual:

4. Apparent azimuth & elevation (AZ-EL) :
Apparent azimuth and elevation of target. Adjusted for light-time, the gravitational deflection of light, stellar aberration, precession and nutation. There is an optional (approximate) adjustment for atmospheric refraction (Earth only).
Azimuth is measured clockwise from north: North(0) -> East(90) -> South(180) -> West(270) Elevation angle is with respect to plane perpendicular to local zenith direction.
TOPOCENTRIC ONLY. Units: DEGREES

Here's a noon analemma for Greenwich observatory, using a 7 day timestep, and a datestep of 4 (so every 4th dot is larger). This plot is for 2023. It doesn't use refraction.

The program can display numeric labels near the datestep points, but they tend to clutter the diagram, so I toggled the labels off for that example.
Here's the code.
""" Retrieve solar azimuth & elevation from Horizons for a specified Earth location
    and plot it in 3D on a globe.

    Written by PM 2Ring 2024.01.13
"""

import re, requests
from functools import lru_cache

url = "https://ssd.jpl.nasa.gov/api/horizons_file.api"
api_version = "1.0"

# Greenwich Lon, Lat, Alt (km)
greenwich = '0.0, 51.4773207, 0.06707'

base_cmd = """
MAKE_EPHEM=YES
EPHEM_TYPE=OBSERVER
QUANTITIES=4
CAL_FORMAT=CAL
CAL_TYPE=GREGORIAN
CSV_FORMAT=YES
OBJ_DATA=NO
COMMAND=10
"""

@lru_cache(maxsize=int(3))
def fetch_data(site, coord, start, stop, step, refracted, verbose=True):       
    if site == "":
        site = "coord"
    apparent = "REFRACTED" if refracted else "AIRLESS"
    
    cmd = f"""
!$$SOF
{base_cmd}
CENTER='{site}@399'
SITE_COORD='{coord}'
START_TIME='{start}'
STOP_TIME='{stop}'
STEP_SIZE='{step}'
APPARENT='{apparent}'
!$$EOF
"""
    #print(cmd)
    req = requests.post(url, data={'format': 'text'}, files={'input': ('cmd', cmd)})
    version = re.search(r"API VERSION:\s*(\S*)", req.text).group(1)
    if version != api_version:
        print(f"Warning: API version is {version}, but this script is for {api_version}")

    m = re.search(r"(?s)\\\$\\\$SOE(.*)\\\$\\\$EOE", req.text)
    if m is None:
        print("NO EPHEMERIS")
        print(req.text)
        return None

    if verbose:
        print(req.text)
    else:
        lines = req.text.splitlines()
        print("\n".join(lines[5:19]), "\n")

    return m.group(1)[1:]

def extract_data(data, datestep):
    data = data.splitlines()
    data = [s.split(',')[:-1] for s in data]
    #for row in data: print(row)
    #return

    pos, dates = [], []
    for i, row in enumerate(data):
        # Extract azimuth & elevation
        az, el = float(row[3]), float(row[4])
        pos.append((az, el))

        # Make timestamps
        if datestep and i % datestep == 0:
            # Extract calendar date
            ymd, hms = row[0].split()
            dates.append(ymd)

    return pos, dates

# Globe with lat-lon grid
var('u,v')
pp = parametric_plot3d
deg = pi / 180
dth = 15 * deg

def xyz(lat, lon):
    r = cos(lat)
    return vector([r*cos(lon), r*sin(lon), sin(lat)])

udom, vdom = (u, -pi/2, pi/2), (v, -pi, pi)
G = pp(xyz(u, v), udom, vdom, color='#007', plot_points=75)
G += sum(pp(xyz(u, lon), udom, color='#00c' if lon != pi else '#0ac') for lon in srange(0, pi*2, dth))
G += sum(pp(xyz(lat, v), vdom, color='#00c' if lat else '#0ac') for lat in srange(dth-pi/2, pi/2, dth))

@interact
def main(txt=HtmlBox("Analemma globe"),
  site='', coord=('lon, lat, alt (km)', greenwich),
  start="2023-Jan-01 12:00", stop="2024-Jan-01", step="7d",
  datestep=4, show_labels=True, refracted=False, color="red",
  verbose=False, auto_update=False):

    result = fetch_data(site.strip(), coord.strip(),
      start.strip(), stop.strip(), step.strip(), refracted, verbose=verbose)
    if result is None:
        print("No ephemeris data found!")
        return

    pos, dates = extract_data(result, datestep)
    #for row in pos: print(row)

    for i, s in enumerate(dates):
        print(i, s)

    # Make a blank globe
    P = G

    for i, (az, el) in enumerate(pos):
        do_label = datestep and i % datestep == 0
        v = xyz(el * deg, -az * deg)
        P += point3d(v, size=5 if do_label else 3, color=color)
        if show_labels and do_label:
            label = str(i // datestep)
            vv = v * 1.1
            P += text3d(label, vv, fontsize=8)

    P.show(frame=False, theme='dark', projection='orthographic', online=True)

The data-fetching part of the code is plain Python. The Sage features are only needed to do the 3D display.
Here's a live version of the program, running on the SageMathCell server.
The program maintains a cache  of its last 3 queries of Horizons data, so if you only change "cosmetic" options the plot is generated without re-fetching the data from Horizons.
Here are the 3D interface controls.

Orbit - right mouse, or left mouse + ctrl/meta/shiftKey

touch: one-finger rotate


Zoom - middle mouse, or mousewheel

touch: two-finger spread or squish


Pan - left mouse, or arrow keys

touch: two-finger move

> **Comments:**
> **username_4**: Awesome. It can be called, the AIO Earth Analemma generator  :). But can you Label the coordinate system and directions in the diagram? Or only numbers allowed? Because it would be easy for layman to understand the diagram.  
> **username_3**: @Jacob I suppose I *could* add coordinates, but they would get in the way, especially for noon analemmas. The grid steps are 15°. The "equator" is 0° elevation. The prime "meridian" is the azimuth=180° (south) line. They're in a lighter, greener shade than the other azimuth & elevation circles, but I guess that may be not so easy to see.  
> **username_4**: So this view is outside the celestial sphere. Then if you made the globe transparent and looked from the opposite side (in the direction of the north at noon), view correspond to how observers on Earth see?  
> **username_3**: @Jacob Right, but you can't actually zoom the camera into the sphere. I guess it's possible by playing with the viewpoint parameter... But it's easy to make the globe transparent: just add something like `opacity=0.5` to the function call that creates the globe: `G = pp(xyz(u, v), ...`  
> **username_4**: I think if there is a way to apply some solid colour below the horizon but still with lines visible that would give a better perspective for a novice to analemma. By the way, Transparency works fine, lines are not too much of a bother.  
> **username_4**: Re: "The "equator" is 0° elevation.":  0° elevation depicted in the diagram (horizontal greenish line in the default view) is the observer's local horizon which will be different from the equator unless the observer is in the North or South poles. If it had been the celestial equator then whichever the time the analemma shall be drawn with its axis orthogonal to the equator.  
> **username_3**: @Jacob Right, it's not the celestial equator, which is why I wrote "equator" in quotes. My globe displays azimuth & elevation, so the poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe). I suppose I should clarify that in my answer...  

---

