Instantly share code, notes, and snippets.

# seba-perez/eclipse_map.py

Last active November 29, 2021 14:30
Script to generate map of the 2021 Total Solar Eclipse
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
 # 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
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
 # 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

### seba-perez commented Nov 29, 2021

my-timelapse-24fps.mp4