Quickstart

Quickstart#

Here is a quick demonstration of some of the main functionality of Invert4Geom. This assumes you are familar with Python and have successfully installed this packaged. See the tutorials for a step-by-step introduction to Invert4Geom, and the how-to guides for more in-depth guides for specific features.

Import the various modules of Invert4Geom and a few other packages

[10]:

import xarray as xr
from polartoolkit import maps
from polartoolkit import utils as polar_utils

from invert4geom import inversion, plotting, regional, synthetic, utils

Load some synthetic gravity data#

[11]:

spacing = 1000
region = (0, 40000, 0, 30000)

true_density_contrast = 2500
true_zref = 500
true_topography, _, _, grav_df = synthetic.load_synthetic_model(
    spacing=spacing,
    region=region,
    density_contrast=true_density_contrast,
    zref=true_zref,
    gravity_noise=0.2,
)
_images/quickstart_4_0.png

Create a starting model#

[12]:

# pick a density contrast (here we use the true value)
density_contrast = true_density_contrast

# pick a reference level (here we use the true value)
zref = true_zref

# create a flat surface at the reference level
starting_topography = xr.ones_like(true_topography) * zref

# prisms above zref have positive density contrast and prisms below zref have negative
# density contrast
density_grid = xr.where(
    starting_topography >= zref, density_contrast, -density_contrast
)

# create layer of prisms
starting_prisms = utils.grids_to_prisms(
    starting_topography,
    zref,
    density=density_grid,
)

Forward gravity of prism layer#

[13]:

grav_df["starting_gravity"] = starting_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=True,
)

grav_df.describe()

[13]:
northing easting upward gravity_anomaly starting_gravity
count 1271.000000 1271.00000 1271.0 1271.000000 1271.0
mean 15000.000000 20000.00000 1000.0 -0.868201 0.0
std 8947.792584 11836.81698 0.0 6.712766 0.0
min 0.000000 0.00000 1000.0 -17.194385 -0.0
25% 7000.000000 10000.00000 1000.0 -5.955907 -0.0
50% 15000.000000 20000.00000 1000.0 -1.790258 0.0
75% 23000.000000 30000.00000 1000.0 2.738776 -0.0
max 30000.000000 40000.00000 1000.0 18.024858 -0.0

Gravity misfit#

The gravity misfit is the difference between the observed gravity anomaly and the forward gravity of the starting topography.

The misfit can be separated into regional and residual components, here we assume the regional component is 0.

[14]:

grav_df = regional.regional_separation(
    method="constant",
    constant=0,
    grav_df=grav_df,
)

grav_df.describe()

[14]:
northing easting upward gravity_anomaly starting_gravity misfit reg res
count 1271.000000 1271.00000 1271.0 1271.000000 1271.0 1271.000000 1271.0 1271.000000
mean 15000.000000 20000.00000 1000.0 -0.868201 0.0 -0.868201 0.0 -0.868201
std 8947.792584 11836.81698 0.0 6.712766 0.0 6.712766 0.0 6.712766
min 0.000000 0.00000 1000.0 -17.194385 -0.0 -17.194385 0.0 -17.194385
25% 7000.000000 10000.00000 1000.0 -5.955907 -0.0 -5.955907 0.0 -5.955907
50% 15000.000000 20000.00000 1000.0 -1.790258 0.0 -1.790258 0.0 -1.790258
75% 23000.000000 30000.00000 1000.0 2.738776 -0.0 2.738776 0.0 2.738776
max 30000.000000 40000.00000 1000.0 18.024858 -0.0 18.024858 0.0 18.024858

Plot the gravity anomalies#

[15]:

grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

fig = maps.plot_grd(
    grav_grid.gravity_anomaly,
    fig_height=10,
    title="Topo free disturbance",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = maps.plot_grd(
    grav_grid.starting_gravity,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Starting gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)


fig = maps.plot_grd(
    grav_grid.misfit,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = maps.plot_grd(
    grav_grid.reg,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Regional misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = maps.plot_grd(
    grav_grid.res,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Residual misfit",
    cmap="balance+h0",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)
fig.show()

makecpt [ERROR]: Option T: min >= max
supplied min value is greater or equal to max value
Grid/points are a constant value, can't make a colorbar histogram!
makecpt [ERROR]: Option T: min >= max
supplied min value is greater or equal to max value
Grid/points are a constant value, can't make a colorbar histogram!
_images/quickstart_12_1.png

Perform inversion#

Now that we have a starting model and residual gravity misfit data we can start the inversion.

[16]:

# run the inversion
results = inversion.run_inversion(
    grav_df=grav_df,
    prism_layer=starting_prisms,
    solver_damping=0.1,
    # set stopping criteria
    max_iterations=200,
    l2_norm_tolerance=0.45,  # gravity error is .2 mGal, L2-norm is sqrt(mGal) so ~0.45
    delta_l2_norm_tolerance=1.005,
    plot_convergence=True,
)

# collect the results
topo_results, grav_results, parameters, elapsed_time = results
_images/quickstart_14_1.png 
[17]:

plotting.plot_inversion_results(
    grav_results,
    topo_results,
    parameters,
    region,
    iters_to_plot=4,
    plot_iter_results=True,
    plot_topo_results=True,
    plot_grav_results=True,
)
_images/quickstart_15_0.png
_images/quickstart_15_1.png
_images/quickstart_15_2.png 
[18]:

final_topography = topo_results.set_index(["northing", "easting"]).to_xarray().topo

_ = polar_utils.grd_compare(
    true_topography,
    final_topography,
    plot=True,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
    reverse_cpt=True,
    cmap="rain",
    diff_lims=(-20, 20),
)
_images/quickstart_16_0.png

As you can see, the inversion successfully recovered the true topography. The root mean square difference between the true and recovered topography is low, but this is not too surprising since we gave the inversion the true density contrast and reference level values.

Contents