Bishop Basement Model

Bishop Basement Model#

Import packages#

[1]:

# set EPSG for plotting functions
import os
import pathlib
import pickle
import string

import numpy as np
import polartoolkit as ptk
import scipy as sp
import verde as vd
import xarray as xr

import invert4geom

os.environ["POLARTOOLKIT_EPSG"] = "10598"

/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/UQpy/__init__.py:6: UserWarning:

pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.

Define model domain parameters#

To account for edge effects (decreasing gravity towards the edge of prism model), we will use a buffer region so the prism model edge is further away from the inversion domain (the gravity data extent).

[4]:

# get full region of gravity data
full_region = vd.get_region((grid.easting.values, grid.northing.values))
full_region

[4]:

(np.float64(6900.0),
 np.float64(376900.0),
 np.float64(145900.0),
 np.float64(535900.0))

[5]:

# zoom in by 10 km from each side
inversion_region = vd.pad_region(full_region, -10e3)
inversion_region

[5]:

(np.float64(16900.0),
 np.float64(366900.0),
 np.float64(155900.0),
 np.float64(525900.0))

[6]:

grav_grid = grid[["gravity_anomaly", "uncert", "upward"]]

# only retain data within the inversion region
grav_grid = grav_grid.sel(
    easting=slice(inversion_region[0], inversion_region[1]),
    northing=slice(inversion_region[2], inversion_region[3]),
)

Create “a-priori” basement measurements#

[8]:

# create grid of 36 points
coords = vd.grid_coordinates(
    region=vd.pad_region(inversion_region, -30e3),
    shape=(6, 6),
)
da = vd.make_xarray_grid(coords, data=np.ones_like(coords[0]), data_names="upward")
constraint_points = vd.grid_to_table(da).drop(columns=["upward"])

# sample true topography at these points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    grid.basement_topo,
    "true_upward",
)

constraint_points["upward"] = constraint_points.true_upward

# re-sample depths with uncertainty to emulate measurement errors
# set each points uncertainty equal to 2% of depth
uncert = np.abs(0.02 * constraint_points.upward)
constraint_points.loc[constraint_points.index, "true_uncert"] = uncert

constraint_points = invert4geom.randomly_sample_data(
    seed=0,
    data_df=constraint_points,
    data_col="upward",
    uncert_col="true_uncert",
)

# create weights column
constraint_points["weight"] = 1 / (constraint_points.true_uncert**2)

# create estimated uncertainty column
uncert = np.abs(0.02 * constraint_points.upward)
constraint_points.loc[constraint_points.index, "uncert"] = uncert

constraint_points.head()

[8]:

	northing	easting	true_upward	upward	true_uncert	weight	uncert
0	185900.0	46900.0	-6209.559082	-6193.944497	124.191182	0.000065	123.878890
1	185900.0	104900.0	-6719.067441	-6736.819871	134.381349	0.000055	134.736397
2	185900.0	162900.0	-7276.991871	-7183.784863	145.539837	0.000047	143.675697
3	185900.0	220900.0	-8420.454391	-8402.788258	168.409088	0.000035	168.055765
4	185900.0	278900.0	-8821.850945	-8916.362853	176.437019	0.000032	178.327257

Create starting basement model#

[9]:

# grid the sampled values using verde
starting_topography_kwargs = dict(
    method="splines",
    region=full_region,
    spacing=grav_data.spacing,
    constraints_df=constraint_points,
    # dampings=np.logspace(-40, 0, 200),
    dampings=None,
    # weights=constraint_points.weight,
)

starting_topography = invert4geom.create_topography(
    **starting_topography_kwargs,
)

[10]:

_ = ptk.grid_compare(
    grid.basement_topo,
    starting_topography.upward,
    region=inversion_region,
    grid1_name="True topography",
    grid2_name="Starting topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

../_images/examples_bishop_basement_model_16_0.png

[11]:

# sample the inverted topography at the constraint points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    starting_topography.upward,
    "starting_topography",
)

rmse = invert4geom.rmse(
    constraint_points.true_upward - constraint_points.starting_topography
)
print(f"RMSE: {rmse:.2f} m")

RMSE: 84.87 m

