8. Combining it all

Here we will combine the approaches of all the past inversions. Here our the main steps we will follow:

1. create a true topography and synthetic observed gravity

2. create a starting model from the points of known topography (constraints)

3. calculate a starting misfit

4. perform a damping parameter cross validation to determine the optimal damping value

5. use this damping value to perform an optimization to find the optimal density and reference level values

6. create the weighting grid to help the inversion adhere to the constraints

7. perform a full final inversion constrained by the weighting grid, using the optimal damping, density, and reference level values.

8.1. Import packages

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

import numpy as np
import pandas as pd
import polartoolkit as ptk
import verde as vd
import xarray as xr

import invert4geom

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

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

8.2. Get Input datasets

true_topography, _, _, _ = invert4geom.load_synthetic_model(
    spacing=1000,
    region=(0, 40000, 0, 30000),
)

# create lower synthetic topography data
lower_topography = invert4geom.synthetic_topography_regional(
    spacing=1000,
    region=(0, 40000, 0, 30000),
    scale=1,
    yoffset=-1000,
)

true_zref = true_topography.to_numpy().mean()
true_density_contrast = 2670 - 1  # sea level
print(f"mean of true topography: {true_zref} m ")

mean of true topography: 492.2704164812973 m

topography_model = invert4geom.create_model(
    zref=true_zref,
    density_contrast=true_density_contrast,
    topography=true_topography.to_dataset(name="upward"),
)
basement_model = invert4geom.create_model(
    zref=lower_topography.to_numpy().mean(),
    density_contrast=3300 - 2800,
    topography=lower_topography.to_dataset(name="upward"),
)

invert4geom.plot_prism_layers(
    [topography_model, basement_model],
    color_by="thickness",
    zscale=20,
)

# make pandas dataframe of locations to calculate gravity
# this represents the station locations of a gravity survey

# create an inside region to reduce gravity edge effects
# region = vd.pad_region(topography_model.region, topography_model.spacing * -5)

# create lists of coordinates
coords = vd.grid_coordinates(
    region=topography_model.region,
    spacing=topography_model.spacing,
    pixel_register=False,
    extra_coords=1001,  # survey elevation
)

# grid the coordinates
observations = vd.make_xarray_grid(
    (coords[0], coords[1]),
    data=coords[2],
    data_names="upward",
    dims=("northing", "easting"),
)

# grav_df = vd.grid_to_table(observations)

data = invert4geom.create_data(observations)

print(f"Gravity region (W,E,S,N): {data.region}")
print(f"Gravity spacing: {data.spacing} m")

data

Gravity region (W,E,S,N): (0.0, 40000.0, 0.0, 30000.0)
Gravity spacing: 1000.0 m

<xarray.Dataset> Size: 11kB
Dimensions:   (northing: 31, easting: 41)
Coordinates:
  * northing  (northing) float64 248B 0.0 1e+03 2e+03 ... 2.8e+04 2.9e+04 3e+04
  * easting   (easting) float64 328B 0.0 1e+03 2e+03 ... 3.8e+04 3.9e+04 4e+04
Data variables:
    upward    (northing, easting) float64 10kB 1.001e+03 1.001e+03 ... 1.001e+03
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       1000.0
    buffer_width:  3000.0
    inner_region:  (3000.0, 37000.0, 3000.0, 27000.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 31
  • easting: 41
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 1e+03 2e+03 ... 2.9e+04 3e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000.])

  • easting
    (easting)
    float64
    0.0 1e+03 2e+03 ... 3.9e+04 4e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000., 31000., 32000., 33000., 34000., 35000.,
           36000., 37000., 38000., 39000., 40000.])

• Data variables: (1)
  • upward
    (northing, easting)
    float64
    1.001e+03 1.001e+03 ... 1.001e+03

    array([[1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           ...,
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.]], shape=(31, 41))

