# Public copy: internal compute identifiers and local environment paths were removed or generalized. Scientific values are retained. Paths may need adapting and some original run inputs are not included. """3-D gravity inversion of a subglacial trough — Manim explainer. Every curve is driven by anim_data.npz, exported from the real SimPEG forward simulation + inversion run on the the compute environment cluster. LaTeX is not functional on this platform: equations go through tex_png(). """ import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from manim import * config.background_color = "#161616" BLUE_3B = "#58C4DD" YELLOW_3B = "#FFFF00" RED_3B = "#FC6255" GREY_3B = "#888888" GREEN_3B = "#83C167" D = np.load("anim_data.npz") def ticks(ax, xvals, yvals, fs=17, color=GREY_3B, xfmt="{:.0f}", yfmt="{:.0f}"): """Manual axis tick labels. Axes(include_numbers=True) builds DecimalNumber, which renders through LaTeX -- unavailable on this platform. Plain Text mobjects instead. """ g = VGroup() for xv in xvals: g.add(Text(xfmt.format(xv), font_size=fs, color=color) .next_to(ax.c2p(xv, ax.y_range[0]), DOWN, buff=0.10)) for yv in yvals: g.add(Text(yfmt.format(yv), font_size=fs, color=color) .next_to(ax.c2p(ax.x_range[0], yv), LEFT, buff=0.10)) return g def tex_png(expr, path, color="white", fs=64, dpi=300): fig = plt.figure(figsize=(0.01, 0.01)) fig.patch.set_alpha(0) fig.text(0, 0, f"${expr}$", color=color, fontsize=fs) fig.savefig(path, dpi=dpi, transparent=True, bbox_inches="tight", pad_inches=0.05) plt.close(fig) return path class GravityInversion(Scene): def construct(self): self.title_card() self.geometry_and_forward() self.the_anomaly() self.inversion_reveal() self.what_is_resolved() self.closing() # ------------------------------------------------------------------ 1 def title_card(self): t = Text("Can gravity see the bottom of a glacier?", font_size=44, color=BLUE_3B) s = Text("A 3-D SimPEG forward model and inversion", font_size=26, color=GREY_3B).next_to(t, DOWN, buff=0.35) self.play(Write(t), run_time=1.6) self.play(FadeIn(s, shift=UP * 0.2), run_time=1.0) self.wait(1.2) self.play(FadeOut(t), FadeOut(s), run_time=0.8) # ------------------------------------------------------------------ 2 def geometry_and_forward(self): cx = D["cx"] / 1e3 bed = D["true_bed"] thk = D["true_thk"] srf = bed + thk ax = Axes( x_range=[-30, 30, 15], y_range=[-1800, 1400, 800], x_length=9.4, y_length=3.5, axis_config={"color": GREY_3B, "include_tip": False, "font_size": 20}, ).shift(DOWN * 1.05) head = Text("A glacier fills a deep bedrock trough", font_size=30, color=BLUE_3B).to_edge(UP, buff=0.35) tk = ticks(ax, [-30, -15, 0, 15, 30], [-1800, -1000, -200, 600, 1400]) xlab = Text("distance across the survey (km)", font_size=20, color=GREY_3B).next_to(tk, DOWN, buff=0.14) ylab = Text("elevation (m)", font_size=20, color=GREY_3B) ylab.rotate(PI / 2).next_to(tk, LEFT, buff=0.14) self.play(FadeIn(head, shift=DOWN * 0.2), Create(ax), FadeIn(xlab), FadeIn(ylab), FadeIn(tk), run_time=1.6) bed_line = ax.plot_line_graph(cx, bed, line_color=RED_3B, add_vertex_dots=False, stroke_width=4)["line_graph"] srf_line = ax.plot_line_graph(cx, srf, line_color=BLUE_3B, add_vertex_dots=False, stroke_width=4)["line_graph"] self.play(Create(bed_line), run_time=1.6) bed_lab = Text("bedrock", font_size=22, color=RED_3B).next_to( ax.c2p(cx[6], bed[6]), DOWN, buff=0.18) self.play(FadeIn(bed_lab), run_time=0.6) self.play(Create(srf_line), run_time=1.6) srf_lab = Text("ice surface", font_size=22, color=BLUE_3B).next_to( ax.c2p(cx[6], srf[6]), UP, buff=0.18) self.play(FadeIn(srf_lab), run_time=0.6) poly = Polygon( *[ax.c2p(x, y) for x, y in zip(cx, srf)], *[ax.c2p(x, y) for x, y in zip(cx[::-1], bed[::-1])], stroke_width=0, fill_color=BLUE_3B, fill_opacity=0.28, ) self.play(FadeIn(poly), run_time=1.0) dens = Text("ice 917 kg/m\u00b3 vs rock 2670 kg/m\u00b3", font_size=23, color=YELLOW_3B).next_to(head, DOWN, buff=0.16) self.play(FadeIn(dens, shift=UP * 0.15), run_time=1.0) self.wait(0.8) tex_png(r"\Delta\rho = -1.753\ \mathrm{g/cc}", "eq_contrast.png", color=YELLOW_3B, fs=54) eqc = ImageMobject("eq_contrast.png").scale(0.34) eqc.next_to(dens, DOWN, buff=0.12) self.play(FadeIn(eqc, shift=UP * 0.15), run_time=1.0) self.wait(1.4) self.play(FadeOut(VGroup(head, dens, bed_lab, srf_lab, xlab, ylab, tk, bed_line, srf_line, poly, ax)), FadeOut(eqc), run_time=1.0) # ------------------------------------------------------------------ 3 def the_anomaly(self): head = Text("A mass deficit makes gravity measurably weaker", font_size=32, color=BLUE_3B).to_edge(UP, buff=0.45) self.play(FadeIn(head, shift=DOWN * 0.2), run_time=1.0) rx = D["rxloc"] gx = np.unique(rx[:, 0]) ny = len(np.unique(rx[:, 1])) nx = len(gx) dcl = D["d_clean"].reshape(ny, nx) dobs = D["dobs"].reshape(ny, nx) mid = ny // 2 xs = gx / 1e3 ax = Axes( x_range=[-30, 30, 15], y_range=[70, 140, 20], x_length=9.2, y_length=3.4, axis_config={"color": GREY_3B, "include_tip": False, "font_size": 20}, ).shift(DOWN * 0.9) tk = ticks(ax, [-30, -15, 0, 15, 30], [70, 90, 110, 130]) ylab = Text("g\u1d63 anomaly (mGal)", font_size=20, color=GREY_3B) ylab.rotate(PI / 2).next_to(tk, LEFT, buff=0.14) xlab = Text("distance (km)", font_size=20, color=GREY_3B).next_to( tk, DOWN, buff=0.14) self.play(Create(ax), FadeIn(ylab), FadeIn(xlab), FadeIn(tk), run_time=1.2) clean = ax.plot_line_graph(xs, dcl[mid], line_color=BLUE_3B, add_vertex_dots=False, stroke_width=4)["line_graph"] self.play(Create(clean), run_time=2.0) note = Text("the trough shows up as a POSITIVE anomaly", font_size=24, color=YELLOW_3B).next_to(head, DOWN, buff=0.2) self.play(FadeIn(note), run_time=0.9) self.wait(0.8) # measurement noise noisy = ax.plot_line_graph(xs, dobs[mid], line_color=YELLOW_3B, add_vertex_dots=False, stroke_width=2.5)["line_graph"] nlab = Text("what an aircraft actually measures (\u00b11.0 mGal)", font_size=21, color=YELLOW_3B).next_to(note, DOWN, buff=0.16) self.play(Create(noisy), FadeIn(nlab), run_time=1.6) self.wait(1.5) self.play(FadeOut(VGroup(head, note, nlab, clean, noisy, ax, xlab, ylab, tk)), run_time=0.9) # ------------------------------------------------------------------ 4 def inversion_reveal(self): head = Text("Now invert: what density model explains that curve?", font_size=31, color=BLUE_3B).to_edge(UP, buff=0.45) self.play(FadeIn(head, shift=DOWN * 0.2), run_time=1.0) cx = D["cx"] / 1e3 cz = D["cz"] sec = D["sec"] bed = D["true_bed"] ax = Axes( x_range=[-30, 30, 15], y_range=[-1800, 1200, 1000], x_length=9.2, y_length=3.2, axis_config={"color": GREY_3B, "include_tip": False, "font_size": 19}, ).shift(DOWN * 0.75) tk = ticks(ax, [-30, -15, 0, 15, 30], [-1800, -800, 200, 1200]) self.play(Create(ax), FadeIn(tk), run_time=1.0) # Recovered density contrast. The model varies smoothly from ~-1.59 g/cc # at depth to 0 at the top, so map magnitude onto a colour ramp # (dark -> bright) rather than opacity alone: a pure-opacity map # saturates and hides exactly the vertical gradient this scene is about. dx = (cx[1] - cx[0]) dz = (cz[1] - cz[0]) ramp = color_gradient(["#0b1d2a", "#1f6f8b", BLUE_3B, "#CFF6FF"], 64) layers = [] # cz ascends (index 0 = deepest); animate from the bottom up for k in range(len(cz)): row = VGroup() for i in range(len(cx)): v = sec[k, i] if v > -0.03: continue frac = float(np.clip(v / -1.753, 0, 1)) c = ax.c2p(cx[i], cz[k]) w = abs(ax.c2p(cx[i] + dx, cz[k])[0] - c[0]) h = abs(ax.c2p(cx[i], cz[k] + dz)[1] - c[1]) row.add(Rectangle(width=w, height=h, stroke_width=0, fill_color=ramp[int(frac * 63)], fill_opacity=0.95).move_to(c)) if len(row): layers.append(row) for row in layers: self.play(FadeIn(row), run_time=0.13) self.wait(0.5) # colour key so the ramp is readable as a quantity key = VGroup() for j, f in enumerate(np.linspace(0, 1, 24)): key.add(Rectangle(width=0.13, height=0.20, stroke_width=0, fill_color=ramp[int(f * 63)], fill_opacity=0.95)) key.arrange(RIGHT, buff=0).to_edge(RIGHT, buff=0.5).shift(UP * 1.75) k0 = Text("0", font_size=17, color=GREY_3B).next_to(key, LEFT, buff=0.12) k1 = Text("-1.75 g/cc", font_size=17, color=GREY_3B).next_to(key, RIGHT, buff=0.12) self.play(FadeIn(key), FadeIn(k0), FadeIn(k1), run_time=0.7) self.key_group = VGroup(key, k0, k1) rec_lab = Text("recovered density contrast", font_size=21, color=BLUE_3B).next_to(head, DOWN, buff=0.14) self.play(FadeIn(rec_lab), run_time=0.7) self.wait(0.8) bed_line = ax.plot_line_graph(cx, bed, line_color=RED_3B, add_vertex_dots=False, stroke_width=5)["line_graph"] self.play(Create(bed_line), run_time=1.8) tl = Text("TRUE bed", font_size=22, color=RED_3B) tl.move_to(ax.c2p(cx[14], bed[14])).shift(DOWN * 0.42) self.play(FadeIn(tl), run_time=0.7) self.wait(1.6) punch = Text("the recovered mass is smeared through the column —\n" "no sharp interface appears at the true bed", font_size=23, color=YELLOW_3B, line_spacing=0.8) punch.to_edge(DOWN, buff=0.18) self.play(FadeIn(punch, shift=UP * 0.2), run_time=1.3) self.wait(2.4) self.grid = VGroup(*layers) self.play(FadeOut(VGroup(head, rec_lab, tl, punch, bed_line, ax, tk)), FadeOut(self.grid), FadeOut(self.key_group), run_time=1.0) # ------------------------------------------------------------------ 5 def what_is_resolved(self): head = Text("What the gravity data does and does not constrain", font_size=32, color=BLUE_3B).to_edge(UP, buff=0.45) self.play(FadeIn(head, shift=DOWN * 0.2), run_time=1.0) cx = D["cx"] / 1e3 thk = D["true_thk"] eq = D["eq_thk"] ax = Axes( x_range=[-30, 30, 15], y_range=[1000, 2600, 400], x_length=9.0, y_length=3.0, axis_config={"color": GREY_3B, "include_tip": False, "font_size": 19}, ).shift(DOWN * 0.55) tk = ticks(ax, [-30, -15, 0, 15, 30], [1000, 1400, 1800, 2200, 2600]) ylab = Text("equivalent ice thickness (m)", font_size=19, color=GREY_3B) ylab.rotate(PI / 2).next_to(tk, LEFT, buff=0.14) xlab = Text("distance (km)", font_size=19, color=GREY_3B).next_to( tk, DOWN, buff=0.14) self.play(Create(ax), FadeIn(ylab), FadeIn(xlab), FadeIn(tk), run_time=1.0) t_line = ax.plot_line_graph(cx, thk, line_color=RED_3B, add_vertex_dots=False, stroke_width=4)["line_graph"] r_line = ax.plot_line_graph(cx, eq, line_color=BLUE_3B, add_vertex_dots=False, stroke_width=4)["line_graph"] self.play(Create(t_line), run_time=1.4) l1 = Text("true", font_size=21, color=RED_3B).move_to( ax.c2p(cx[-9], thk[-9])).shift(UP * 0.32) self.play(FadeIn(l1), run_time=0.5) self.play(Create(r_line), run_time=1.4) l2 = Text("recovered", font_size=21, color=BLUE_3B).move_to( ax.c2p(cx[9], eq[9])).shift(DOWN * 0.36) self.play(FadeIn(l2), run_time=0.5) self.wait(1.0) good = Text("column-integrated mass: recovered to 5%", font_size=23, color=GREEN_3B) good.to_edge(DOWN, buff=0.22) self.play(FadeIn(good, shift=UP * 0.15), run_time=1.0) self.wait(1.8) self.play(FadeOut(VGroup(t_line, r_line, l1, l2, ax, ylab, xlab, tk, good)), run_time=0.8) rows = VGroup( Text("trough position 86 m error", font_size=26, color=GREEN_3B), Text("trough width +3.1%", font_size=26, color=GREEN_3B), Text("integrated mass +5.0%", font_size=26, color=GREEN_3B), Text("absolute bed depth NOT RESOLVED", font_size=26, color=RED_3B), ).arrange(DOWN, aligned_edge=LEFT, buff=0.34).shift(DOWN * 0.2) for r in rows: self.play(FadeIn(r, shift=RIGHT * 0.25), run_time=0.65) self.wait(2.2) self.play(FadeOut(rows), FadeOut(head), run_time=0.9) # ------------------------------------------------------------------ 6 def closing(self): tex_png(r"g_z \propto \int \Delta\rho\, \mathrm{d}z", "eq_final.png", color=BLUE_3B, fs=60) eq = ImageMobject("eq_final.png").scale(0.8).shift(UP * 1.1) self.play(FadeIn(eq, shift=UP * 0.2), run_time=1.2) msg = Text("Gravity measures the integral of density along a column.", font_size=27, color=WHITE).shift(DOWN * 0.25) msg2 = Text("Infinitely many depth distributions give the same anomaly.", font_size=27, color=YELLOW_3B).next_to(msg, DOWN, buff=0.3) self.play(FadeIn(msg), run_time=1.1) self.wait(0.5) self.play(FadeIn(msg2), run_time=1.1) self.wait(1.6) foot = Text("SimPEG 3-D gravity inversion \u00b7 529 stations \u00b7 49,376 cells \u00b7 \u03c7 = 1.01", font_size=21, color=GREY_3B).to_edge(DOWN, buff=0.5) self.play(FadeIn(foot), run_time=1.0) self.wait(2.2) self.play(FadeOut(VGroup(msg, msg2, foot)), FadeOut(eq), run_time=1.2)