Regional field hyperparameters

The gravity misfit in Invert4Geom refers to the difference between the observed gravity anomaly and the forward calculated gravity effect of a topographic surface represented by a layer of prisms. The component of this misfit we are interested in is that which results from the difference between the true topography and the starting topographic model. This is referred to as the residual misfit, while the remaining misfit is the regional misfit. Prior to the inversion, we need to estimate this regional misfit and remove it to isolate the residual misfit.

Here we will present several methods we provide for estimating the regional component of gravity misfit. We will calculate the forward gravity effect of two topographic surfaces, the basement and Moho, from the Bishop model. Combined, these will make a synthetic observed gravity dataset. In this scenario, we are interested in inverting for the basement topography, and therefore need to estimate and remove the gravity effects of the Moho. To start, we will create a starting model for the basement topography from a series of known basement elevations. We will forward calculate the gravity of this starting model, and subtract it from the observed gravity to get a gravity misfit. From this misfit, we will then showcase several methods for estimated the regional component of the misfit.

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"] = "3857"

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

Get data

Here we will load a commonly used synthetic gravity and basement topography model. It includes topography of the Moho and the crystalline basement. The gravity effect from the Moho will represent the regional which we are aiming to isolate. We will forward model the gravity effects of both layers, and add some noise, to create an observed gravity dataset. We will then use a series of point where we know the basement topography to create a starting model, forward calculate its gravity effect, and remove it from the observed gravity to get a gravity misfit. We will then demonstrate the range of techniques implemented within Invert4Geom for isolated the regional component of this gravity misfit.

grav_ds = xr.load_dataset("../tmp/regional_sep_grav.nc")
grav_data = invert4geom.create_data(grav_ds)
grav_data.inv.df.describe()

	northing	easting	upward	basement_grav	moho_grav	gravity_anomaly_no_noise	gravity_anomaly	uncert	forward_gravity	true_res
count	7650.000000	7650.000000	7650.0	7650.000000	7650.000000	7650.000000	7650.000000	7650.0	7650.000000	7650.000000
mean	340900.000000	191900.000000	1000.0	0.302843	100.854428	101.157271	101.157271	2.0	1.548577	-1.245734
std	103923.425856	98149.168359	0.0	23.604988	27.683439	9.703771	9.919851	0.0	24.926461	9.140039
min	162900.000000	23900.000000	1000.0	-33.028677	37.680354	75.617509	70.861220	2.0	-36.414267	-37.589976
25%	250900.000000	107900.000000	1000.0	-18.377589	82.066231	96.218487	96.092524	2.0	-17.408451	-6.242037
50%	340900.000000	191900.000000	1000.0	-4.927040	109.034686	102.149812	102.057202	2.0	-3.040309	0.182452
75%	430900.000000	275900.000000	1000.0	21.431379	123.718316	106.787344	107.006336	2.0	19.823536	4.467848
max	518900.000000	359900.000000	1000.0	58.138421	139.560469	126.329916	129.736453	2.0	59.698559	19.038660

constraint_points = pd.read_csv(
    "../tmp/regional_sep_constraint_points.csv", index_col=0
)
constraint_points.describe()

	easting	northing	true_upward	upward	uncert	weight
count	15.000000	15.000000	15.000000	15.000000	15.000000	15.000000
mean	176074.067983	374901.920809	-5347.056802	-5377.150137	106.941136	0.000192
std	93521.094221	105430.197097	2146.323675	2193.119014	42.926474	0.000320
min	27437.029770	165086.148419	-8882.480862	-9295.521690	27.675267	0.000032
25%	95949.329071	307923.492200	-6589.406942	-6598.243116	77.771252	0.000058
50%	165200.779866	407905.723662	-5230.910817	-5244.731393	104.618216	0.000091
75%	265171.796613	437026.517185	-3888.562584	-3840.304455	131.788139	0.000166
max	312585.151503	514943.513495	-1383.763355	-1401.012642	177.649617	0.001306

grids = xr.load_dataset("../tmp/regional_sep_grids.nc")
basement_topo = grids.basement
moho_topo = grids.moho
starting_topography = grids.starting