• Attributes: (7)

  region :

  (0.0, 40000.0, 0.0, 30000.0)

  spacing :

  1000.0

  buffer_width :

  3000.0

  inner_region :

  (3000.0, 37000.0, 3000.0, 27000.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

data.inv.forward_gravity(topography_model, name="upper_surface_grav")
data.inv.forward_gravity(basement_model, name="lower_surface_grav")
data.inv.df.describe()

	northing	easting	upward	upper_surface_grav	lower_surface_grav
count	1271.000000	1271.00000	1271.0	1271.000000	1271.000000
mean	15000.000000	20000.00000	1001.0	-0.129619	0.004659
std	8947.792584	11836.81698	0.0	7.147140	0.830306
min	0.000000	0.00000	1001.0	-17.180401	-1.201774
25%	7000.000000	10000.00000	1001.0	-5.537833	-0.699901
50%	15000.000000	20000.00000	1001.0	-1.136301	-0.177953
75%	23000.000000	30000.00000	1001.0	3.743994	0.629891
max	30000.000000	40000.00000	1001.0	19.581120	1.860259

# contaminate gravity with 0.1 mGal of random noise
data["gravity_anomaly_no_noise"] = data.upper_surface_grav + data.lower_surface_grav
data["gravity_anomaly"], stddev = invert4geom.contaminate(
    data.gravity_anomaly_no_noise,
    stddev=0.1,
    percent=False,
    seed=0,
)
data["uncert"] = xr.full_like(data.gravity_anomaly, stddev)
data

<xarray.Dataset> Size: 62kB
Dimensions:                   (northing: 31, easting: 41)
Coordinates:
  * northing                  (northing) float64 248B 0.0 1e+03 ... 3e+04
  * easting                   (easting) float64 328B 0.0 1e+03 ... 3.9e+04 4e+04
Data variables:
    upward                    (northing, easting) float64 10kB 1.001e+03 ... ...
    upper_surface_grav        (northing, easting) float64 10kB 9.526 ... 2.855
    lower_surface_grav        (northing, easting) float64 10kB -0.04793 ... 0...
    gravity_anomaly_no_noise  (northing, easting) float64 10kB 9.478 ... 3.368
    gravity_anomaly           (northing, easting) float64 10kB 9.495 ... 3.425
    uncert                    (northing, easting) float64 10kB 0.1 0.1 ... 0.1
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       1000.0
    buffer_width:  3000.0
    inner_region:  (3000.0, 37000.0, 3000.0, 27000.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 31
  • easting: 41
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 1e+03 2e+03 ... 2.9e+04 3e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000.])

  • easting
    (easting)
    float64
    0.0 1e+03 2e+03 ... 3.9e+04 4e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000., 31000., 32000., 33000., 34000., 35000.,
           36000., 37000., 38000., 39000., 40000.])

