seba-perez/eclipse_map.py

## 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.

# Seba Perez (CIRAS/USACH)

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 (copy/paste from website).
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 2021 \n%s (UTC) \nCIRAS USACH" %(t0), color="white")

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

    # free memory
    x=0; y=0; pp=0
    return True


# ****************************************
# Total Eclipse Chile 2021 (Antarctica)
# ****************************************
lat_0 = -76.776
lon_0 = -46.52
date = (2021,12,4,7,30,0)

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

# Make animation:
minute=35
seconds=0

# Single snapshot.
t0=ephem.date((2021,12,4,7,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, totname='totality_2021.txt')


# # Generate snapshots in parallel for animation.
# from multiprocessing import Pool
# from contextlib import closing

# def generate_snapshot(seconds):
#     seconds *= 10
#     t0=ephem.date((2021,12,4,7,minute,seconds))
#     (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
#     print("get_eclipse done ",seconds)
#     plot_map(p,lons,lats,t0,alt=0.,lat_0=lat_0,lon_0=lon_0, width=10e6, height=10e6, totname='totality_2021.txt')
#     # free memory
#     p = 0; dts = 0; times = 0
#     print("plot_map done ",seconds)
#     return True

# if __name__ == '__main__':

#     with closing(Pool(processes=4)) as o:
#         print(o.map(generate_snapshot, range(180), chunksize=4))


# 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

## totality_2021.txt
# Total Solar Eclipse of 2021 Dec 04
# UT1	Northern Limit	Southern Limit	Central Line	M:S	Sun	Sun	Path	Central
# Latitude	Longitude	Latitude	Longitude	Latitude	Longitude	Ratio	Alt.	Azm.	Width	Duration
51°52.7'S	048°57.8'W	54°22.7'S	053°46.1'W	53°05.0'S	051°12.3'W	1.031	0°	-	414 km	01m27.8s
# 58°26.7'S	038°53.9'W	-	-	56°17.5'S	045°48.3'W	1.033	5°	124°	440 km	01m34.4s
60°15.5'S	037°07.1'W	55°24.4'S	051°55.7'W	58°47.2'S	042°48.3'W	1.034	7°	121°	450 km	01m38.8s
61°52.2'S	035°48.3'W	58°48.6'S	047°24.5'W	60°42.1'S	041°02.6'W	1.034	9°	119°	452 km	01m41.8s
63°21.3'S	034°48.3'W	60°54.4'S	045°32.2'W	62°22.1'S	039°50.7'W	1.035	11°	118°	452 km	01m44.1s
64°45.3'S	034°02.7'W	62°39.6'S	044°25.3'W	63°53.5'S	039°00.7'W	1.035	12°	116°	450 km	01m46.1s
66°05.5'S	033°29.1'W	64°13.9'S	043°45.5'W	65°19.0'S	038°27.6'W	1.035	13°	115°	447 km	01m47.8s
67°22.7'S	033°06.4'W	65°41.0'S	043°25.9'W	66°40.1'S	038°08.7'W	1.036	14°	114°	444 km	01m49.2s
68°37.6'S	032°54.0'W	67°03.0'S	043°23.5'W	67°57.9'S	038°03.0'W	1.036	15°	113°	441 km	01m50.4s
69°50.6'S	032°51.8'W	68°21.0'S	043°37.3'W	69°13.1'S	038°10.2'W	1.036	15°	112°	437 km	01m51.4s
71°02.0'S	033°00.5'W	69°35.7'S	044°07.1'W	70°25.9'S	038°30.7'W	1.036	16°	112°	434 km	01m52.2s
72°12.0'S	033°21.0'W	70°47.5'S	044°53.8'W	71°36.8'S	039°05.5'W	1.036	16°	112°	431 km	01m52.9s
73°20.8'S	033°54.7'W	71°56.7'S	045°58.8'W	72°45.8'S	039°56.2'W	1.037	17°	112°	428 km	01m53.5s
74°28.4'S	034°43.6'W	73°03.2'S	047°24.2'W	73°53.1'S	041°05.1'W	1.037	17°	112°	425 km	01m53.9s
75°34.8'S	035°50.7'W	74°07.1'S	049°12.9'W	74°58.5'S	042°35.2'W	1.037	17°	113°	422 km	01m54.2s
76°40.0'S	037°19.9'W	75°07.9'S	051°28.7'W	76°01.8'S	044°30.6'W	1.037	17°	114°	420 km	01m54.4s
77°43.8'S	039°16.4'W	76°05.4'S	054°15.9'W	77°02.8'S	046°56.6'W	1.037	17°	115°	418 km	01m54.4s
78°45.7'S	041°47.3'W	76°58.6'S	057°39.7'W	78°00.8'S	049°59.5'W	1.037	17°	117°	417 km	01m54.3s
79°45.2'S	045°02.2'W	77°46.7'S	061°45.6'W	78°55.0'S	053°47.5'W	1.037	17°	120°	415 km	01m54.1s
80°41.3'S	049°13.4'W	78°28.4'S	066°38.4'W	79°44.3'S	058°29.4'W	1.037	17°	124°	414 km	01m53.7s
81°32.7'S	054°36.1'W	79°01.8'S	072°20.9'W	80°26.9'S	064°14.0'W	1.037	16°	129°	414 km	01m53.2s
82°17.2'S	061°26.9'W	79°25.1'S	078°51.6'W	81°00.8'S	071°06.7'W	1.036	16°	135°	413 km	01m52.6s
82°51.8'S	069°57.8'W	79°35.9'S	086°02.2'W	81°23.2'S	079°05.4'W	1.036	16°	142°	413 km	01m51.8s
83°13.0'S	080°05.3'W	79°32.0'S	093°36.6'W	81°31.4'S	087°54.6'W	1.036	15°	150°	413 km	01m50.9s
83°17.1'S	091°17.8'W	79°11.5'S	101°12.6'W	81°23.0'S	097°05.0'W	1.036	14°	158°	414 km	01m49.7s
83°01.5'S	102°35.4'W	78°32.6'S	108°26.8'W	80°56.3'S	105°59.0'W	1.036	13°	166°	415 km	01m48.4s
82°25.6'S	112°55.0'W	77°33.5'S	114°59.6'W	80°10.7'S	114°03.6'W	1.035	13°	173°	416 km	01m46.8s
81°30.4'S	121°36.4'W	76°10.4'S	120°37.9'W	79°05.6'S	120°58.3'W	1.035	11°	180°	418 km	01m45.0s
80°16.9'S	128°29.3'W	74°14.6'S	125°14.0'W	77°39.4'S	126°35.5'W	1.034	10°	184°	420 km	01m42.8s
78°44.8'S	133°41.4'W	71°09.1'S	128°35.3'W	75°46.7'S	130°55.2'W	1.034	8°	188°	422 km	01m40.0s
# 76°50.5'S	137°24.7'W	-	-	73°09.9'S	133°54.7'W	1.033	6°	190°	424 km	01m36.2s
67°03.1'S	138°44.1'W	67°35.9'S	129°04.9'W	67°21.9'S	134°10.7'W	1.031	0°	-	418 km	01m27.8s
	# 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.

	# Seba Perez (CIRAS/USACH)

	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 (copy/paste from website).
	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 + yy2.0) < r0] = 1.0
	n_sun = np.sum(A1)
	A2[np.sqrt(xx2.0 + yy2.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 2021 \n%s (UTC) \nCIRAS USACH" %(t0), color="white")

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

	# free memory
	x=0; y=0; pp=0
	return True




	# ****************************************
	# Total Eclipse Chile 2021 (Antarctica)
	# ****************************************
	lat_0 = -76.776
	lon_0 = -46.52
	date = (2021,12,4,7,30,0)

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

	# Make animation:
	minute=35
	seconds=0

	# Single snapshot.
	t0=ephem.date((2021,12,4,7,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, totname='totality_2021.txt')


	# # Generate snapshots in parallel for animation.
	# from multiprocessing import Pool
	# from contextlib import closing

	# def generate_snapshot(seconds):
	# seconds *= 10
	# t0=ephem.date((2021,12,4,7,minute,seconds))
	# (p, times, dts) = get_eclipse(t0, n_t=1, dt=60.0, alts=[0.], lats=lats, lons=lons)
	# print("get_eclipse done ",seconds)
	# plot_map(p,lons,lats,t0,alt=0.,lat_0=lat_0,lon_0=lon_0, width=10e6, height=10e6, totname='totality_2021.txt')
	# # free memory
	# p = 0; dts = 0; times = 0
	# print("plot_map done ",seconds)
	# return True

	# if __name__ == '__main__':

	# with closing(Pool(processes=4)) as o:
	# print(o.map(generate_snapshot, range(180), chunksize=4))



	# 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
	# Total Solar Eclipse of 2021 Dec 04
	# UT1 Northern Limit Southern Limit Central Line M:S Sun Sun Path Central
	# Latitude Longitude Latitude Longitude Latitude Longitude Ratio Alt. Azm. Width Duration
	51°52.7'S 048°57.8'W 54°22.7'S 053°46.1'W 53°05.0'S 051°12.3'W 1.031 0° - 414 km 01m27.8s
	# 58°26.7'S 038°53.9'W - - 56°17.5'S 045°48.3'W 1.033 5° 124° 440 km 01m34.4s
	60°15.5'S 037°07.1'W 55°24.4'S 051°55.7'W 58°47.2'S 042°48.3'W 1.034 7° 121° 450 km 01m38.8s
	61°52.2'S 035°48.3'W 58°48.6'S 047°24.5'W 60°42.1'S 041°02.6'W 1.034 9° 119° 452 km 01m41.8s
	63°21.3'S 034°48.3'W 60°54.4'S 045°32.2'W 62°22.1'S 039°50.7'W 1.035 11° 118° 452 km 01m44.1s
	64°45.3'S 034°02.7'W 62°39.6'S 044°25.3'W 63°53.5'S 039°00.7'W 1.035 12° 116° 450 km 01m46.1s
	66°05.5'S 033°29.1'W 64°13.9'S 043°45.5'W 65°19.0'S 038°27.6'W 1.035 13° 115° 447 km 01m47.8s
	67°22.7'S 033°06.4'W 65°41.0'S 043°25.9'W 66°40.1'S 038°08.7'W 1.036 14° 114° 444 km 01m49.2s
	68°37.6'S 032°54.0'W 67°03.0'S 043°23.5'W 67°57.9'S 038°03.0'W 1.036 15° 113° 441 km 01m50.4s
	69°50.6'S 032°51.8'W 68°21.0'S 043°37.3'W 69°13.1'S 038°10.2'W 1.036 15° 112° 437 km 01m51.4s
	71°02.0'S 033°00.5'W 69°35.7'S 044°07.1'W 70°25.9'S 038°30.7'W 1.036 16° 112° 434 km 01m52.2s
	72°12.0'S 033°21.0'W 70°47.5'S 044°53.8'W 71°36.8'S 039°05.5'W 1.036 16° 112° 431 km 01m52.9s
	73°20.8'S 033°54.7'W 71°56.7'S 045°58.8'W 72°45.8'S 039°56.2'W 1.037 17° 112° 428 km 01m53.5s
	74°28.4'S 034°43.6'W 73°03.2'S 047°24.2'W 73°53.1'S 041°05.1'W 1.037 17° 112° 425 km 01m53.9s
	75°34.8'S 035°50.7'W 74°07.1'S 049°12.9'W 74°58.5'S 042°35.2'W 1.037 17° 113° 422 km 01m54.2s
	76°40.0'S 037°19.9'W 75°07.9'S 051°28.7'W 76°01.8'S 044°30.6'W 1.037 17° 114° 420 km 01m54.4s
	77°43.8'S 039°16.4'W 76°05.4'S 054°15.9'W 77°02.8'S 046°56.6'W 1.037 17° 115° 418 km 01m54.4s
	78°45.7'S 041°47.3'W 76°58.6'S 057°39.7'W 78°00.8'S 049°59.5'W 1.037 17° 117° 417 km 01m54.3s
	79°45.2'S 045°02.2'W 77°46.7'S 061°45.6'W 78°55.0'S 053°47.5'W 1.037 17° 120° 415 km 01m54.1s
	80°41.3'S 049°13.4'W 78°28.4'S 066°38.4'W 79°44.3'S 058°29.4'W 1.037 17° 124° 414 km 01m53.7s
	81°32.7'S 054°36.1'W 79°01.8'S 072°20.9'W 80°26.9'S 064°14.0'W 1.037 16° 129° 414 km 01m53.2s
	82°17.2'S 061°26.9'W 79°25.1'S 078°51.6'W 81°00.8'S 071°06.7'W 1.036 16° 135° 413 km 01m52.6s
	82°51.8'S 069°57.8'W 79°35.9'S 086°02.2'W 81°23.2'S 079°05.4'W 1.036 16° 142° 413 km 01m51.8s
	83°13.0'S 080°05.3'W 79°32.0'S 093°36.6'W 81°31.4'S 087°54.6'W 1.036 15° 150° 413 km 01m50.9s
	83°17.1'S 091°17.8'W 79°11.5'S 101°12.6'W 81°23.0'S 097°05.0'W 1.036 14° 158° 414 km 01m49.7s
	83°01.5'S 102°35.4'W 78°32.6'S 108°26.8'W 80°56.3'S 105°59.0'W 1.036 13° 166° 415 km 01m48.4s
	82°25.6'S 112°55.0'W 77°33.5'S 114°59.6'W 80°10.7'S 114°03.6'W 1.035 13° 173° 416 km 01m46.8s
	81°30.4'S 121°36.4'W 76°10.4'S 120°37.9'W 79°05.6'S 120°58.3'W 1.035 11° 180° 418 km 01m45.0s
	80°16.9'S 128°29.3'W 74°14.6'S 125°14.0'W 77°39.4'S 126°35.5'W 1.034 10° 184° 420 km 01m42.8s
	78°44.8'S 133°41.4'W 71°09.1'S 128°35.3'W 75°46.7'S 130°55.2'W 1.034 8° 188° 422 km 01m40.0s
	# 76°50.5'S 137°24.7'W - - 73°09.9'S 133°54.7'W 1.033 6° 190° 424 km 01m36.2s
	67°03.1'S 138°44.1'W 67°35.9'S 129°04.9'W 67°21.9'S 134°10.7'W 1.031 0° - 418 km 01m27.8s