Total Eclipse visualization

eclipse_map.py

# Script to plot the path of totality of an eclipse. 
# Uses ephem to calculate the sun eclipse fraction for a
# grid of latitudes, longitudes and altitudes.

import ephem
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import scipy.ndimage
import sys
import re


# Routines to read coordinates from Fred Espenak's tracks.
def dms2dd(degrees, minutes, seconds, direction):
    dd = float(degrees) + float(minutes)/60 + float(seconds)/(60*60);
    if direction == 'W' or direction == 'S':
        dd *= -1
    return dd;

def dd2dms(deg):
    d = int(deg)
    md = abs(deg - d) * 60
    m = int(md)
    sd = (md - m) * 60
    return [d, m, sd]

# Example: dd = parse_dms("78°55'N")
def parse_dms(dms):
    parts = re.split('[°\'"]+', dms)
    # values in Fred's table have no seconds info
    lat = dms2dd(parts[0], parts[1], 0, parts[2]) 
    return (lat)



# Approximate eclipse fraction borrowed from Juha Vierinen. 
def intersection(r0, r1, d, n_s=100):
    A1 = np.zeros([n_s, n_s])
    A2 = np.zeros([n_s, n_s])
    I = np.zeros([n_s, n_s])
    x = np.linspace(-2.0 * r0, 2.0 * r0, num=n_s)
    y = np.linspace(-2.0 * r0, 2.0 * r0, num=n_s)
    xx, yy = np.meshgrid(x, y)
    A1[np.sqrt((xx + d)**2.0 + yy**2.0) < r0] = 1.0
    n_sun = np.sum(A1)
    A2[np.sqrt(xx**2.0 + yy**2.0) < r1] = 1.0
    S = A1 + A2
    I[S > 1] = 1.0
    eclipse = np.sum(I) / n_sun
    return(eclipse)

# Calculate the eclipsed totality for a grid of points, adapted from a
# code by Juha Vierinen.
def get_eclipse(
        t0,
        n_t=1 * 10,
        dt=60.0,
        alts=np.linspace(
            0,
            600e3,
            num=600),
        lats=[69.0],
        lons=[16.02]):
    # Location
    obs = ephem.Observer()
    n_alts = len(alts)
    n_lats = len(lats)
    n_lons = len(lons)

    p = np.zeros([n_t, n_alts, n_lats, n_lons])
    times = np.arange(n_t) * dt
    dts = []
    for ti, t in enumerate(times):
        # print("Time %1.2f (s)" % (t))
        for ai, alt in enumerate(alts):
            obs.elevation = alt
            # (ephem.date(ephem.date(t0)+t*ephem.second))
            obs.date = t0
            sun, moon = ephem.Sun(), ephem.Moon()
            for lai, lat in enumerate(lats):
                for loi, lon in enumerate(lons):
                    obs.lon, obs.lat = '%1.2f' % (lon), '%1.2f' % (lat)  # ESR

                    # Output list
                    results = []
                    seps = []
                    sun.compute(obs)
                    moon.compute(obs)
                    r_sun = (sun.size / 2.0) / 3600.0
                    r_moon = (moon.size / 2.0) / 3600.0
                    s = np.degrees(
                        ephem.separation(
                            (sun.az, sun.alt), (moon.az, moon.alt)))
                    percent_eclipse = 0.0

                    if s < (r_moon + r_sun):
                        # print("eclipsed")
                        if s < 0.85e-2:
                            percent_eclipse = 1.0
                        else:
                            percent_eclipse = intersection(
                                r_moon, r_sun, s, n_s=200)

                    if np.degrees(sun.alt) <= r_sun:
                        if np.degrees(sun.alt) <= -r_sun:
                            percent_eclipse = 1.0
                        else:
                            percent_eclipse = 1.0 - \
                                ((np.degrees(sun.alt) + r_sun) / (2.0 * r_sun)) * (1.0 - percent_eclipse)

                    p[ti, ai, lai, loi] = percent_eclipse
        dts.append(obs.date)
    return(p, times, dts)