[12]:

# create model object
model = invert4geom.create_model(
    zref=constraint_points.upward.mean(),
    density_contrast=2800 - 2300,
    topography=starting_topography,
)
model.inv.plot_model(
    color_by="density",
    zscale=20,
)

../_images/examples_bishop_basement_model_18_0.png

Gravity misfit#

All inversions in Invert4Geom are based on a gravity misfit, not a gravity anomaly. This means before the inversion, we must create a starting prism model, forward model it’s gravity effect, remove it from the gravity anomaly, and get a gravity misfit.

However, if we know nothing about the starting model, it can simply be a flat layer of zero thickness, as we will use here. In this case, the forward gravity would just be zero so there is no need to perform the forward modelling. The misfit is therefore just equal to the topo-free disturbance.

Forward gravity of starting prism layer#

[13]:

# calculate forward gravity of starting prism layer
grav_data.inv.forward_gravity(model)

[14]:

_ = ptk.grid_compare(
    grav_data.gravity_anomaly,
    grav_data.forward_gravity,
    grid1_name="Observed gravity",
    grid2_name="Starting gravity",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="Gravity misfit",
    rmse_in_title=False,
    points=constraint_points,
    points_style="x.3c",
)

../_images/examples_bishop_basement_model_22_0.png

Regional estimation - constraint point minimization#

[15]:

# use the constraints to find the best regional field
regional_grav_kwargs = dict(
    method="constraints",
    grid_method="eq_sources",
    constraints_df=constraint_points,
    # constraints_weights_column="weight",
    cv=True,
    cv_kwargs=dict(
        n_trials=50,
        damping_limits=(1e-60, 10),
        depth_limits=(100, 1000e3),
        progressbar=False,
        fname="../tmp/bishop_model_regional_sep",
    ),
    # damping=None,
    # depth="default",
    block_size=None,
)

grav_data.inv.regional_separation(
    **regional_grav_kwargs,
)

grav_data.inv.df.describe()

[15]:

	northing	easting	gravity_anomaly	uncert	upward	forward_gravity	misfit	reg	res	starting_forward_gravity	starting_misfit	starting_reg	starting_res
count	1368.000000	1368.000000	1368.000000	1.368000e+03	1368.0	1368.000000	1368.000000	1368.000000	1368.000000	1368.000000	1368.000000	1368.000000	1368.000000
mean	340900.000000	191900.000000	115.021219	2.000000e-01	10.0	0.896367	114.124852	114.529292	-0.404440	0.896367	114.124852	114.529292	-0.404440
std	109698.662868	103920.936691	8.668170	8.329718e-17	0.0	40.280643	39.807129	39.976472	5.720895	40.280643	39.807129	39.976472	5.720895
min	155900.000000	16900.000000	91.608025	2.000000e-01	10.0	-59.217044	32.942960	31.710159	-21.857700	-59.217044	32.942960	31.710159	-21.857700
25%	245900.000000	104400.000000	109.730411	2.000000e-01	10.0	-27.678731	80.084974	81.831830	-2.748203	-27.678731	80.084974	81.831830	-2.748203
50%	340900.000000	191900.000000	115.379669	2.000000e-01	10.0	-8.689017	122.087862	122.950898	-0.421042	-8.689017	122.087862	122.950898	-0.421042
75%	435900.000000	279400.000000	120.316522	2.000000e-01	10.0	35.669196	144.293965	142.981497	1.877564	35.669196	144.293965	142.981497	1.877564
max	525900.000000	366900.000000	142.955701	2.000000e-01	10.0	82.655013	179.574375	177.935318	34.954172	82.655013	179.574375	177.935318	34.954172

[16]:

grav_data.inv.plot_anomalies()

../_images/examples_bishop_basement_model_25_0.png

Single inversion#

Perform a single inversion to experiment with values of stopping criteria.

[17]:

# setup the inversion
inv = invert4geom.Inversion(
    grav_data,
    model,
    solver_damping=0.01,
    # set stopping criteria
    l2_norm_tolerance=0.5,
    delta_l2_norm_tolerance=1.008,
)

[18]:

# run the inversion
inv.invert(
    plot_dynamic_convergence=True,
    results_fname="../tmp/bishop_model",
)