• Data variables: (6)
  • upward
    (northing, easting)
    float64
    1.001e+03 1.001e+03 ... 1.001e+03

    array([[1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           ...,
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.],
           [1001., 1001., 1001., ..., 1001., 1001., 1001.]], shape=(31, 41))

  • upper_surface_grav
    (northing, easting)
    float64
    9.526 10.42 9.945 ... 3.298 2.855

    array([[ 9.52599368, 10.41680021,  9.94512583, ..., -5.16203178,
            -4.90305804, -4.03920067],
           [10.79727554, 11.95143824, 11.44352326, ..., -6.60127418,
            -6.26251749, -5.12226524],
           [10.78881866, 11.97658685, 11.48107397, ..., -7.30892236,
            -6.92371638, -5.65237603],
           ...,
           [17.18205561, 19.34004567, 19.11977053, ...,  2.23705848,
             2.22502307,  1.91091003],
           [17.45173684, 19.58111981, 19.33936793, ...,  3.13975322,
             3.10672962,  2.66116012],
           [15.7716986 , 17.46425943, 17.1940276 , ...,  3.33313286,
             3.29819575,  2.85527623]], shape=(31, 41))

  • lower_surface_grav
    (northing, easting)
    float64
    -0.04793 -0.06061 ... 0.6343 0.5129

    array([[-0.04793395, -0.06060717, -0.07998517, ..., -0.58762998,
            -0.52686104, -0.42462671],
           [-0.07434473, -0.09583946, -0.12326178, ..., -0.75965007,
            -0.68059402, -0.54456622],
           [-0.11580004, -0.14961758, -0.18551007, ..., -0.88593077,
            -0.79223608, -0.63211419],
           ...,
           [ 0.14108053,  0.1822873 ,  0.20599875, ...,  1.13107715,
             1.01120178,  0.80647036],
           [ 0.15041094,  0.19327919,  0.21659307, ...,  0.93672405,
             0.83967634,  0.67284909],
           [ 0.13218094,  0.16879522,  0.18876405, ...,  0.70658651,
             0.63427963,  0.51286141]], shape=(31, 41))

  • gravity_anomaly_no_noise
    (northing, easting)
    float64
    9.478 10.36 9.865 ... 3.932 3.368

    array([[ 9.47805973, 10.35619304,  9.86514066, ..., -5.74966176,
            -5.42991908, -4.46382738],
           [10.72293081, 11.85559878, 11.32026147, ..., -7.36092425,
            -6.94311151, -5.66683146],
           [10.67301862, 11.82696926, 11.2955639 , ..., -8.19485312,
            -7.71595246, -6.28449022],
           ...,
           [17.32313614, 19.52233297, 19.32576928, ...,  3.36813562,
             3.23622485,  2.71738039],
           [17.60214778, 19.77439901, 19.555961  , ...,  4.07647727,
             3.94640596,  3.33400921],
           [15.90387955, 17.63305465, 17.38279165, ...,  4.03971937,
             3.93247538,  3.36813764]], shape=(31, 41))

  • gravity_anomaly
    (northing, easting)
    float64
    9.495 10.35 9.933 ... 3.842 3.425

    array([[ 9.49467327, 10.34702308,  9.93322345, ..., -5.66722369,
            -5.27653543, -4.58569341],
           [10.87836371, 11.99422685, 11.40243314, ..., -7.15664447,
            -6.92021906, -5.72611034],
           [10.63930279, 11.72189518, 11.17183641, ..., -8.11205378,
            -7.62750407, -6.27289033],
           ...,
           [17.4219381 , 19.58892558, 19.29977445, ...,  3.36891674,
             3.25747048,  2.87259276],
           [17.55865491, 19.94461239, 19.41905626, ...,  4.1924685 ,
             3.97231975,  3.47507535],
           [15.81399477, 17.75016929, 17.47774022, ...,  4.0730052 ,
             3.84194476,  3.42512672]], shape=(31, 41))

  • uncert
    (northing, easting)
    float64
    0.1 0.1 0.1 0.1 ... 0.1 0.1 0.1 0.1

    array([[0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1],
           [0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1],
           [0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1],
           ...,
           [0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1],
           [0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1],
           [0.1, 0.1, 0.1, ..., 0.1, 0.1, 0.1]], shape=(31, 41))