# Plot cartography, projected tracks and umbra + penumbra contour maps. 
# Totality_path.txt can be copied from Fred Espenak's website or NASA's eclipse database. 
def plot_map(p, lons, lats, t0, alt=100e3, lat_0=69, lon_0=16, width=8e6, height=8e6, totname = 'totality_path.txt', scale = 1.0, cities=False):

    fig = plt.figure(figsize=(9, 9))
    plt.style.use('dark_background')

    # create polar stereographic Basemap instance.
    m = Basemap(projection='stere', lon_0=lon_0, lat_0=lat_0, resolution='h', width=width, height=height)

    # Choose from etopo, shaderelief or bluemarble texture for the map.
    # m.etopo()
    m.shadedrelief(scale = scale)
    # m.bluemarble(scale = scale)

    # Draw coastlines, state and country boundaries, edge of map.
    countriescolor = (47/255, 79/255, 79/255)
    countriescolor = (169/255, 169/255, 169/255)
    countriescolor = (0/255, 0/255, 0/255)
    m.drawcoastlines(color=countriescolor, linewidth=0.5)
    m.drawcountries(color=countriescolor)
    m.drawstates(color=countriescolor)
    # m.drawrivers(color='tab:blue')

    # Draw parallels and meridians. 
    parallels = np.arange(-90.,0., 20.)
    # m.drawparallels(parallels,labels=[1,0,0,1])
    m.drawparallels(parallels)
    meridians = np.arange(10.,360., 30.)
    # m.drawmeridians(meridians,labels=[1,0,0,1])
    m.drawmeridians(meridians)
    
    lon, lat = np.meshgrid(lons, lats)
    # Compute map proj coordinates.
    x, y = m(lon, lat) 

    # Draw the umbra and penumbra as filled contours.
    pp = p[0, 0, :, :]
    # pp[pp<0.49] = np.nan
    # cs = m.pcolormesh(x, y, pp, alpha = .5, edgecolors=None, 
    #                   vmin=0, vmax=1.0, cmap="inferno_r", zorder=1)

    levels = np.linspace(0.1, 1.0, num=10, endpoint=True)
    cs = m.contourf(x, y, pp, levels, cmap = "magma_r", alpha=0.6, zorder=1)
    m.contour(x, y, pp, levels, colors='black', linewidths=0.5, zorder=1)

    # Draw totality.
    levels = [0.997, 1.0]
    m.contourf(x, y, pp, levels, colors='black', alpha=0.66, zorder=1)
    # m.contour(x, y, pp, levels, colors='white', linewidths=1.0)

    # Read in totality path lines from Fred Espenak at NASA.
    (nlat0, nlon0, slat0, slon0, clat0, clon0) = np.loadtxt(totname, dtype=np.dtype('str'), skiprows=3, usecols=(0,1,2,3,4,5), unpack=True)

    nlat = np.array([parse_dms(xi) for xi in nlat0])
    nlon = np.array([parse_dms(xi) for xi in nlon0])
    slat = np.array([parse_dms(xi) for xi in slat0])
    slon = np.array([parse_dms(xi) for xi in slon0])
    clat = np.array([parse_dms(xi) for xi in clat0])
    clon = np.array([parse_dms(xi) for xi in clon0])

    # Northern limit of totality.
    m.plot(nlon, nlat,  color='white', linestyle='dashed', latlon=True, zorder=11, linewidth=0.8)
    # Southern limit of totality.
    m.plot(slon, slat,  color='white', linestyle='dashed', latlon=True, zorder=11, linewidth=0.8)
    # Central line of totality.
    m.plot(clon, clat, linewidth=1., color='tab:red', latlon=True, zorder=11)


    if cities:
        plt.text(width/2 + width/50, height/2, 'Villarrica', zorder=12)
        plt.plot([width/2], [height/2], "wo")
        
    # add colorbar.
    cbar = m.colorbar(cs, location='bottom', pad="4%")

    # Plot title.
    # plt.title("Total Solar Eclipse\n%s (UTC)" % (t0))
    fname = "eclipse-snapshots/eclipse-%f.jpg" % (float(t0))
    # plt.text(width/20, height/20,"HagoCiencia.cl \nUSACH",color="black")
    plt.text(width/20, height/20, "Total Solar Eclipse \n%s (UTC) \nUSACH" %(t0), color="black")

    print("saving %s" % (fname))
    plt.savefig(fname, bbox_inches='tight', dpi=300, trasparent=True)
    plt.close()