../_images/examples_bishop_basement_model_28_0.png

[19]:

inv.stats_df

[19]:

	iteration	rmse	l2_norm	delta_l2_norm	iter_time_sec
0	0.0	4.269698	2.066325	inf	NaN
1	1.0	1.474512	1.214295	1.701666	1.454722
2	2.0	0.721347	0.849321	1.429723	0.253134
3	3.0	0.447640	0.669059	1.269427	0.262714
4	4.0	0.325153	0.570222	1.173330	0.245501
5	5.0	0.260119	0.510019	1.118042	0.249381
6	6.0	0.220256	0.469315	1.086731	0.252319

[20]:

inv.plot_inversion_results(
    iters_to_plot=2,
)

../_images/examples_bishop_basement_model_30_0.png

../_images/examples_bishop_basement_model_30_1.png

../_images/examples_bishop_basement_model_30_2.png

[21]:

_ = ptk.grid_compare(
    grid.basement_topo,
    inv.model.topography,
    region=inversion_region,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

../_images/examples_bishop_basement_model_31_0.png

Damping parameter cross validation#

[22]:

inv.reinitialize_inversion()

# resample data at 1/2 spacing to include test points for cross-validation
inv.data = invert4geom.add_test_points(inv.data)

[23]:

damping_cv_obj = inv.optimize_inversion_damping(
    damping_limits=(0.001, 1),
    n_trials=8,
    fname="../tmp/bishop_model_damping_cv",
)
inv.solver_damping

[23]:

0.004113679284993101

../_images/examples_bishop_basement_model_34_3.png

[24]:

inv.plot_convergence()

inv.plot_inversion_results(iters_to_plot=2)

_ = ptk.grid_compare(
    grid.basement_topo,
    inv.model.topography,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

../_images/examples_bishop_basement_model_35_0.png

../_images/examples_bishop_basement_model_35_1.png

../_images/examples_bishop_basement_model_35_2.png

../_images/examples_bishop_basement_model_35_3.png

../_images/examples_bishop_basement_model_35_4.png

[25]:

# sample the inverted topography at the constraint points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    inv.model.topography,
    "inverted_topography",
)

rmse = invert4geom.rmse(
    constraint_points.true_upward - constraint_points.inverted_topography
)
print(f"RMSE: {rmse:.2f} m")

RMSE: 130.38 m

Density contrast / zref optimization#

Since this optimization uses the inversion error at constraint points, we can’t use the constraint point minimization technique to estimate the regional field since that inherently sets the inversion error at constraints to zero, invalidating the optimization scores.

There are two options for how to get around this issue:

use a different regional estimation technique while finding the optimal density contrast and zref values, then use the found optimal values with the constraint point minimization regional estimation technique afterwards.
separate the constraints into testing and training sets, so that only the training set is used during the regional separation, and only the testing set is used for scoring the density contrast and zref optimization.

Well just use the 1st option below.

Density optimization#

[26]:

inv.reinitialize_inversion()

[27]:

# run the optimization for the density contrast
density_optimization_obj = inv.optimize_inversion_zref_density_contrast(
    constraints_df=constraint_points,
    density_contrast_limits=(100, 600),
    n_trials=10,
    regional_grav_kwargs={
        "method": "constant",
        "constant": inv.data.reg.mean(),
    },
    starting_topography=starting_topography,
    plot_scores=True,
    fname="../tmp/bishop_model_density_optimization",
)

Best density_contrast value (600) is at the limit of provided values (600, 100) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
'reg' already a column of `grav_df`, but is being overwritten since calculate_regional_misfit is True

../_images/examples_bishop_basement_model_40_3.png

[28]:

# The optimal value has been saved as the Model's density contrast for future inversions
print(inv.model.density_contrast)

Zref optimization#

[29]:

inv.reinitialize_inversion()

[30]:

# run the optimization for the zref
zref_optimization_obj = inv.optimize_inversion_zref_density_contrast(
    constraints_df=constraint_points,
    zref_limits=(-8e3, 0),
    n_trials=10,
    regional_grav_kwargs={
        "method": "constant",
        "constant": inv.data.reg.mean(),
    },
    starting_topography=starting_topography,
    plot_scores=True,
    fname="../tmp/bishop_model_zref_optimization",
)

'reg' already a column of `grav_df`, but is being overwritten since calculate_regional_misfit is True

../_images/examples_bishop_basement_model_44_3.png

[31]:

# The optimal value has been saved as the model's density contrast for future inversions
print(inv.model.zref)

-6090.645590613878

Redo regional separation and inversion with optimal values#

[32]:

# recreate model with optimal zref and density contrast
model = invert4geom.create_model(
    zref=inv.model.zref,
    density_contrast=inv.model.density_contrast,
    topography=starting_topography,
)

# calculate forward gravity of starting prism layer
grav_data.inv.forward_gravity(model)

# calculate regional gravity
grav_data.inv.regional_separation(
    **regional_grav_kwargs,
)

grav_data.inv.plot_anomalies()

../_images/examples_bishop_basement_model_47_0.png

[33]:

# setup the inversion
optimal_inv = invert4geom.Inversion(
    grav_data,
    model,
    solver_damping=inv.solver_damping,
    # set stopping criteria
    l2_norm_tolerance=0.5,
    delta_l2_norm_tolerance=1.008,
)
optimal_inv.solver_damping, optimal_inv.model.zref, optimal_inv.model.density_contrast

[33]:

(0.004113679284993101, -6090.645590613878, 600)

[34]:

# run the inversion
optimal_inv.invert(
    plot_dynamic_convergence=True,
    results_fname="../tmp/bishop_model_optimal",
)

../_images/examples_bishop_basement_model_49_0.png

[35]:

_ = ptk.grid_compare(
    grid.basement_topo,
    optimal_inv.model.topography,
    region=inversion_region,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

../_images/examples_bishop_basement_model_50_0.png

[36]:

# sample the inverted topography at the constraint points
constraint_points = ptk.sample_grids(
    constraint_points,
    optimal_inv.model.topography,
    "inverted_topography",
)

rmse = ptk.rmse(constraint_points.true_upward - constraint_points.inverted_topography)
print(f"RMSE: {rmse:.2f} m")

RMSE: 121.55 m

Absolute value of inversion error#

[37]:

inversion_error = np.abs(grid.basement_topo - optimal_inv.model.topography)
fig = ptk.plot_grid(
    inversion_error,
    hist=True,
    cmap="thermal",
    title="Absolute value of inversion error",
    robust=True,
    points=constraint_points,
    points_style="x.4c",
    # points_pen="1p",
    points_fill="white",
)
fig.show()

../_images/examples_bishop_basement_model_53_0.png

[38]:

# plotting functions for uncertainty results
def uncert_plots(
    results,
    inversion_region,
    spacing,
    true_topography,
    constraints_df=None,
    weight_by=None,
):
    if (weight_by == "constraints") & (constraints_df is None):
        msg = "must provide constraints_df if weighting by constraints"
        raise ValueError(msg)

    stats_ds = invert4geom.merged_stats(
        results=results,
        plot=True,
        constraints_df=constraints_df,
        weight_by=weight_by,
        region=inversion_region,
    )

    try:
        mean = stats_ds.weighted_mean
        stdev = stats_ds.weighted_stdev
    except AttributeError:
        mean = stats_ds.z_mean
        stdev = stats_ds.z_stdev

    _ = ptk.grid_compare(
        true_topography,
        mean,
        region=vd.pad_region(inversion_region, -spacing),
        grid1_name="True topography",
        grid2_name="Inverted topography",
        robust=True,
        hist=True,
        inset=False,
        title="difference",
        reverse_cpt=True,
        cmap="rain",
        points=constraints_df,
        points_style="x.3c",
    )
    _ = ptk.grid_compare(
        np.abs(true_topography - mean),
        stdev,
        region=vd.pad_region(inversion_region, -spacing),
        grid1_name="True error",
        grid2_name="Stochastic uncertainty",
        cmap="thermal",
        robust=True,
        hist=True,
        inset=False,
        title="difference",
        points=constraints_df,
        points_style="x.3c",
        points_fill="white",
    )
    return stats_ds