fig = ptk.plot_grid(
    basement_topo,
    region=grav_data.inner_region,
    fig_height=10,
    title="Basement topography",
    reverse_cpt=True,
    cmap="rain",
    cbar_label="m",
    hist=True,
)

fig = ptk.plot_grid(
    moho_topo,
    region=grav_data.inner_region,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Moho topography",
    reverse_cpt=True,
    cmap="rain",
    cbar_label="m",
    hist=True,
)

fig = ptk.plot_grid(
    grav_data.basement_grav,
    fig=fig,
    origin_shift="both",
    xshift_amount=-1.5,
    # yshift_amount=-1.2,
    fig_height=10,
    title="Basement gravity",
    cmap="viridis",
    cbar_label="mGal",
    hist=True,
)

fig = ptk.plot_grid(
    grav_data.moho_grav,
    fig=fig,
    region=grav_data.inner_region,
    origin_shift="x",
    fig_height=10,
    title="Moho gravity",
    cmap="viridis",
    cbar_label="mGal",
    hist=True,
)

fig = ptk.plot_grid(
    grav_data.gravity_anomaly,
    fig=fig,
    region=grav_data.inner_region,
    origin_shift="x",
    fig_height=10,
    title="Observed gravity",
    cmap="viridis",
    cbar_label="mGal",
    hist=True,
)

fig.show()

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

_ = 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",
)

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

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

fig.show()

# function for plotting true vs estimated regional and residual fields
def regional_comparison(grav_data):
    # compare with true regional
    _ = ptk.grid_compare(
        grav_data.moho_grav,
        grav_data.reg,
        robust=True,
        grid1_name="True regional field",
        grid2_name="Estimated regional field",
        hist=True,
        inset=False,
        verbose="q",
        title="difference",
        points=constraint_points,
        points_style="x.3c",
    )
    # compare with true residual
    _ = ptk.grid_compare(
        grav_data.true_res,
        grav_data.res,
        robust=True,
        grid1_name="True residual field",
        grid2_name="Estimated residual field",
        hist=True,
        inset=False,
        verbose="q",
        title="difference",
        points=constraint_points,
        points_style="x.3c",
    )

Each of these regional separation techniques requires the user input to pick the value of a certain parameter, each of which has a large effect on the outcome. Below these hyperparameters for each technique are explained:

Filter technique: the width of the gaussian filter

Trend technique: the degree-order of the trend fit to the data

Equivalent Sources technique: the source depth and optionally a damping parameter and cell-size for a block reduction

Constraint Point Minimization technique: the tension factor if using minimum curvature gridding, the damping value if using bi-harmonic spline gridding, or the source depth and damping parameters if using equivalent sources gridding.

Hyperparameter estimation

In the estimating a regional field notebook we just picked these hyperparameter values arbitrarily. Below we show a more informed technique for choosing these hyperparameter values. This technique relies on having a priori information on the topography of interest to the inversion, in this case the basement topography. Therefore, we need to know the elevation of the basement at some locations to use these hyperparameter estimation techniques. These are referred to as constraints, and they would typically come from outcrops, boreholes, or other geophysical techniques such as seismic surveys. At the top of this notebook we sampled the basement topography at a series of random locations to simulate this constraints.

Recall the the gravity misfit can be divided into regional and residual components: misfit = regional + residual.

In theory, the optimal regional misfit should result in residual misfit which is 0 at constraint points (we know the topography at constraints, so the starting model is exactly equal to the true topography and therefore the residual misfit should be very small at these points). This is achieved with a regional misfit equal to the full misfit at the constraints. This can obviously be achieved by setting the regional misfit entirely equal to the total misfit (i.e. the residual misfit is 0), but this would result in no change in the topography during this inversion. Therefore, the optimal regional misfit is the model which 1) minimizes the residual misfit at constraint points, while maximizing the overall amplitude of the residual misfit.

This optimization can be performed as a multi-objective optimization (minimizing score 1 and maximizing score 2), or the scores can be combined (score1/score2) and the new resulting metric is optimized as a single-objective optimization.

We will only show this hyperparameter optimization for the constraint point minimization technique.

Constraint point minimization, uses the constraints to estimate the regional field, and therefore we can’t use the constraint points to determine the optimal hyperparameter values. For this, we must use a cross-validation scheme, where some of the constraint points are omitted, and these are used to evaluate the effectiveness of the regional estimation. This works by 1) splitting the constraints into testing and training sets, 2) using only the testing set to performing the constraint point minimization regional estimation for a range of possible hyperparameter values, and 3) assigning scores to each hyperparameter value. The hyperparameter value which produces the lowest value for the first score, and the highest values for the second score is the optimal one.

Constraint point minimization hyperparameters

Split constraint points into testing and training sets

There are a few methods we can use to split the constraints into testing and training sets:

1) Random Splitting

We can randomly split the constraints into testing and training with the function cross_validation.random_split_test_train(), but certain constraints are likely very important and greatly effect the constraint point minimization technique. This is normally not the optimal technique.

2) Leave One Out

Instead of a single split, we can perform many sets of splits, each called a Fold. A simple technique for this is to create 1 fold per number of constraint points, where each fold leaves one point out of the training set to make up the testing set. If you have many points, this results in a lot of folds which can be very computationally demanding since we will need to perform a hyperparameter cross-validation for each fold. You can use function cross_validation.split_test_train with method="LeaveOneOut" for this.

3) Block K-Fold

The third, and in most circumstances best approach is a Block K-Fold. For this the constraints are split into testing and training sets K times, and the cross-validation is performed for each of these K folds. Since spatial data is common auto-correlated (nearby points have similar values), to ensure we don’t overestimate cross-validation scores, it is best to split the data into spatial blocks, then split the data within these blocks into the training and testing sets.

block_kfold_split_df = invert4geom.split_test_train(
    constraint_points,
    method="KFold",
    shape=(20, 20),
    coord_names=("easting", "northing"),
    plot=True,
)

block_kfold_split_df

	easting	northing	true_upward	upward	uncert	weight	fold_0	fold_1	fold_2	fold_3	fold_4
0	93942.740553	165086.148419	-6166.691339	-6151.184549	123.333827	0.000066	train	train	train	test	train
1	185744.836752	437747.618234	-5230.910817	-5244.731393	104.618216	0.000091	train	train	test	train	train
2	165200.779866	503888.251841	-4551.466895	-4493.169645	91.029338	0.000121	train	train	train	test	train
3	312585.151503	412789.886738	-6178.090981	-6165.129332	123.561820	0.000065	train	train	train	train	test
4	81410.282014	268837.863004	-4562.000830	-4610.875312	91.240017	0.000120	train	test	train	train	train
5	137758.290763	436305.416136	-3739.300871	-3712.258617	74.786017	0.000179	train	train	test	train	train
6	114799.031988	407905.723662	-4037.824297	-3932.517835	80.756486	0.000153	test	train	train	train	train
7	256089.893751	300737.319016	-8556.923537	-8394.841550	171.138471	0.000034	test	train	train	train	train
8	97955.917590	381913.273733	-3696.069980	-3748.091074	73.921400	0.000183	train	test	train	train	train
9	296715.509301	315109.665385	-7771.865470	-7968.559179	155.437309	0.000041	train	train	train	train	test
10	27437.029770	370907.112964	-1383.763355	-1401.012642	27.675267	0.001306	train	train	test	train	train
11	212464.442079	413038.244716	-5835.762415	-5830.939043	116.715248	0.000073	train	test	train	train	train
12	297311.998972	202734.681191	-8882.480862	-9295.521690	177.649617	0.000032	train	train	train	test	train
13	274253.699476	491584.093601	-7000.722903	-7031.356900	140.014458	0.000051	test	train	train	train	train
14	87441.415368	514943.513495	-2611.977474	-2677.063300	52.239549	0.000366	train	train	train	train	test

Report the scores for constraint point minimization

Now that the constraint points are split into testing and training sets, we can perform a constrain point minimization regional field estimation and calculate the scores for how well it performed.

The constraint point minimization technique works by sampling the gravity misfit values at the constraints points, interpolating across the entire grid region with just the points to create the regional misfit grid, then subtracting this from the total misfit grid to get the residual misfit grid.

Once the gravity misfit values are sampled, there are 3 methods for how to perform the interpolation:

1. gridding with tensioned minimum curvature using PyGMT. This includes a hyperparameter of the tension factor (0-1).

2. gridding with bi-harmonic splines using Verde. This includes the hyperparameter of the damping factor.

3. gridding with equivalent sources using Harmonica. This includes several hyperparameters; source depth, damping factor, block size, and synthetic gravity observation elevation.

In the regional_separation_score function, the argument constraints_df is used for the constraint point minimization, while the argument testing_df is used for calculating the scores.

# values needed in all the regional separation techniques
kwargs = {
    "true_regional": grav_data.moho_grav,
    "constraints_weights_column": "weight",
}

kwargs2 = {k: v for k, v in kwargs.items() if k not in ["true_regional"]}  # noqa: FURB171

grav_data

<xarray.Dataset> Size: 491kB
Dimensions:                   (northing: 90, easting: 85)
Coordinates:
  * northing                  (northing) float64 720B 1.629e+05 ... 5.189e+05
  * easting                   (easting) float64 680B 2.39e+04 ... 3.599e+05
Data variables:
    upward                    (northing, easting) float64 61kB 1e+03 ... 1e+03
    basement_grav             (northing, easting) float64 61kB -2.757 ... -16.32
    moho_grav                 (northing, easting) float64 61kB 96.96 ... 114.7
    gravity_anomaly_no_noise  (northing, easting) float64 61kB 94.2 ... 98.37
    gravity_anomaly           (northing, easting) float64 61kB 94.46 ... 98.23
    uncert                    (northing, easting) float64 61kB 2.0 2.0 ... 2.0
    forward_gravity           (northing, easting) float64 61kB 23.65 ... -10.12
    true_res                  (northing, easting) float64 61kB -26.4 ... -6.2
Attributes:
    region:        (21900.0, 361900.0, 160900.0, 520900.0)
    spacing:       4000.0
    buffer_width:  32000.0
    inner_region:  (53900.0, 329900.0, 192900.0, 488900.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 90
  • easting: 85
• Coordinates: (2)
  • northing
    (northing)
    float64
    1.629e+05 1.669e+05 ... 5.189e+05

    array([162900., 166900., 170900., 174900., 178900., 182900., 186900., 190900.,
           194900., 198900., 202900., 206900., 210900., 214900., 218900., 222900.,
           226900., 230900., 234900., 238900., 242900., 246900., 250900., 254900.,
           258900., 262900., 266900., 270900., 274900., 278900., 282900., 286900.,
           290900., 294900., 298900., 302900., 306900., 310900., 314900., 318900.,
           322900., 326900., 330900., 334900., 338900., 342900., 346900., 350900.,
           354900., 358900., 362900., 366900., 370900., 374900., 378900., 382900.,
           386900., 390900., 394900., 398900., 402900., 406900., 410900., 414900.,
           418900., 422900., 426900., 430900., 434900., 438900., 442900., 446900.,
           450900., 454900., 458900., 462900., 466900., 470900., 474900., 478900.,
           482900., 486900., 490900., 494900., 498900., 502900., 506900., 510900.,
           514900., 518900.])

  • easting
    (easting)
    float64
    2.39e+04 2.79e+04 ... 3.599e+05

    array([ 23900.,  27900.,  31900.,  35900.,  39900.,  43900.,  47900.,  51900.,
            55900.,  59900.,  63900.,  67900.,  71900.,  75900.,  79900.,  83900.,
            87900.,  91900.,  95900.,  99900., 103900., 107900., 111900., 115900.,
           119900., 123900., 127900., 131900., 135900., 139900., 143900., 147900.,
           151900., 155900., 159900., 163900., 167900., 171900., 175900., 179900.,
           183900., 187900., 191900., 195900., 199900., 203900., 207900., 211900.,
           215900., 219900., 223900., 227900., 231900., 235900., 239900., 243900.,
           247900., 251900., 255900., 259900., 263900., 267900., 271900., 275900.,
           279900., 283900., 287900., 291900., 295900., 299900., 303900., 307900.,
           311900., 315900., 319900., 323900., 327900., 331900., 335900., 339900.,
           343900., 347900., 351900., 355900., 359900.])

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

    array([[1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           ...,
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.]], shape=(90, 85))

  • basement_grav
    (northing, easting)
    float64
    -2.757 -2.269 ... -15.08 -16.32

    array([[ -2.75728452,  -2.26912037,  -1.54295406, ..., -30.04938468,
            -29.72766243, -29.26519392],
           [ -3.19591508,  -2.84957243,  -2.26214604, ..., -30.44009774,
            -30.11138714, -29.64642291],
           [ -3.68273706,  -3.38788798,  -2.85800467, ..., -30.69676301,
            -30.37177371, -29.91630505],
           ...,
           [ 49.15593025,  49.72647897,  51.64220592, ..., -14.0350275 ,
            -15.56835694, -16.77758749],
           [ 50.12599246,  51.14187619,  54.31705249, ..., -13.65483287,
            -15.27561236, -16.55138008],
           [ 51.06923308,  52.46450044,  56.40852397, ..., -13.47186317,
            -15.07652003, -16.32296732]], shape=(90, 85))

  • moho_grav
    (northing, easting)
    float64
    96.96 97.46 98.04 ... 115.1 114.7

    array([[ 96.96009401,  97.45508342,  98.04199782, ..., 136.5651284 ,
            136.05302845, 135.25139904],
           [ 96.61602396,  97.12166895,  97.72420799, ..., 137.56380391,
            137.03881464, 136.21582586],
           [ 96.25139293,  96.76233549,  97.3746887 , ..., 138.27344225,
            137.73920133, 136.90136679],
           ...,
           [ 47.87149365,  46.33750476,  45.21040943, ..., 117.12786899,
            116.84037412, 116.42052156],
           [ 49.45769985,  47.97233442,  46.88136976, ..., 116.24281777,
            115.97890977, 115.5889103 ],
           [ 51.35664376,  49.92969985,  48.8815118 , ..., 115.29378884,
            115.05359499, 114.69434154]], shape=(90, 85))

  • gravity_anomaly_no_noise
    (northing, easting)
    float64
    94.2 95.19 96.5 ... 99.98 98.37

    array([[ 94.20280949,  95.18596305,  96.49904376, ..., 106.51574373,
            106.32536602, 105.98620512],
           [ 93.42010888,  94.27209653,  95.46206195, ..., 107.12370617,
            106.9274275 , 106.56940295],
           [ 92.56865586,  93.3744475 ,  94.51668402, ..., 107.57667924,
            107.36742763, 106.98506174],
           ...,
           [ 97.0274239 ,  96.06398373,  96.85261535, ..., 103.0928415 ,
            101.27201718,  99.64293407],
           [ 99.58369231,  99.11421061, 101.19842225, ..., 102.5879849 ,
            100.70329741,  99.03753022],
           [102.42587684, 102.39420029, 105.29003577, ..., 101.82192567,
             99.97707497,  98.37137422]], shape=(90, 85))

  • gravity_anomaly
    (northing, easting)
    float64
    94.46 94.93 97.79 ... 98.29 98.23

    array([[ 94.46018307,  94.92766646,  97.7858022 , ..., 105.76652985,
            104.14898692, 103.43675793],
           [ 94.686845  ,  95.44034129,  98.05709272, ..., 104.15547357,
            110.22801877, 108.41029206],
           [ 94.70843873,  93.4757061 ,  96.35590674, ..., 109.71524149,
            105.53007292, 108.60040874],
           ...,
           [ 96.85126287,  93.03372872,  95.55059834, ..., 103.83128619,
            100.74902571, 100.43442709],
           [ 98.74700119,  98.24134546, 102.1598126 , ..., 103.29816871,
            101.23321073,  95.94763095],
           [104.40343433, 102.72981529, 107.21056807, ..., 104.59482396,
             98.28614629,  98.23329747]], shape=(90, 85))

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

    array([[2., 2., 2., ..., 2., 2., 2.],
           [2., 2., 2., ..., 2., 2., 2.],
           [2., 2., 2., ..., 2., 2., 2.],
           ...,
           [2., 2., 2., ..., 2., 2., 2.],
           [2., 2., 2., ..., 2., 2., 2.],
           [2., 2., 2., ..., 2., 2., 2.]], shape=(90, 85))

  • forward_gravity
    (northing, easting)
    float64
    23.65 22.47 21.19 ... -10.28 -10.12

    array([[ 23.64738954,  22.4715767 ,  21.19180649, ..., -34.37783292,
            -33.98912037, -33.43151595],
           [ 24.19381782,  23.01517569,  21.72814676, ..., -34.85990868,
            -34.45410995, -33.87567805],
           [ 24.68568337,  23.50636922,  22.21525878, ..., -35.15293357,
            -34.73316668, -34.1386564 ],
           ...,
           [ 59.38519692,  58.4353564 ,  57.33507902, ..., -10.40462055,
            -10.28481513, -10.12461428],
           [ 59.15900768,  58.19207652,  57.07765541, ..., -10.44045858,
            -10.32016942, -10.15897672],
           [ 58.79544578,  57.8090524 ,  56.67965543, ..., -10.40126237,
            -10.28282641, -10.12344275]], shape=(90, 85))

  • true_res
    (northing, easting)
    float64
    -26.4 -24.74 -22.73 ... -4.794 -6.2

    array([[-26.40467406, -24.74069707, -22.73476054, ...,   4.32844824,
              4.26145794,   4.16632203],
           [-27.3897329 , -25.86474811, -23.9902928 , ...,   4.41981094,
              4.34272281,   4.22925514],
           [-28.36842043, -26.8942572 , -25.07326345, ...,   4.45617055,
              4.36139298,   4.22235135],
           ...,
           [-10.22926668,  -8.70887743,  -5.69287311, ...,  -3.63040694,
             -5.28354181,  -6.65297321],
           [ -9.03301522,  -7.05020032,  -2.76060292, ...,  -3.21437429,
             -4.95544295,  -6.39240336],
           [ -7.7262127 ,  -5.34455196,  -0.27113146, ...,  -3.0706008 ,
             -4.79369362,  -6.19952457]], shape=(90, 85))

• Attributes: (7)

  region :

  (21900.0, 361900.0, 160900.0, 520900.0)

  spacing :

  4000.0

  buffer_width :

  32000.0

  inner_region :

  (53900.0, 329900.0, 192900.0, 488900.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

# get just a single fold of the splits
training_df = block_kfold_split_df[block_kfold_split_df.fold_0 == "train"]
testing_df = block_kfold_split_df[block_kfold_split_df.fold_0 == "test"]

constraint_RMS_score, residual_RMS_score, _, _ = invert4geom.regional_separation_score(
    grav_ds=grav_data,
    method="constraints",
    grid_method="pygmt",
    tension_factor=0.3,
    constraints_df=training_df,
    testing_df=testing_df,
    **kwargs,
)
print(
    "Scores with tensioned minimum curvature gridding:\n"
    f"\tRMS of residual at constraint points: {constraint_RMS_score}\n",
    f"\tRMS of residual everywhere: {residual_RMS_score}",
)

constraint_RMS_score, residual_RMS_score, _, _ = invert4geom.regional_separation_score(
    grav_ds=grav_data,
    method="constraints",
    grid_method="verde",
    spline_dampings=0.1,
    constraints_df=training_df,
    testing_df=testing_df,
    **kwargs,
)
print(
    "Scores with bi-harmonic splines gridding:\n"
    f"\tRMS of residual at constraint points: {constraint_RMS_score}\n",
    f"\tRMS of residual everywhere: {residual_RMS_score}",
)

constraint_RMS_score, residual_RMS_score, _, _ = invert4geom.regional_separation_score(
    grav_ds=grav_data,
    method="constraints",
    grid_method="eq_sources",
    damping=10,
    block_size=10e3,
    depth=10e3,
    grav_obs_height=1e3,
    constraints_df=training_df,
    testing_df=testing_df,
    **kwargs,
)
print(
    "Scores with equivalent sources gridding:\n"
    f"\tRMS of residual at constraint points: {constraint_RMS_score}\n",
    f"\tRMS of residual everywhere: {residual_RMS_score}",
)

Scores with tensioned minimum curvature gridding:
        RMS of residual at constraint points: 5.500162994294296
        RMS of residual everywhere: 10.595938188669912
Scores with bi-harmonic splines gridding:
        RMS of residual at constraint points: 84.09071521705687
        RMS of residual everywhere: 84.72647150208118
Scores with equivalent sources gridding:
        RMS of residual at constraint points: 101.75307885138996
        RMS of residual everywhere: 104.16577380436921

Optimization for hyperparameter values

We can also run the optimization on the constraints which have been split into several folds of testing and training data. Each fold will require it’s own optimization, and the results from each fold will be merged to see which hyperparameter values are optimal across all the folds.We are just using the Block K-Folds approach here. Since we are working with synthetic data which we know the true regional field for, for each method we first optimize on minimizing the error between thetrue regional and the estimated regional. Then we perform the optimization with now knowledge of the true regional to see how the optimizations performed.

For each method below we have plot Slice Plots which show parameter values on the x axis and scores on the y axis. The scores shown on the y axis are the RMSE between the true and estimated regional fields (blue dots), the RMS of the residual at constraint points (red dots) and the RMS of the residual amplitude (greed dots).

For each method we also show create and ensemble of regional fields using a range or parameter values and calculate a cell-wise standard deviation to give an idea of the possible variation (uncertainty) each of these methods has.

Tensioned Minimum Curvature with PyGMT

Optimize on true regional

test_dfs, train_dfs = invert4geom.kfold_df_to_lists(block_kfold_split_df)

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="pygmt",
        tension_factor_limits=(0, 1),
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)

regional_comparison(grav_data)

`depth` only used if `grid_method` is 'eq_sources'

Optimize without true regional

test_dfs, train_dfs = invert4geom.kfold_df_to_lists(block_kfold_split_df)

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="pygmt",
        tension_factor_limits=(0, 1),
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        # optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)

regional_comparison(grav_data)

`depth` only used if `grid_method` is 'eq_sources'

parameter_sampling_dict = {
    "tension_factor": {
        "distribution": "uniform",
        "loc": 0,  # lower bound
        "scale": 1,  # range
    },
}
regional_misfit_stats, _ = invert4geom.regional_misfit_uncertainty(
    grav_ds=grav_data,
    runs=40,
    parameter_dict=parameter_sampling_dict,
    method="constraints",
    grid_method="pygmt",
    constraints_df=constraint_points,
    **kwargs,
)

Bi-Harmonic Splines with Verde

Optimize on true regional

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="verde",
        spline_damping_limits=(1e-15, 1),
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)
regional_comparison(grav_data)

`depth` only used if `grid_method` is 'eq_sources'

Optimize without true regional

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="verde",
        spline_damping_limits=(1e-15, 1),
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        # optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)
regional_comparison(grav_data)

`depth` only used if `grid_method` is 'eq_sources'

Alternatively, we can use an internal cross-validation (CV) which finds the parameters which create the best fit to the sampled residual values at the constraints. This is faster and simpler, as the only step that is repeated for the CV is the fitting of splines to the data, and not the entire process of sample the misfit and remove the regional from this misfit to get the residual. However, this is less informed than the above K-fold CV as it isn’t attempting to retain a high amplitude residual field, just trying to fit the data as best as possible.

grav_data.inv.regional_separation(
    method="constraints",
    grid_method="verde",
    constraints_df=constraint_points,
    spline_dampings=[*list(np.logspace(-20, 1, 100)), None],
    **kwargs2,
)
regional_comparison(grav_data)

best_trial.params["spline_dampings"]

2.6358979721550307e-09

np.log10(best_trial.params["spline_dampings"])

np.float64(-8.579071404012632)

parameter_sampling_dict = {
    "spline_dampings": {
        # "distribution": "uniform",
        # "loc": -15,  # lower bound, 1e-15
        # "scale": 5,  # range, i.e. 1e-15 to 1e-10
        "distribution": "normal",
        "loc": np.log10(best_trial.params["spline_dampings"]),  # mean
        "scale": 2,  # stdev
        "log": True,
    },
}
regional_misfit_stats, _ = invert4geom.regional_misfit_uncertainty(
    grav_ds=grav_data,
    runs=40,
    parameter_dict=parameter_sampling_dict,
    method="constraints",
    grid_method="verde",
    constraints_df=constraint_points,
    **kwargs,
)

Equivalent Sources with Harmonica

# mean_constraint_distance = np.mean(
#     vd.median_distance(
#         (constraint_points.easting, constraint_points.northing),
#         k_nearest=1,
#         ),
#     )

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="eq_sources",
        damping_limits=(1e-15, 1),
        # depth_limits=(2.5 * mean_constraint_distance, 6 *mean_constraint_distance),
        # block_size_limits=[spacing, spacing * 20],
        # grav_obs_height_limits=(1, 100e3),
        # if limits are omitted, must provide constant values
        damping=None,
        depth="default",
        block_size=None,
        # grav_obs_height=1e3,
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)

regional_comparison(grav_data)

study, grav_data, best_trial = (
    invert4geom.optimize_regional_constraint_point_minimization(
        grav_ds=grav_data,
        testing_training_df=block_kfold_split_df,
        grid_method="eq_sources",
        damping_limits=(1e-15, 1),
        # depth_limits=(2.5 * mean_constraint_distance, 6 *mean_constraint_distance),
        # block_size_limits=[spacing, spacing * 20],
        # grav_obs_height_limits=(1, 100e3),
        # if limits are omitted, must provide constant values
        damping=None,
        depth="default",
        block_size=None,
        # grav_obs_height=1e3,
        n_trials=20,
        plot=True,
        fname="../tmp/study",
        # optimize_on_true_regional_misfit=True,
        **kwargs,
    )
)

