Reference level cross validation

In scenarios where you have no prior knowledge of the elevation of the density contrast of interest, we use a flat starting model with an arbitrary reference level. The reference level (zref) is very important to the inversion as different values can vertically shift the inverted topography. Here we present a cross-validation approach to determine the optimal value for zref. This follows the same approach as the past notebook (density_cross_validation.ipynb). To simplify, we assume we know the appropiate density constrast value to use.

Import packages

from __future__ import annotations

%load_ext autoreload
%autoreload 2


import logging

import numpy as np
import pandas as pd
import verde as vd
import xarray as xr
from polartoolkit import maps
from polartoolkit import utils as polar_utils
from tqdm.autonotebook import tqdm

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

Create observed gravity data

True topography

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

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

Constraint points

Sample the starting topography at 10 random locations to create a set of constraints points, simulating locations where the topography is known.

# 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")
)

# plot the true topography
fig = maps.plot_grd(
    true_topography,
    fig_height=10,
    title="True topography",
    cmap="rain",
    reverse_cpt=True,
    grd2_cpt=True,
    cbar_label="elevation (m)",
    frame=["nSWe", "xaf10000", "yaf10000"],
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
)
fig.show()

True topography prism layer

# 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()
print(f"mean of true topography: {zref} m ")

# 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,
)

mean of true topography: 492.2704164812973 m

Forward gravity of prism layer

# 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

	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

# 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()

<matplotlib.collections.QuadMesh at 0x7f65584551c0>

For simplicity here we assume that we know the optimal density contrast value and use this when creating our starting model.

What we don’t know in this scenario is what reference elevation zref to use. We will use a cross-validation of a range of elevations to find the one which is optimal.

Each zref will give a cross validation score, and the lowest score will show which zref is optimal.

The cross validation score is calculated as the root mean square error between the points of known topography (constraints) and the inverted topography at those points, same as in the past nootebook.

Get Cross Validation Score

First we need to use the zref value to create the starting model. We will then use this in an inversion to calculate a score.

# set a zref value
zref = 300  # instead of the correct value of ~490 m

# create flat topography grid with a constant height
starting_topography = xr.full_like(true_topography, zref)

# 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,
)

# gravity of starting model is 0 since its flat, so observed_grav = misfit
grav_df["misfit"] = grav_df["observed_grav"]
grav_df["reg"] = 0
grav_df["res"] = grav_df.misfit

# set kwargs to pass to the inversion
kwargs = {
    "input_grav_column": "observed_grav",
    "prism_layer": starting_prisms,
    "deriv_type": "annulus",
    "solver_damping": 0.1,
    "zref": zref,
    "density_contrast": density_contrast,
    # set stopping criteria
    "max_iterations": 30,
    "l2_norm_tolerance": 0.5,
    "delta_l2_norm_tolerance": 1.005,
}

# run inversion, calculate the score
score = cross_validation.constraints_cv_score(
    grav=grav_df,
    constraints=constraint_points,
    **kwargs,
)
score

194.8080323346519

Cross Validation

Lets see if we can improve the score with other values for zref.

Now we can repeat this with a range of zref values to find the optimal (lowest) score. But remember we need to recreate the starting model with each zref.

# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.WARNING)

# set which zref values to include
zrefs = np.linspace(400, 600, 8)

# run inversions and collect scores
scores = []
for zref in tqdm(zrefs, desc="Reference levels"):
    # create flat topography grid with a constant height
    starting_topography = xr.full_like(true_topography, zref)

    # re-calculate density grid with new 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
    # this isn't exactly necessary since the prism's thickness is 0 so gravity will be 0
    grav_df["starting_grav"] = starting_prisms.prism_layer.gravity(
        coordinates=(
            grav_df.easting,
            grav_df.northing,
            grav_df.upward,
        ),
        field="g_z",
        progressbar=False,
    )

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

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

    # remove the regional from the misfit to get the residual
    grav_df["res"] = grav_df.misfit - grav_df.reg

    # update zref value in kwargs
    kwargs["zref"] = zref

    # update starting model in kwargs
    kwargs["prism_layer"] = starting_prisms

    # run cross validation
    score = cross_validation.constraints_cv_score(
        grav=grav_df,
        constraints=constraint_points,
        **kwargs,
    )
    scores.append(score)

# print zref and score pairs
for zref, score in zip(zrefs, scores):
    print(f"Reference level: {zref} -> Score: {score}")

best_idx = np.argmin(scores)
best_score = scores[best_idx]
best_zref = zrefs[best_idx]
print(f"Best score of {best_score} with reference level={best_zref}")

Reference level: 400.0 -> Score: 95.28058460003426
Reference level: 428.57142857142856 -> Score: 66.79420158764773
Reference level: 457.14285714285717 -> Score: 38.39424789324149
Reference level: 485.7142857142857 -> Score: 10.840430866668894
Reference level: 514.2857142857143 -> Score: 19.748427566026876
Reference level: 542.8571428571429 -> Score: 47.84885673118
Reference level: 571.4285714285714 -> Score: 76.26691479957566
Reference level: 600.0 -> Score: 104.7602426080746
Best score of 10.840430866668894 with reference level=485.7142857142857

plotting.plot_cv_scores(
    scores,
    zrefs,
    param_name="Reference level",
    # logx=True,
    # logy=True,
)

Run inversion with optimal value

# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.INFO)

# use the best zref value
zref = best_zref

# create flat topography grid with a constant height
starting_topography = xr.full_like(true_topography, zref)

# re-calculate density grid with the best zref
density = xr.where(
    starting_topography >= zref,
    kwargs.get("density_contrast"),
    -kwargs.get("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=False,
)

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

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

# remove the regional from the misfit to get the residual
grav_df["res"] = grav_df.misfit - grav_df.reg

# update starting model in kwargs
kwargs["prism_layer"] = starting_prisms

# make new kwargs without zref or density
new_kwargs = {
    key: value
    for key, value in kwargs.items()
    if key
    not in [
        "zref",
        "density_contrast",
    ]
}
results = inversion.run_inversion(
    input_grav=grav_df,
    plot_convergence=True,
    zref=zref,
    density_contrast=kwargs.get("density_contrast"),
    **new_kwargs,
)

# 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.8638 m, RMSE:33.4603 m
INFO:root:updated misfit RMSE: 4.1037
INFO:root:updated L2-norm: 2.0258, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.3203, tolerance: 1.005
INFO:root:
 ####################################
 iteration 2
INFO:root:Layer correction median: 12.4942 m, RMSE:18.7286 m
INFO:root:updated misfit RMSE: 2.4567
INFO:root:updated L2-norm: 1.5674, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2924, tolerance: 1.005
INFO:root:
 ####################################
 iteration 3
INFO:root:Layer correction median: 6.9882 m, RMSE:10.8116 m
INFO:root:updated misfit RMSE: 1.5601
INFO:root:updated L2-norm: 1.249, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2549, tolerance: 1.005
INFO:root:
 ####################################
 iteration 4
INFO:root:Layer correction median: 4.0601 m, RMSE:6.5069 m
INFO:root:updated misfit RMSE: 1.0644
INFO:root:updated L2-norm: 1.0317, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2107, tolerance: 1.005
INFO:root:
 ####################################
 iteration 5
INFO:root:Layer correction median: 2.2559 m, RMSE:4.1179 m
INFO:root:updated misfit RMSE: 0.7832
INFO:root:updated L2-norm: 0.885, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1658, tolerance: 1.005
INFO:root:
 ####################################
 iteration 6
INFO:root:Layer correction median: 1.2827 m, RMSE:2.7533 m
INFO:root:updated misfit RMSE: 0.6173
INFO:root:updated L2-norm: 0.7857, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1264, tolerance: 1.005
INFO:root:
 ####################################
 iteration 7
INFO:root:Layer correction median: 0.7539 m, RMSE:1.9435 m
INFO:root:updated misfit RMSE: 0.5144
INFO:root:updated L2-norm: 0.7172, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0955, tolerance: 1.005
INFO:root:
 ####################################
 iteration 8
INFO:root:Layer correction median: 0.4439 m, RMSE:1.4409 m
INFO:root:updated misfit RMSE: 0.447
INFO:root:updated L2-norm: 0.6686, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0728, tolerance: 1.005
INFO:root:
 ####################################
 iteration 9
INFO:root:Layer correction median: 0.2756 m, RMSE:1.1138 m
INFO:root:updated misfit RMSE: 0.4004
INFO:root:updated L2-norm: 0.6328, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0566, tolerance: 1.005
INFO:root:
 ####################################
 iteration 10
INFO:root:Layer correction median: 0.1832 m, RMSE:0.8914 m
INFO:root:updated misfit RMSE: 0.3667
INFO:root:updated L2-norm: 0.6055, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.045, tolerance: 1.005
INFO:root:
 ####################################
 iteration 11
INFO:root:Layer correction median: 0.1276 m, RMSE:0.7343 m
INFO:root:updated misfit RMSE: 0.3412
INFO:root:updated L2-norm: 0.5841, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0366, tolerance: 1.005
INFO:root:
 ####################################
 iteration 12
INFO:root:Layer correction median: 0.0874 m, RMSE:0.6198 m
INFO:root:updated misfit RMSE: 0.3213
INFO:root:updated L2-norm: 0.5668, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0305, tolerance: 1.005
INFO:root:
 ####################################
 iteration 13
INFO:root:Layer correction median: 0.0643 m, RMSE:0.5341 m
INFO:root:updated misfit RMSE: 0.3052
INFO:root:updated L2-norm: 0.5524, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.026, tolerance: 1.005
INFO:root:
 ####################################
 iteration 14
INFO:root:Layer correction median: 0.0585 m, RMSE:0.4686 m
INFO:root:updated misfit RMSE: 0.2918
INFO:root:updated L2-norm: 0.5402, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0227, tolerance: 1.005
INFO:root:
 ####################################
 iteration 15
INFO:root:Layer correction median: 0.054 m, RMSE:0.4174 m
INFO:root:updated misfit RMSE: 0.2804
INFO:root:updated L2-norm: 0.5296, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0201, tolerance: 1.005
INFO:root:
 ####################################
 iteration 16
INFO:root:Layer correction median: 0.0441 m, RMSE:0.3768 m
INFO:root:updated misfit RMSE: 0.2706
INFO:root:updated L2-norm: 0.5202, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0181, tolerance: 1.005
INFO:root:
 ####################################
 iteration 17
INFO:root:Layer correction median: 0.0346 m, RMSE:0.3439 m
INFO:root:updated misfit RMSE: 0.2619
INFO:root:updated L2-norm: 0.5117, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0165, tolerance: 1.005
INFO:root:
 ####################################
 iteration 18
INFO:root:Layer correction median: 0.0274 m, RMSE:0.3169 m
INFO:root:updated misfit RMSE: 0.2541
INFO:root:updated L2-norm: 0.504, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0152, tolerance: 1.005
INFO:root:
 ####################################
 iteration 19
INFO:root:Layer correction median: 0.0283 m, RMSE:0.2945 m
INFO:root:updated misfit RMSE: 0.247
INFO:root:updated L2-norm: 0.497, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0142, tolerance: 1.005
INFO:root:
Inversion terminated after 19 iterations because L2-norm (0.4969812889678823) was less then set tolerance: 0.5
Change parameter 'l2_norm_tolerance' if desired.

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,
)

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),
)

Run inversion with poor choice of density contrast

# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.INFO)

# use the best zref value
zref = 300

# create flat topography grid with a constant height
starting_topography = xr.full_like(true_topography, zref)

# re-calculate density grid with the best zref
density = xr.where(
    starting_topography >= zref,
    kwargs.get("density_contrast"),
    -kwargs.get("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=False,
)

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

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

# remove the regional from the misfit to get the residual
grav_df["res"] = grav_df.misfit - grav_df.reg

# update starting model in kwargs
kwargs["prism_layer"] = starting_prisms

# make new kwargs without zref or density
new_kwargs = {
    key: value
    for key, value in kwargs.items()
    if key
    not in [
        "zref",
        "density_contrast",
    ]
}
results = inversion.run_inversion(
    input_grav=grav_df,
    plot_convergence=True,
    zref=zref,
    density_contrast=kwargs.get("density_contrast"),
    **new_kwargs,
)

# 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.9009 m, RMSE:33.584 m
INFO:root:updated misfit RMSE: 4.2877
INFO:root:updated L2-norm: 2.0707, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2917, tolerance: 1.005
INFO:root:
 ####################################
 iteration 2
INFO:root:Layer correction median: 13.0884 m, RMSE:19.6149 m
INFO:root:updated misfit RMSE: 2.6947
INFO:root:updated L2-norm: 1.6415, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2614, tolerance: 1.005
INFO:root:
 ####################################
 iteration 3
INFO:root:Layer correction median: 7.6334 m, RMSE:11.834 m
INFO:root:updated misfit RMSE: 1.8038
INFO:root:updated L2-norm: 1.3431, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.2222, tolerance: 1.005
INFO:root:
 ####################################
 iteration 4
INFO:root:Layer correction median: 4.4461 m, RMSE:7.4625 m
INFO:root:updated misfit RMSE: 1.297
INFO:root:updated L2-norm: 1.1389, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1793, tolerance: 1.005
INFO:root:
 ####################################
 iteration 5
INFO:root:Layer correction median: 2.5975 m, RMSE:4.961 m
INFO:root:updated misfit RMSE: 0.9991
INFO:root:updated L2-norm: 0.9995, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1394, tolerance: 1.005
INFO:root:
 ####################################
 iteration 6
INFO:root:Layer correction median: 1.4566 m, RMSE:3.4853 m
INFO:root:updated misfit RMSE: 0.8153
INFO:root:updated L2-norm: 0.9029, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.107, tolerance: 1.005
INFO:root:
 ####################################
 iteration 7
INFO:root:Layer correction median: 0.8325 m, RMSE:2.5772 m
INFO:root:updated misfit RMSE: 0.6955
INFO:root:updated L2-norm: 0.834, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0827, tolerance: 1.005
INFO:root:
 ####################################
 iteration 8
INFO:root:Layer correction median: 0.5254 m, RMSE:1.9897 m
INFO:root:updated misfit RMSE: 0.6133
INFO:root:updated L2-norm: 0.7831, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.065, tolerance: 1.005
INFO:root:
 ####################################
 iteration 9
INFO:root:Layer correction median: 0.3399 m, RMSE:1.59 m
INFO:root:updated misfit RMSE: 0.5542
INFO:root:updated L2-norm: 0.7444, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.052, tolerance: 1.005
INFO:root:
 ####################################
 iteration 10
INFO:root:Layer correction median: 0.2372 m, RMSE:1.3057 m
INFO:root:updated misfit RMSE: 0.5101
INFO:root:updated L2-norm: 0.7142, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0423, tolerance: 1.005
INFO:root:
 ####################################
 iteration 11
INFO:root:Layer correction median: 0.1744 m, RMSE:1.0958 m
INFO:root:updated misfit RMSE: 0.4763
INFO:root:updated L2-norm: 0.6902, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0349, tolerance: 1.005
INFO:root:
 ####################################
 iteration 12
INFO:root:Layer correction median: 0.1277 m, RMSE:0.9362 m
INFO:root:updated misfit RMSE: 0.4497
INFO:root:updated L2-norm: 0.6706, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0292, tolerance: 1.005
INFO:root:
 ####################################
 iteration 13
INFO:root:Layer correction median: 0.0991 m, RMSE:0.812 m
INFO:root:updated misfit RMSE: 0.4283
INFO:root:updated L2-norm: 0.6544, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0247, tolerance: 1.005
INFO:root:
 ####################################
 iteration 14
INFO:root:Layer correction median: 0.0747 m, RMSE:0.7135 m
INFO:root:updated misfit RMSE: 0.4107
INFO:root:updated L2-norm: 0.6409, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0211, tolerance: 1.005
INFO:root:
 ####################################
 iteration 15
INFO:root:Layer correction median: 0.0702 m, RMSE:0.6341 m
INFO:root:updated misfit RMSE: 0.3961
INFO:root:updated L2-norm: 0.6294, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0183, tolerance: 1.005
INFO:root:
 ####################################
 iteration 16
INFO:root:Layer correction median: 0.0607 m, RMSE:0.5692 m
INFO:root:updated misfit RMSE: 0.3837
INFO:root:updated L2-norm: 0.6194, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.016, tolerance: 1.005
INFO:root:
 ####################################
 iteration 17
INFO:root:Layer correction median: 0.0461 m, RMSE:0.5155 m
INFO:root:updated misfit RMSE: 0.373
INFO:root:updated L2-norm: 0.6108, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0142, tolerance: 1.005
INFO:root:
 ####################################
 iteration 18
INFO:root:Layer correction median: 0.04 m, RMSE:0.4707 m
INFO:root:updated misfit RMSE: 0.3638
INFO:root:updated L2-norm: 0.6031, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0127, tolerance: 1.005
INFO:root:
 ####################################
 iteration 19
INFO:root:Layer correction median: 0.0346 m, RMSE:0.4329 m
INFO:root:updated misfit RMSE: 0.3556
INFO:root:updated L2-norm: 0.5963, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0114, tolerance: 1.005
INFO:root:
 ####################################
 iteration 20
INFO:root:Layer correction median: 0.0327 m, RMSE:0.4008 m
INFO:root:updated misfit RMSE: 0.3483
INFO:root:updated L2-norm: 0.5902, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0104, tolerance: 1.005
INFO:root:
 ####################################
 iteration 21
INFO:root:Layer correction median: 0.0342 m, RMSE:0.3733 m
INFO:root:updated misfit RMSE: 0.3418
INFO:root:updated L2-norm: 0.5846, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0095, tolerance: 1.005
INFO:root:
 ####################################
 iteration 22
INFO:root:Layer correction median: 0.0368 m, RMSE:0.3495 m
INFO:root:updated misfit RMSE: 0.3359
INFO:root:updated L2-norm: 0.5796, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0088, tolerance: 1.005
INFO:root:
 ####################################
 iteration 23
INFO:root:Layer correction median: 0.0322 m, RMSE:0.3289 m
INFO:root:updated misfit RMSE: 0.3305
INFO:root:updated L2-norm: 0.5749, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0082, tolerance: 1.005
INFO:root:
 ####################################
 iteration 24
INFO:root:Layer correction median: 0.03 m, RMSE:0.3108 m
INFO:root:updated misfit RMSE: 0.3255
INFO:root:updated L2-norm: 0.5705, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0076, tolerance: 1.005
INFO:root:
 ####################################
 iteration 25
INFO:root:Layer correction median: 0.0284 m, RMSE:0.2949 m
INFO:root:updated misfit RMSE: 0.3209
INFO:root:updated L2-norm: 0.5665, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0072, tolerance: 1.005
INFO:root:
 ####################################
 iteration 26
INFO:root:Layer correction median: 0.0278 m, RMSE:0.2807 m
INFO:root:updated misfit RMSE: 0.3166
INFO:root:updated L2-norm: 0.5627, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0068, tolerance: 1.005
INFO:root:
 ####################################
 iteration 27
INFO:root:Layer correction median: 0.0261 m, RMSE:0.2682 m
INFO:root:updated misfit RMSE: 0.3126
INFO:root:updated L2-norm: 0.5591, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0064, tolerance: 1.005
INFO:root:
 ####################################
 iteration 28
INFO:root:Layer correction median: 0.0233 m, RMSE:0.2569 m
INFO:root:updated misfit RMSE: 0.3088
INFO:root:updated L2-norm: 0.5557, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0061, tolerance: 1.005
INFO:root:
 ####################################
 iteration 29
INFO:root:Layer correction median: 0.0209 m, RMSE:0.2467 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.0058, tolerance: 1.005
INFO:root:
 ####################################
 iteration 30
INFO:root:Layer correction median: 0.0193 m, RMSE:0.2375 m
INFO:root:updated misfit RMSE: 0.3019
INFO:root:updated L2-norm: 0.5494, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0056, tolerance: 1.005
INFO:root:
Inversion terminated after 30 iterations with L2-norm=0.55 because maximum number of iterations (30) reached.

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,
)

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),
)