# Examples.

# ****************************************
# Total Eclipse Chile 2020 (Temuco coords)
# ****************************************
lat_0 = -39.28
lon_0 = -72.23
date = (2020,12,14,16,0,0)
date = (2020,12,14,15,50,0)
 
# # Latitude and longitude domains for calculation.
# nlats = 100
# nlons = 200
# lats = np.linspace(lat_0-36.,lat_0+36.,num=nlats)
# lons = np.linspace(lon_0-66.,lon_0+66.,num=nlons)

# # Make animation:
# minute=35
# seconds=0
# for i in range(20):
#     t0=ephem.date((2020,12,14,15,minute,seconds))
#     (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
#     plot_map(p, lons, lats, t0, alt=0., lat_0=lat_0,lon_0=lon_0, width=10e6, height=10e6)
#     seconds += 30
    
# Latitude and longitude domains for calculation.
nlats = 200
nlons = 300
lats = np.linspace(lat_0-15.,lat_0+15.,num=nlats)
lons = np.linspace(lon_0-22.,lon_0+22.,num=nlons)
 
# Zoom in:
minute=45
seconds=0
for i in range(20):
    t0=ephem.date((2020,12,14,15,minute,seconds))
    (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
    plot_map(p, lons, lats, t0, alt=0., lat_0=lat_0,lon_0=lon_0, width=3e6, height=3e6)
    seconds += 30

# Latitude and longitude domains for calculation.
nlats = 400
nlons = 600
lats = np.linspace(lat_0-4.,lat_0+4.,num=nlats)
lons = np.linspace(lon_0-4.,lon_0+4.,num=nlons)

# Zoom in:
minute=55
seconds=0
for i in range(60):
    t0=ephem.date((2020,12,14,15,minute,seconds))
    (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
    plot_map(p, lons, lats, t0, alt=0., lat_0=lat_0,lon_0=lon_0, width=0.6e6, height=.6e6, scale = 2.0, cities=True)
    seconds += 15
    

# Latitude and longitude domains for calculation.
nlats = 200
nlons = 300
lats = np.linspace(lat_0-36.,lat_0+36.,num=nlats)
lons = np.linspace(lon_0-66.,lon_0+66.,num=nlons)

# Make animation:
minute=70
seconds=0
for i in range(15):
    t0=ephem.date((2020,12,14,15,minute,seconds))
    (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
    plot_map(p, lons, lats, t0, alt=0., lat_0=lat_0,lon_0=lon_0, width=10e6, height=10e6)
    seconds += 15



# # ****************************************
# # Total Eclipse Chile 2019 (La Serena)
# # ****************************************
# lat_0 = -29.9027
# lon_0 = -71.2519
# date = (2019,7,2,20,40,0)
# # date = (2019,7,2,20,28,0)

# # Latitude and longitude domains for calculation.
# nlats = 300
# nlons = 600
# lats = np.linspace(lat_0-40.,lat_0+40.,num=nlats)
# lons = np.linspace(lon_0-66.,lon_0+66.,num=nlons)

# # Make animation:
# minute = 38
# for i in range(1):
#     minute += 1
#     t0=ephem.date((2019,7,2,20,minute,0))
#     (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
#     plot_map(p, lons, lats, t0, alt=0., lat_0=lat_0,lon_0=lon_0, width=9e6, height=9e6, totname = 'totality_path_2019.txt')
    

# # Total Eclipse USA 2017
# lat_0 = 38.496077
# lon_0 = 360 - 90.152751
# date = (2017,8,21,17,30,0)


# Animation.
# install ffmpeg with brew
# ffmpeg -framerate 6 -pattern_type glob -i "eclipse-snapshots/*.jpg" -s:v 1440x1080 -c:v libx264 -crf 17 -pix_fmt yuv420p my-timelapse.mp4

Snapshot example: