Ctrl+K

Adhering to constraints

Adhering to constraints#

This notebook explains our approach to smallness regularization, with the goal of adhering to prior knowledge of topography we’re aiming to recover. In the scenario of a sediment-basement contact inversion, if we know the basement depth at a few points from drill holes or seismic data, we want the inverted results to adhere to these points, which we refer to as constraints.

Again, we will use the same synthetic data from the past two user guides.

Import packages#

[1]:

from __future__ import annotations

%load_ext autoreload
%autoreload 2

import logging

import pandas as pd
import verde as vd
import xarray as xr
from polartoolkit import utils as polar_utils

from invert4geom import inversion, plotting, synthetic, utils

Create observed gravity data#

To run the inversion, we need to have observed gravity data. In this simple example, we will first create a synthetic topography, which represents the true Earth topography which we hope to recover during the inverison. From this topography, we will create a layer of vertical right-rectangular prisms, which allows us to calculated the gravity effect of the topography. This will act as our observed gravity data.

True topography#

[2]:

# set grid parameters
spacing = 1000
region = [0, 40000, 0, 30000]

# create synthetic topography data
true_topography = synthetic.synthetic_topography_simple(
    spacing,
    region,
)

Starting topography#

Sample the starting topography at 10 random locations and regrid with those sampled values. This sumulates only knowing the depth to this topography at 10 boreholes.

[3]:

# create 10 random point withing the region
num_constraints = 10
coords = vd.scatter_points(region=region, 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 = utils.sample_grids(
    constraint_points, true_topography, "upward", coord_names=("easting", "northing")
)

# grid the sampled values using verde
grd = vd.Spline()
coords = (constraint_points.easting, constraint_points.northing)
grd.fit(coords, constraint_points.upward)
starting_topography = grd.grid(
    region=region,
    spacing=spacing,
).scalars

# re-sample the starting topography at the constraint points to see how the gridded did
constraint_points = utils.sample_grids(
    constraint_points,
    starting_topography,
    "starting_topography",
    coord_names=("easting", "northing"),
)

rmse = utils.rmse(constraint_points.upward - constraint_points.starting_topography)
print(f"RMSE at the constraints between the starting and true topography: {rmse:.2f} m")

RMSE at the constraints between the starting and true topography: 0.14 m

[4]:

_ = polar_utils.grd_compare(
    true_topography,
    starting_topography,
    # plot_type="xarray",
    plot=True,
    grid1_name="True topography",
    grid2_name="Starting topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
)
../_images/user_guide_adhering_to_constraints_8_0.png

Prism layer#

[5]:

# the density contrast is between rock (~2670 kg/m3) and air (~1 kg/m3)
density_contrast = 2670 - 1

# prisms are created between the mean topography value and the height of the topography
zref = true_topography.values.mean()

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

# create layer of prisms
prisms = utils.grids_to_prisms(
    true_topography,
    zref,
    density=density,
)

Forward gravity of prism layer#

[6]:

# make pandas dataframe of locations to calculate gravity
# this represents the station locations of a gravity survey
# create lists of coordinates
coords = vd.grid_coordinates(
    region=region,
    spacing=spacing,
    pixel_register=False,
    extra_coords=1000,  # survey elevation
)

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

grav_df = vd.grid_to_table(observations)

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

[6]:
northing easting upward grav
0 0.0 0.0 1000.0 9.534643
1 0.0 1000.0 1000.0 10.422834
2 0.0 2000.0 1000.0 9.949973
3 0.0 3000.0 1000.0 9.269279
4 0.0 4000.0 1000.0 8.532160
... ... ... ... ...
1266 30000.0 36000.0 1000.0 3.332716
1267 30000.0 37000.0 1000.0 3.330307
1268 30000.0 38000.0 1000.0 3.335438
1269 30000.0 39000.0 1000.0 3.300721
1270 30000.0 40000.0 1000.0 2.858299

1271 rows Ă— 4 columns

[7]:

# contaminate gravity with 0.5 mGal of random noise
grav_df["observed_grav"], stddev = synthetic.contaminate(
    grav_df.grav,
    stddev=0.5,
    percent=False,
    seed=0,
)

grav_df.set_index(["northing", "easting"]).to_xarray().observed_grav.plot()

[7]:

<matplotlib.collections.QuadMesh at 0x7f0d26922300>
../_images/user_guide_adhering_to_constraints_13_1.png

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.

[8]:

# for a reference level we use the mean value of the true topography
# zref = starting_topography.values.mean()

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

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

# calculate forward gravity of starting prism layer
grav_df["starting_grav"] = starting_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=True,
)

# calculate misfit as observed - starting
grav_df["misfit"] = grav_df["observed_grav"] - grav_df["starting_grav"]

# in many cases, we want to remove a regional signal from the misfit to isolate the
# residual signal. In this simple case, we assume there is no regional misfit and the
# full misfit is equal to the residual misfit.

# set regional misfit to 0
grav_df["reg"] = 0

# set the residual misfit to the full misfit
grav_df["res"] = grav_df.misfit

grav_df

[8]:
northing easting upward grav observed_grav starting_grav misfit reg res
0 0.0 0.0 1000.0 9.534643 9.617711 -6.928812 16.546523 0 16.546523
1 0.0 1000.0 1000.0 10.422834 10.376985 -8.207987 18.584972 0 18.584972
2 0.0 2000.0 1000.0 9.949973 10.290387 -8.333770 18.624158 0 18.624158
3 0.0 3000.0 1000.0 9.269279 9.341932 -8.213033 17.554965 0 17.554965
4 0.0 4000.0 1000.0 8.532160 8.284528 -7.985477 16.270005 0 16.270005
... ... ... ... ... ... ... ... ... ...
1266 30000.0 36000.0 1000.0 3.332716 2.794223 -0.512233 3.306456 0 3.306456
1267 30000.0 37000.0 1000.0 3.330307 3.683446 -0.787084 4.470530 0 4.470530
1268 30000.0 38000.0 1000.0 3.335438 3.501867 -1.039166 4.541033 0 4.541033
1269 30000.0 39000.0 1000.0 3.300721 2.848068 -1.241752 4.089819 0 4.089819
1270 30000.0 40000.0 1000.0 2.858299 3.143244 -1.209681 4.352925 0 4.352925

1271 rows Ă— 9 columns

Weighting grid#

To force the invesion to adhere to the starting model we need to supply a weighting grid. At each iteration, the correction grid is multiplied by this weighting grid to alter the iteration’s correction. Therefore, this weighting grid should be ~0 at the constraints, so that they aren’t altered from the starting model. These values should increase to ~1 at a distance to allow the inversion to be un-affected at locations far from constraints.

[9]:

# calculate the distance between each grid cell and the nearest constraint, then
# normalize those values between 0 and 1
min_dist = utils.normalized_mindist(
    constraint_points,
    starting_prisms,
    low=0,
    high=1,
)
starting_prisms["weights"] = min_dist
starting_prisms.weights.plot()

[9]:

<matplotlib.collections.QuadMesh at 0x7f0d1fcaaae0>
../_images/user_guide_adhering_to_constraints_17_1.png

Perform inversion#

Now we can perform the inversion, supplied the argument weights_after_solving=True and ensuring that the weighting grid is included as the weights variable to the argument prism_layer. Note that we have increased the max_iterations from 10 to 20. This is because the weighting grid reduces the correction values at each iterations, resulting in the need for more iterations.

[11]:

# set Python's logging level to get information about the inversion\s progress
logger = logging.getLogger()
logger.setLevel(logging.INFO)

# run the inversion
results = inversion.run_inversion(
    input_grav=grav_df,
    input_grav_column="observed_grav",
    prism_layer=starting_prisms,
    zref=zref,
    density_contrast=density_contrast,
    # display the convergence of the inversion
    plot_convergence=True,
    # choose the small prism approximation method for calculating the vertical
    # derivative of gravity
    deriv_type="annulus",
    solver_damping=0.1,
    # set stopping criteria
    max_iterations=100,
    l2_norm_tolerance=0.5,
    delta_l2_norm_tolerance=1.005,
    # enable the use of weights
    weights_after_solving=True,
)

# collect the results
topo_results, grav_results, parameters, elapsed_time = results

INFO:root:starting inversion
INFO:root:extracted prism spacing is 1000.0

INFO:root:
 ####################################
 iteration 1
INFO:root:Layer correction median: -5.068 m, RMSE:33.278 m
INFO:root:updated misfit RMSE: 5.4494
INFO:root:updated L2-norm: 2.3344, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.141, tolerance: 1.005
INFO:root:
 ####################################
 iteration 2
INFO:root:Layer correction median: -3.0173 m, RMSE:25.4069 m
INFO:root:updated misfit RMSE: 4.297
INFO:root:updated L2-norm: 2.0729, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1261, tolerance: 1.005
INFO:root:
 ####################################
 iteration 3
INFO:root:Layer correction median: 0.0667 m, RMSE:19.8508 m
INFO:root:updated misfit RMSE: 3.4811
INFO:root:updated L2-norm: 1.8658, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.111, tolerance: 1.005
INFO:root:
 ####################################
 iteration 4
INFO:root:Layer correction median: 1.5637 m, RMSE:15.9003 m
INFO:root:updated misfit RMSE: 2.891
INFO:root:updated L2-norm: 1.7003, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0973, tolerance: 1.005
INFO:root:
 ####################################
 iteration 5
INFO:root:Layer correction median: 1.9424 m, RMSE:13.04 m
INFO:root:updated misfit RMSE: 2.4529
INFO:root:updated L2-norm: 1.5662, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0856, tolerance: 1.005
INFO:root:
 ####################################
 iteration 6
INFO:root:Layer correction median: 2.0415 m, RMSE:10.9158 m
INFO:root:updated misfit RMSE: 2.1194
INFO:root:updated L2-norm: 1.4558, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0758, tolerance: 1.005
INFO:root:
 ####################################
 iteration 7
INFO:root:Layer correction median: 1.8222 m, RMSE:9.2986 m
INFO:root:updated misfit RMSE: 1.8594
INFO:root:updated L2-norm: 1.3636, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0676, tolerance: 1.005
INFO:root:
 ####################################
 iteration 8
INFO:root:Layer correction median: 1.5765 m, RMSE:8.0393 m
INFO:root:updated misfit RMSE: 1.6527
INFO:root:updated L2-norm: 1.2856, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0607, tolerance: 1.005
INFO:root:
 ####################################
 iteration 9
INFO:root:Layer correction median: 1.2876 m, RMSE:7.0392 m
INFO:root:updated misfit RMSE: 1.4854
INFO:root:updated L2-norm: 1.2188, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0548, tolerance: 1.005
INFO:root:
 ####################################
 iteration 10
INFO:root:Layer correction median: 1.0493 m, RMSE:6.2312 m
INFO:root:updated misfit RMSE: 1.3479
INFO:root:updated L2-norm: 1.161, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0498, tolerance: 1.005
INFO:root:
 ####################################
 iteration 11
INFO:root:Layer correction median: 0.8518 m, RMSE:5.5684 m
INFO:root:updated misfit RMSE: 1.2333
INFO:root:updated L2-norm: 1.1105, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0454, tolerance: 1.005
INFO:root:
 ####################################
 iteration 12
INFO:root:Layer correction median: 0.7535 m, RMSE:5.0174 m
INFO:root:updated misfit RMSE: 1.1365
INFO:root:updated L2-norm: 1.0661, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0417, tolerance: 1.005
INFO:root:
 ####################################
 iteration 13
INFO:root:Layer correction median: 0.6356 m, RMSE:4.5537 m
INFO:root:updated misfit RMSE: 1.0539
INFO:root:updated L2-norm: 1.0266, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0384, tolerance: 1.005
INFO:root:
 ####################################
 iteration 14
INFO:root:Layer correction median: 0.5434 m, RMSE:4.1591 m
INFO:root:updated misfit RMSE: 0.9827
INFO:root:updated L2-norm: 0.9913, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0356, tolerance: 1.005
INFO:root:
 ####################################
 iteration 15
INFO:root:Layer correction median: 0.4579 m, RMSE:3.8198 m
INFO:root:updated misfit RMSE: 0.9208
INFO:root:updated L2-norm: 0.9596, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0331, tolerance: 1.005
INFO:root:
 ####################################
 iteration 16
INFO:root:Layer correction median: 0.4052 m, RMSE:3.5253 m
INFO:root:updated misfit RMSE: 0.8664
INFO:root:updated L2-norm: 0.9308, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0309, tolerance: 1.005
INFO:root:
 ####################################
 iteration 17
INFO:root:Layer correction median: 0.3523 m, RMSE:3.2675 m
INFO:root:updated misfit RMSE: 0.8183
INFO:root:updated L2-norm: 0.9046, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.029, tolerance: 1.005
INFO:root:
 ####################################
 iteration 18
INFO:root:Layer correction median: 0.3064 m, RMSE:3.0401 m
INFO:root:updated misfit RMSE: 0.7755
INFO:root:updated L2-norm: 0.8806, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0272, tolerance: 1.005
INFO:root:
 ####################################
 iteration 19
INFO:root:Layer correction median: 0.2735 m, RMSE:2.8381 m
INFO:root:updated misfit RMSE: 0.7372
INFO:root:updated L2-norm: 0.8586, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0257, tolerance: 1.005
INFO:root:
 ####################################
 iteration 20
INFO:root:Layer correction median: 0.2482 m, RMSE:2.6574 m
INFO:root:updated misfit RMSE: 0.7027
INFO:root:updated L2-norm: 0.8382, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0243, tolerance: 1.005
INFO:root:
 ####################################
 iteration 21
INFO:root:Layer correction median: 0.2348 m, RMSE:2.4949 m
INFO:root:updated misfit RMSE: 0.6714
INFO:root:updated L2-norm: 0.8194, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.023, tolerance: 1.005
INFO:root:
 ####################################
 iteration 22
INFO:root:Layer correction median: 0.2075 m, RMSE:2.3479 m
INFO:root:updated misfit RMSE: 0.643
INFO:root:updated L2-norm: 0.8019, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0218, tolerance: 1.005
INFO:root:
 ####################################
 iteration 23
INFO:root:Layer correction median: 0.188 m, RMSE:2.2145 m
INFO:root:updated misfit RMSE: 0.6172
INFO:root:updated L2-norm: 0.7856, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0208, tolerance: 1.005
INFO:root:
 ####################################
 iteration 24
INFO:root:Layer correction median: 0.1646 m, RMSE:2.0928 m
INFO:root:updated misfit RMSE: 0.5934
INFO:root:updated L2-norm: 0.7704, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0198, tolerance: 1.005
INFO:root:
 ####################################
 iteration 25
INFO:root:Layer correction median: 0.1533 m, RMSE:1.9814 m
INFO:root:updated misfit RMSE: 0.5717
INFO:root:updated L2-norm: 0.7561, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0189, tolerance: 1.005
INFO:root:
 ####################################
 iteration 26
INFO:root:Layer correction median: 0.1431 m, RMSE:1.879 m
INFO:root:updated misfit RMSE: 0.5516
INFO:root:updated L2-norm: 0.7427, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.018, tolerance: 1.005
INFO:root:
 ####################################
 iteration 27
INFO:root:Layer correction median: 0.1292 m, RMSE:1.7847 m
INFO:root:updated misfit RMSE: 0.5331
INFO:root:updated L2-norm: 0.7302, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0172, tolerance: 1.005
INFO:root:
 ####################################
 iteration 28
INFO:root:Layer correction median: 0.1249 m, RMSE:1.6976 m
INFO:root:updated misfit RMSE: 0.516
INFO:root:updated L2-norm: 0.7183, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0165, tolerance: 1.005
INFO:root:
 ####################################
 iteration 29
INFO:root:Layer correction median: 0.1145 m, RMSE:1.6169 m
INFO:root:updated misfit RMSE: 0.5001
INFO:root:updated L2-norm: 0.7072, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0158, tolerance: 1.005
INFO:root:
 ####################################
 iteration 30
INFO:root:Layer correction median: 0.1095 m, RMSE:1.5421 m
INFO:root:updated misfit RMSE: 0.4853
INFO:root:updated L2-norm: 0.6967, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0151, tolerance: 1.005
INFO:root:
 ####################################
 iteration 31
INFO:root:Layer correction median: 0.1047 m, RMSE:1.4724 m
INFO:root:updated misfit RMSE: 0.4716
INFO:root:updated L2-norm: 0.6867, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0145, tolerance: 1.005
INFO:root:
 ####################################
 iteration 32
INFO:root:Layer correction median: 0.101 m, RMSE:1.4076 m
INFO:root:updated misfit RMSE: 0.4587
INFO:root:updated L2-norm: 0.6773, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0139, tolerance: 1.005
INFO:root:
 ####################################
 iteration 33
INFO:root:Layer correction median: 0.0959 m, RMSE:1.347 m
INFO:root:updated misfit RMSE: 0.4467
INFO:root:updated L2-norm: 0.6684, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0133, tolerance: 1.005
INFO:root:
 ####################################
 iteration 34
INFO:root:Layer correction median: 0.0837 m, RMSE:1.2904 m
INFO:root:updated misfit RMSE: 0.4355
INFO:root:updated L2-norm: 0.6599, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0128, tolerance: 1.005
INFO:root:
 ####################################
 iteration 35
INFO:root:Layer correction median: 0.074 m, RMSE:1.2375 m
INFO:root:updated misfit RMSE: 0.425
INFO:root:updated L2-norm: 0.6519, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0123, tolerance: 1.005
INFO:root:
 ####################################
 iteration 36
INFO:root:Layer correction median: 0.0651 m, RMSE:1.1878 m
INFO:root:updated misfit RMSE: 0.4151
INFO:root:updated L2-norm: 0.6443, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0118, tolerance: 1.005
INFO:root:
 ####################################
 iteration 37
INFO:root:Layer correction median: 0.0572 m, RMSE:1.1413 m
INFO:root:updated misfit RMSE: 0.4058
INFO:root:updated L2-norm: 0.637, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0114, tolerance: 1.005
INFO:root:
 ####################################
 iteration 38
INFO:root:Layer correction median: 0.0542 m, RMSE:1.0975 m
INFO:root:updated misfit RMSE: 0.3971
INFO:root:updated L2-norm: 0.6301, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.011, tolerance: 1.005
INFO:root:
 ####################################
 iteration 39
INFO:root:Layer correction median: 0.0506 m, RMSE:1.0563 m
INFO:root:updated misfit RMSE: 0.3888
INFO:root:updated L2-norm: 0.6235, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0106, tolerance: 1.005
INFO:root:
 ####################################
 iteration 40
INFO:root:Layer correction median: 0.0505 m, RMSE:1.0176 m
INFO:root:updated misfit RMSE: 0.381
INFO:root:updated L2-norm: 0.6173, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0102, tolerance: 1.005
INFO:root:
 ####################################
 iteration 41
INFO:root:Layer correction median: 0.0432 m, RMSE:0.981 m
INFO:root:updated misfit RMSE: 0.3737
INFO:root:updated L2-norm: 0.6113, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0098, tolerance: 1.005
INFO:root:
 ####################################
 iteration 42
INFO:root:Layer correction median: 0.0388 m, RMSE:0.9466 m
INFO:root:updated misfit RMSE: 0.3667
INFO:root:updated L2-norm: 0.6056, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0095, tolerance: 1.005
INFO:root:
 ####################################
 iteration 43
INFO:root:Layer correction median: 0.0383 m, RMSE:0.9141 m
INFO:root:updated misfit RMSE: 0.3601
INFO:root:updated L2-norm: 0.6001, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0091, tolerance: 1.005
INFO:root:
 ####################################
 iteration 44
INFO:root:Layer correction median: 0.0356 m, RMSE:0.8833 m
INFO:root:updated misfit RMSE: 0.3538
INFO:root:updated L2-norm: 0.5948, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0088, tolerance: 1.005
INFO:root:
 ####################################
 iteration 45
INFO:root:Layer correction median: 0.0381 m, RMSE:0.8542 m
INFO:root:updated misfit RMSE: 0.3479
INFO:root:updated L2-norm: 0.5898, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0085, tolerance: 1.005
INFO:root:
 ####################################
 iteration 46
INFO:root:Layer correction median: 0.0365 m, RMSE:0.8267 m
INFO:root:updated misfit RMSE: 0.3422
INFO:root:updated L2-norm: 0.585, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0082, tolerance: 1.005
INFO:root:
 ####################################
 iteration 47
INFO:root:Layer correction median: 0.0351 m, RMSE:0.8007 m
INFO:root:updated misfit RMSE: 0.3368
INFO:root:updated L2-norm: 0.5804, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.008, tolerance: 1.005
INFO:root:
 ####################################
 iteration 48
INFO:root:Layer correction median: 0.0346 m, RMSE:0.776 m
INFO:root:updated misfit RMSE: 0.3317
INFO:root:updated L2-norm: 0.5759, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0077, tolerance: 1.005
INFO:root:
 ####################################
 iteration 49
INFO:root:Layer correction median: 0.0335 m, RMSE:0.7525 m
INFO:root:updated misfit RMSE: 0.3268
INFO:root:updated L2-norm: 0.5717, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0075, tolerance: 1.005
INFO:root:
 ####################################
 iteration 50
INFO:root:Layer correction median: 0.0323 m, RMSE:0.7303 m
INFO:root:updated misfit RMSE: 0.3221
INFO:root:updated L2-norm: 0.5675, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0072, tolerance: 1.005
INFO:root:
 ####################################
 iteration 51
INFO:root:Layer correction median: 0.0323 m, RMSE:0.7091 m
INFO:root:updated misfit RMSE: 0.3176
INFO:root:updated L2-norm: 0.5636, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.007, tolerance: 1.005
INFO:root:
 ####################################
 iteration 52
INFO:root:Layer correction median: 0.0299 m, RMSE:0.689 m
INFO:root:updated misfit RMSE: 0.3133
INFO:root:updated L2-norm: 0.5598, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0068, tolerance: 1.005
INFO:root:
 ####################################
 iteration 53
INFO:root:Layer correction median: 0.0292 m, RMSE:0.6699 m
INFO:root:updated misfit RMSE: 0.3092
INFO:root:updated L2-norm: 0.5561, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0066, tolerance: 1.005
INFO:root:
 ####################################
 iteration 54
INFO:root:Layer correction median: 0.0285 m, RMSE:0.6517 m
INFO:root:updated misfit RMSE: 0.3053
INFO:root:updated L2-norm: 0.5525, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0064, tolerance: 1.005
INFO:root:
 ####################################
 iteration 55
INFO:root:Layer correction median: 0.0274 m, RMSE:0.6344 m
INFO:root:updated misfit RMSE: 0.3015
INFO:root:updated L2-norm: 0.5491, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0063, tolerance: 1.005
INFO:root:
 ####################################
 iteration 56
INFO:root:Layer correction median: 0.0282 m, RMSE:0.6179 m
INFO:root:updated misfit RMSE: 0.2979
INFO:root:updated L2-norm: 0.5458, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0061, tolerance: 1.005
INFO:root:
 ####################################
 iteration 57
INFO:root:Layer correction median: 0.0266 m, RMSE:0.6022 m
INFO:root:updated misfit RMSE: 0.2944
INFO:root:updated L2-norm: 0.5426, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0059, tolerance: 1.005
INFO:root:
 ####################################
 iteration 58
INFO:root:Layer correction median: 0.026 m, RMSE:0.5871 m
INFO:root:updated misfit RMSE: 0.291
INFO:root:updated L2-norm: 0.5395, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0058, tolerance: 1.005
INFO:root:
 ####################################
 iteration 59
INFO:root:Layer correction median: 0.0278 m, RMSE:0.5728 m
INFO:root:updated misfit RMSE: 0.2878
INFO:root:updated L2-norm: 0.5364, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0056, tolerance: 1.005
INFO:root:
 ####################################
 iteration 60
INFO:root:Layer correction median: 0.0239 m, RMSE:0.5591 m
INFO:root:updated misfit RMSE: 0.2846
INFO:root:updated L2-norm: 0.5335, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0055, tolerance: 1.005
INFO:root:
 ####################################
 iteration 61
INFO:root:Layer correction median: 0.0248 m, RMSE:0.546 m
INFO:root:updated misfit RMSE: 0.2816
INFO:root:updated L2-norm: 0.5307, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0054, tolerance: 1.005
INFO:root:
 ####################################
 iteration 62
INFO:root:Layer correction median: 0.0234 m, RMSE:0.5334 m
INFO:root:updated misfit RMSE: 0.2787
INFO:root:updated L2-norm: 0.5279, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0052, tolerance: 1.005
INFO:root:
 ####################################
 iteration 63
INFO:root:Layer correction median: 0.023 m, RMSE:0.5214 m
INFO:root:updated misfit RMSE: 0.2759
INFO:root:updated L2-norm: 0.5252, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0051, tolerance: 1.005
INFO:root:
 ####################################
 iteration 64
INFO:root:Layer correction median: 0.0221 m, RMSE:0.5099 m
INFO:root:updated misfit RMSE: 0.2731
INFO:root:updated L2-norm: 0.5226, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.005, tolerance: 1.005
INFO:root:
 ####################################
 iteration 65
INFO:root:Layer correction median: 0.0223 m, RMSE:0.4989 m
INFO:root:updated misfit RMSE: 0.2705
INFO:root:updated L2-norm: 0.5201, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0049, tolerance: 1.005
INFO:root:
Inversion terminated after 65 iterations because there was no significant variation in the L2-norm over 2 iterations
Change parameter 'delta_l2_norm_tolerance' if desired.
../_images/user_guide_adhering_to_constraints_19_3.png 
[12]:

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/user_guide_adhering_to_constraints_20_0.png
../_images/user_guide_adhering_to_constraints_20_1.png
../_images/user_guide_adhering_to_constraints_20_2.png 
[13]:

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

_ = polar_utils.grd_compare(
    true_topography,
    final_topography,
    # plot_type="xarray",
    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),
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
)
../_images/user_guide_adhering_to_constraints_21_0.png

As we can see by the low error values surround the constraints (black crosses), include the weighting grid has help the inversion adhere to the constraints. We can sample the inverted topography at the constraints and compare with the constraints true values.

[14]:

# sample the inverted topography at the constraint points
constraint_points = utils.sample_grids(
    constraint_points,
    final_topography,
    "inverted_topography",
    coord_names=("easting", "northing"),
)

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

RMSE: 1.73 m

While this inversion’s results are very similar as the past example (including_starting_model.ipynb), the RMSE at the constraint points is much more, showing the inversion is correctly adhering to these constraints.

Note at we increaded the max iterations allowed from 30 to 100. This is necessary since with the weighting grid we require many more iterations to get to a similar outcome.

Contents