• Attributes: (7)

  region :

  (0.0, 40000.0, 0.0, 30000.0)

  spacing :

  1000.0

  buffer_width :

  3000.0

  inner_region :

  (3000.0, 37000.0, 3000.0, 27000.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

fig = ptk.plot_grid(
    data.lower_surface_grav,
    fig_height=10,
    title="lower prisms gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    data.upper_surface_grav,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="upper prisms gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    data.gravity_anomaly,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Observed gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig.show()

8.3. Create “a-priori” topography measurements

# create 15 random point within the outcropping basement region
num_constraints = 15
coords = vd.scatter_points(
    region=vd.pad_region(data.region, -data.spacing),
    size=num_constraints,
    random_state=7,
)
constraint_points = pd.DataFrame(data={"easting": coords[0], "northing": coords[1]})

# sample true topography at these points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    true_topography,
    "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, "uncert"] = uncert

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

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

constraint_points.head()

	easting	northing	true_upward	upward	uncert	weight
0	3899.714996	26468.618776	655.102400	656.749723	13.102048	0.005825
1	30636.914105	6974.789900	407.011208	405.935844	8.140224	0.015091
2	17659.550795	13659.470931	470.155613	476.177579	9.403112	0.011310
3	28491.676758	27073.768551	515.239653	516.320627	10.304793	0.009417
4	38163.601456	1697.178371	419.197817	414.706788	8.383956	0.014227

8.4. Create starting model

# grid the sampled values using verde
starting_topography_kwargs = dict(
    method="splines",
    region=topography_model.region,
    spacing=topography_model.spacing,
    constraints_df=constraint_points,
    dampings=np.logspace(-40, 0, 100),
    weights_col="weight",
)

starting_topography = invert4geom.create_topography(**starting_topography_kwargs)

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

model = invert4geom.create_model(
    zref=constraint_points.upward.mean(),
    density_contrast=2500,
    topography=starting_topography,
)
model.inv.plot_model(
    color_by="density",
    zscale=20,
)

8.5. Calculate starting gravity misfit

Now we need to calculate the forward gravity of the starting topography. We then can subtract it from our observed gravity to get a starting gravity misfit.

We don’t know the optimal values for the density contrast or the reference level, so we will make a guess at appropriate values for the starting model.

data.inv.forward_gravity(model)

# calculate the true residual misfit
# true misfit is difference between noise-free gravity and starting gravity
# true regional misfit is the moho gravity
# so true residual is misfit - true_regional
data["true_res"] = (
    data.gravity_anomaly_no_noise - data.forward_gravity - data.lower_surface_grav
)
data.inv.df.describe()

	northing	easting	upward	upper_surface_grav	lower_surface_grav	gravity_anomaly_no_noise	gravity_anomaly	uncert	forward_gravity	true_res
count	1271.000000	1271.00000	1271.0	1271.000000	1271.000000	1271.000000	1271.000000	1.271000e+03	1271.000000	1271.000000
mean	15000.000000	20000.00000	1001.0	-0.129619	0.004659	-0.124959	-0.124959	1.000000e-01	1.079930	-1.209548
std	8947.792584	11836.81698	0.0	7.147140	0.830306	7.314752	7.313577	2.776650e-17	6.599125	4.456801
min	0.000000	0.00000	1001.0	-17.180401	-1.201774	-16.908338	-17.061563	1.000000e-01	-7.883841	-19.046853
25%	7000.000000	10000.00000	1001.0	-5.537833	-0.699901	-5.826630	-5.827131	1.000000e-01	-3.944309	-2.622186
50%	15000.000000	20000.00000	1001.0	-1.136301	-0.177953	-0.557795	-0.572175	1.000000e-01	-0.765707	-1.088926
75%	23000.000000	30000.00000	1001.0	3.743994	0.629891	4.171648	4.139840	1.000000e-01	4.823461	0.332004
max	30000.000000	40000.00000	1001.0	19.581120	1.860259	19.774399	19.944612	1.000000e-01	22.159800	12.693226

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

fig = ptk.plot_grid(
    data.lower_surface_grav,
    fig_height=10,
    title="True regional",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    data.true_res,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="True residual",
    hist=True,
    cbar_label="mgal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig.show()

8.6. Regional separation

def regional_comparison(data):
    # compare with true regional
    _ = ptk.grid_compare(
        data.lower_surface_grav,
        data.reg,
        robust=True,
        grid1_name="True regional field",
        grid2_name="Estimated regional field",
        hist=True,
        inset=False,
        title="difference",
        points=constraint_points,
        points_style="x.3c",
    )
    # compare with true residual
    _ = ptk.grid_compare(
        data.true_res,
        data.res,
        robust=True,
        grid1_name="True residual field",
        grid2_name="Estimated residual field",
        hist=True,
        inset=False,
        title="difference",
        points=constraint_points,
        points_style="x.3c",
    )

# estimate regional with the mean misfit at constraints
regional_grav_kwargs = dict(
    method="constraints",
    grid_method="eq_sources",
    constraints_weights_column="weight",
    constraints_df=constraint_points,
    cv=True,
    cv_kwargs=dict(
        n_trials=50,
        damping_limits=(1e-10, 1),
        # depth_limits=(1e3,100e3),
        progressbar=False,
        fname="../tmp/08_combining_it_all_regional_sep",
        # plot=True,
    ),
    depth="default",
    block_size=None,
)

data.inv.regional_separation(
    **regional_grav_kwargs,
)

regional_comparison(data)
data.inv.df.describe()

	northing	easting	upward	upper_surface_grav	lower_surface_grav	gravity_anomaly_no_noise	gravity_anomaly	uncert	forward_gravity	true_res	misfit	reg	res	starting_forward_gravity	starting_misfit	starting_reg	starting_res
count	1271.000000	1271.00000	1271.0	1271.000000	1271.000000	1271.000000	1271.000000	1.271000e+03	1271.000000	1271.000000	1271.000000	1271.000000	1271.000000	1271.000000	1271.000000	1271.000000	1271.000000
mean	15000.000000	20000.00000	1001.0	-0.129619	0.004659	-0.124959	-0.124959	1.000000e-01	1.079930	-1.209548	-1.204889	-1.050441	-0.154448	1.079930	-1.204889	-1.050441	-0.154448
std	8947.792584	11836.81698	0.0	7.147140	0.830306	7.314752	7.313577	2.776650e-17	6.599125	4.456801	4.644997	0.899712	4.423033	6.599125	4.644997	0.899712	4.423033
min	0.000000	0.00000	1001.0	-17.180401	-1.201774	-16.908338	-17.061563	1.000000e-01	-7.883841	-19.046853	-19.654062	-2.393489	-18.180347	-7.883841	-19.654062	-2.393489	-18.180347
25%	7000.000000	10000.00000	1001.0	-5.537833	-0.699901	-5.826630	-5.827131	1.000000e-01	-3.944309	-2.622186	-2.931387	-1.856200	-1.829092	-3.944309	-2.931387	-1.856200	-1.829092
50%	15000.000000	20000.00000	1001.0	-1.136301	-0.177953	-0.557795	-0.572175	1.000000e-01	-0.765707	-1.088926	-1.578812	-1.142989	-0.050248	-0.765707	-1.578812	-1.142989	-0.050248
75%	23000.000000	30000.00000	1001.0	3.743994	0.629891	4.171648	4.139840	1.000000e-01	4.823461	0.332004	0.296782	-0.286862	1.744463	4.823461	0.296782	-0.286862	1.744463
max	30000.000000	40000.00000	1001.0	19.581120	1.860259	19.774399	19.944612	1.000000e-01	22.159800	12.693226	13.114356	0.689128	14.012389	22.159800	13.114356	0.689128	14.012389

data.inv.plot_anomalies()

8.7. Single inversion

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

# setup the inversion
inv = invert4geom.Inversion(
    data,
    model,
    solver_damping=0.005,
    # set stopping criteria
    l2_norm_tolerance=0.3,  # grav error of .1mGal -> sqrt(.1) = 0.3
    delta_l2_norm_tolerance=1.008,
)

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

inv.stats_df

	iteration	rmse	l2_norm	delta_l2_norm	iter_time_sec
0	0.0	4.636082	2.153156	inf	NaN
1	1.0	1.360583	1.166440	1.845920	0.712336
2	2.0	0.478946	0.692059	1.685464	0.143883
3	3.0	0.246154	0.496139	1.394890	0.157124
4	4.0	0.213645	0.462217	1.073389	0.167241
5	5.0	0.276410	0.525747	0.879163	0.174846
6	6.0	0.432148	0.657380	0.799761	0.176366

inv.plot_inversion_results(
    iters_to_plot=2,
)

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

8.8. Damping parameter cross validation

8.8.1. Resample to create testing / training data sets

data.inv.df

	northing	easting	upward	upper_surface_grav	lower_surface_grav	gravity_anomaly_no_noise	gravity_anomaly	uncert	forward_gravity	true_res	misfit	reg	res	starting_forward_gravity	starting_misfit	starting_reg	starting_res
0	0.0	0.0	1001.0	9.525994	-0.047934	9.478060	9.494673	0.1	10.740517	-1.214524	-1.245844	-2.206690	0.960846	10.740517	-1.245844	-2.206690	0.960846
1	0.0	1000.0	1001.0	10.416800	-0.060607	10.356193	10.347023	0.1	11.470352	-1.053551	-1.123328	-2.237354	1.114026	11.470352	-1.123328	-2.237354	1.114026
2	0.0	2000.0	1001.0	9.945126	-0.079985	9.865141	9.933223	0.1	10.699837	-0.754711	-0.766613	-2.265615	1.499002	10.699837	-0.766613	-2.265615	1.499002
3	0.0	3000.0	1001.0	9.265122	-0.102636	9.162486	9.177017	0.1	9.740384	-0.475262	-0.563367	-2.291301	1.727934	9.740384	-0.563367	-2.291301	1.727934
4	0.0	4000.0	1001.0	8.528509	-0.126515	8.401994	8.352468	0.1	8.762056	-0.233547	-0.409588	-2.314268	1.904680	8.762056	-0.409588	-2.314268	1.904680
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1266	30000.0	36000.0	1001.0	3.330540	0.765378	4.095918	3.988219	0.1	5.931613	-2.601073	-1.943393	0.556764	-2.500158	5.931613	-1.943393	0.556764	-2.500158
1267	30000.0	37000.0	1001.0	3.328091	0.745878	4.073969	4.144597	0.1	6.177712	-2.849622	-2.033116	0.528620	-2.561735	6.177712	-2.033116	0.528620	-2.561735
1268	30000.0	38000.0	1001.0	3.333133	0.706587	4.039719	4.073005	0.1	6.390391	-3.057258	-2.317386	0.499154	-2.816540	6.390391	-2.317386	0.499154	-2.816540
1269	30000.0	39000.0	1001.0	3.298196	0.634280	3.932475	3.841945	0.1	6.478970	-3.180774	-2.637025	0.468717	-3.105742	6.478970	-2.637025	0.468717	-3.105742
1270	30000.0	40000.0	1001.0	2.855276	0.512861	3.368138	3.425127	0.1	5.730724	-2.875448	-2.305597	0.437636	-2.743233	5.730724	-2.305597	0.437636	-2.743233

1271 rows × 17 columns

inv.data.inv.df.columns

Index(['northing', 'easting', 'upward', 'upper_surface_grav',
       'lower_surface_grav', 'gravity_anomaly_no_noise', 'gravity_anomaly',
       'uncert', 'forward_gravity', 'true_res', 'misfit', 'reg', 'res',
       'starting_forward_gravity', 'starting_misfit', 'starting_reg',
       'starting_res', 'iter_1_initial_residual', 'iter_2_initial_residual',
       'iter_3_initial_residual', 'iter_4_initial_residual',
       'iter_5_initial_residual', 'iter_6_initial_residual'],
      dtype='str')

inv.reinitialize_inversion()

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

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

0.022462769883683984

damping_cv_obj.study.trials_dataframe()

	number	value	datetime_start	datetime_complete	duration	params_damping	user_attrs_fname	system_attrs_fixed_params	state
0	0	0.483118	2026-02-17 02:20:15.599064	2026-02-17 02:20:16.368142	0 days 00:00:00.769078	0.001000	../tmp/08_combining_it_all_damping_CV_trial_0	{'damping': 0.001}	COMPLETE
1	1	4.592549	2026-02-17 02:20:16.369218	2026-02-17 02:20:16.921685	0 days 00:00:00.552467	1.000000	../tmp/08_combining_it_all_damping_CV_trial_1	{'damping': 1}	COMPLETE
2	2	0.606439	2026-02-17 02:20:16.922760	2026-02-17 02:20:18.047964	0 days 00:00:01.125204	0.001983	../tmp/08_combining_it_all_damping_CV_trial_2	NaN	COMPLETE
3	3	0.073751	2026-02-17 02:20:18.049014	2026-02-17 02:20:19.624230	0 days 00:00:01.575216	0.062714	../tmp/08_combining_it_all_damping_CV_trial_3	NaN	COMPLETE
4	4	0.033137	2026-02-17 02:20:19.631559	2026-02-17 02:20:20.402218	0 days 00:00:00.770659	0.022463	../tmp/08_combining_it_all_damping_CV_trial_4	NaN	COMPLETE
5	5	0.051635	2026-02-17 02:20:20.403196	2026-02-17 02:20:21.298904	0 days 00:00:00.895708	0.037994	../tmp/08_combining_it_all_damping_CV_trial_5	NaN	COMPLETE
6	6	0.080030	2026-02-17 02:20:21.299891	2026-02-17 02:20:22.238563	0 days 00:00:00.938672	0.010451	../tmp/08_combining_it_all_damping_CV_trial_6	NaN	COMPLETE
7	7	0.065969	2026-02-17 02:20:22.239557	2026-02-17 02:20:22.980530	0 days 00:00:00.740973	0.016081	../tmp/08_combining_it_all_damping_CV_trial_7	NaN	COMPLETE

inv.plot_convergence()

inv.plot_inversion_results(iters_to_plot=2)

_ = ptk.grid_compare(
    true_topography,
    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",
)

# 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.upward - constraint_points.inverted_topography
)
print(f"RMSE: {rmse:.2f} m")

RMSE: 10.02 m

8.9. Z-ref and density contrast optimizations

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:

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

2. use an internal cross-validation during the optimization by separating 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 during the density contrast and zref optimization.

Well perform just the first method below. For the second method, see the following notebook:

inv.reinitialize_inversion()

# run the optimization for the zref and density
density_zref_optimization_obj = inv.optimize_inversion_zref_density_contrast(
    constraints_df=constraint_points,
    density_contrast_limits=(2300, 3000),
    zref_limits=(0, 1e3),
    n_trials=15,
    regional_grav_kwargs={
        "method": "constant",
        "constant": 0,
    },
    starting_topography=starting_topography,
    plot_scores=False,
    fname="../tmp/08_combining_it_all_density_and_zref_optimization",
)

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

# to re-load the inv object from the saved pickle file
with pathlib.Path(
    "../tmp/08_combining_it_all_density_and_zref_optimization.pickle"
).open("rb") as f:
    inv = pickle.load(f)

# to re-load the study from the saved pickle file
with pathlib.Path(
    "../tmp/08_combining_it_all_density_and_zref_optimization_study.pickle"
).open("rb") as f:
    study = pickle.load(f)

inv.solver_damping, inv.model.density_contrast, inv.model.zref

(0.022462769883683984, 2671, 487.5964339699131)

fig = invert4geom.plot_2_parameter_scores_uneven(
    study,
    param_names=(
        "params_zref",
        "params_density_contrast",
    ),
    plot_param_names=(
        "Reference level (m)",
        "Density contrast (kg/m$^3$)",
    ),
)

inv.plot_convergence()

inv.plot_inversion_results(iters_to_plot=2)

_ = ptk.grid_compare(
    true_topography,
    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",
)

8.10. Save all data

We’ll use the results from the first method of density constrast estimation since it performed better.

The notebook Uncertainty analysis will use all of the results from this notebook.

fname = "../tmp/08_combining_it_all_results.pickle"
pathlib.Path(fname).unlink(missing_ok=True)

with pathlib.Path(fname).open("ab") as file:
    pickle.dump(inv, file, protocol=pickle.HIGHEST_PROTOCOL)
    pickle.dump(regional_grav_kwargs, file, protocol=pickle.HIGHEST_PROTOCOL)
    pickle.dump(starting_topography_kwargs, file, protocol=pickle.HIGHEST_PROTOCOL)
    pickle.dump(true_topography, file, protocol=pickle.HIGHEST_PROTOCOL)
    pickle.dump(constraint_points, file, protocol=pickle.HIGHEST_PROTOCOL)