regional_comparison(grav_data)

Alternatively, we can use an internal CV which finds the parameters which create the best fit to the sampled residual values at the constraits. This is faster and simplier, as the only step that is repeated for the CV is the fitting of equivalent sources to the data, and not the entire process of sample the misfit and remove the regional from this misfit to get the residual. However, this is less informed than the above K-fold CV as it isn’t attempting to retain a high amplitude residual field, just trying to fit the data as best as possible.

grav_data.inv.regional_separation(
    method="constraints",
    grid_method="eq_sources",
    constraints_df=constraint_points,
    cv=True,
    cv_kwargs=dict(
        n_trials=40, damping_limits=(1e-15, 1), plot=True, fname="../tmp/regional_sep"
    ),
    grav_obs_height=1e3,
    depth="default",
    block_size=None,
    **kwargs2,
)

regional_comparison(grav_data)

# re-load the study from the saved pickle file
with pathlib.Path("../tmp/regional_sep.pickle").open("rb") as f:
    study = pickle.load(f)
eq_damping = min(study.best_trials, key=lambda t: t.values[0]).params["damping"]  # noqa: PD011
eq_damping

6.3567182976308735e-09

np.log10(eq_damping)

np.float64(-8.196767034203528)

mean_constraint_distance = np.mean(
    vd.median_distance(
        (constraint_points.easting, constraint_points.northing), k_nearest=1
    )
)

parameter_sampling_dict = {
    "depth": {
        "distribution": "uniform",
        "loc": 2.5 * mean_constraint_distance,  # lower bound
        "scale": 3.5 * mean_constraint_distance,  # range
    },
    "grav_obs_height": {
        "distribution": "uniform",
        "loc": 0,  # lower bound
        "scale": 10e3,  # range
    },
    # "damping": {
    #     "distribution": "uniform",
    #     "loc": -15, # lower bound
    #     "scale": 10, # range
    #     "log":True,
    # },
    "damping": {
        "distribution": "normal",
        "loc": np.log10(eq_damping),  # mean
        "scale": 2,  # stdev
        "log": True,
    },
}
regional_misfit_stats, _ = invert4geom.regional_misfit_uncertainty(
    grav_ds=grav_data,
    runs=40,
    parameter_dict=parameter_sampling_dict,
    method="constraints",
    grid_method="eq_sources",
    constraints_df=constraint_points,
    damping=eq_damping,
    depth="default",
    block_size=None,
    **kwargs,
)