Ensemble 2: gravity spacings vs noise; no regional field¶

*with density estimation

This notebooks performs 100 synthetic inversions, with all combinations of 10 values of gravity flight line spacing, and 10 values of the gravity data noise. There is no regional gravity field. For each inversion, we estimate the optimal density contrast value for the seafloor with a cross-validation of the constraint points.

import packages

[1]:

%load_ext autoreload
%autoreload 2
import itertools
import logging
import os
import pathlib
import pickle

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pygmt
import scipy
import verde as vd
import xarray as xr
from invert4geom import inversion, optimization, plotting, utils
from invert4geom import synthetic as inv_synthetic
from polartoolkit import maps
from polartoolkit import utils as polar_utils
from tqdm.autonotebook import tqdm

import synthetic_bathymetry_inversion.plotting as synth_plotting
from synthetic_bathymetry_inversion import synthetic

os.environ["POLARTOOLKIT_HEMISPHERE"] = "south"

logging.getLogger().setLevel(logging.INFO)

/home/sungw937/synthetic_bathymetry_inversion/.pixi/envs/default/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.

[2]:

ensemble_path = "../results/Ross_Sea/ensembles/Ross_Sea_ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation"
ensemble_fname = f"{ensemble_path}.csv"

[3]:

# set grid parameters
spacing = 2e3
inversion_region = (-40e3, 110e3, -1600e3, -1400e3)

true_density_contrast = 1476

bathymetry, _, original_grav_df = synthetic.load_synthetic_model(
    spacing=spacing,
    inversion_region=inversion_region,
    buffer=spacing * 10,
    basement=False,
    zref=0,
    bathymetry_density_contrast=true_density_contrast,
)
buffer_region = polar_utils.get_grid_info(bathymetry)[1]

requested spacing (2000.0) is smaller than the original (5000.0).

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_4_1.png

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_4_3.png

[4]:

bathymetry = bathymetry.sel(
    easting=slice(inversion_region[0], inversion_region[1]),
    northing=slice(inversion_region[2], inversion_region[3]),
)

[5]:

original_grav_df = original_grav_df.drop(columns=["basement_grav", "disturbance"])
original_grav_df = original_grav_df.rename(
    columns={"gravity_anomaly": "gravity_anomaly_full_res_no_noise"},
)
original_grav_df

[5]:

	northing	easting	upward	bathymetry_grav	gravity_anomaly_full_res_no_noise
0	-1600000.0	-40000.0	1000.0	-35.551085	-35.551085
1	-1600000.0	-38000.0	1000.0	-36.054683	-36.054683
2	-1600000.0	-36000.0	1000.0	-36.473168	-36.473168
3	-1600000.0	-34000.0	1000.0	-36.755627	-36.755627
4	-1600000.0	-32000.0	1000.0	-36.951045	-36.951045
...	...	...	...	...	...
7671	-1400000.0	102000.0	1000.0	-25.760090	-25.760090
7672	-1400000.0	104000.0	1000.0	-25.911429	-25.911429
7673	-1400000.0	106000.0	1000.0	-26.032814	-26.032814
7674	-1400000.0	108000.0	1000.0	-26.121903	-26.121903
7675	-1400000.0	110000.0	1000.0	-26.206160	-26.206160

7676 rows × 5 columns

[6]:

# semi-regularly spaced
constraint_points = synthetic.constraint_layout_number(
    shape=(6, 7),
    region=inversion_region,
    padding=-spacing,
    shapefile="../results/Ross_Sea/Ross_Sea_outline.shp",
    add_outside_points=True,
    grid_spacing=spacing,
)

# sample true topography at these points
constraint_points = utils.sample_grids(
    constraint_points,
    bathymetry,
    "true_upward",
    coord_names=("easting", "northing"),
)
constraint_points["upward"] = constraint_points.true_upward
constraint_points.head()

[6]:

	northing	easting	inside	true_upward	upward
0	-1600000.0	-40000.0	False	-601.093994	-601.093994
1	-1600000.0	-38000.0	False	-609.216919	-609.216919
2	-1600000.0	-36000.0	False	-616.355957	-616.355957
3	-1600000.0	-34000.0	False	-621.262268	-621.262268
4	-1600000.0	-32000.0	False	-625.510925	-625.510925

[7]:

# grid the sampled values using verde
starting_topography_kwargs = dict(
    method="splines",
    region=buffer_region,
    spacing=spacing,
    constraints_df=constraint_points,
    dampings=None,
)

starting_bathymetry = utils.create_topography(**starting_topography_kwargs)

starting_bathymetry

[7]:

<xarray.DataArray 'scalars' (northing: 121, easting: 96)> Size: 93kB
array([[-541.24413869, -544.57181187, -547.92293689, ..., -360.00006254,
        -357.06767408, -354.19957767],
       [-543.34402687, -546.81675803, -550.35256332, ..., -362.90253226,
        -359.96874159, -357.11431886],
       [-545.05533622, -548.66036837, -552.37518162, ..., -365.66137905,
        -362.73269531, -359.90052825],
       ...,
       [-591.95335283, -595.51882199, -599.06869705, ..., -440.89315875,
        -440.6944619 , -440.40553782],
       [-590.53134833, -594.09076637, -597.64079287, ..., -440.69158328,
        -440.42525249, -440.07197234],
       [-589.1663267 , -592.73504777, -596.30209679, ..., -440.51760947,
        -440.1713932 , -439.74434038]], shape=(121, 96))
Coordinates:
  * northing  (northing) float64 968B -1.62e+06 -1.618e+06 ... -1.38e+06
  * easting   (easting) float64 768B -6e+04 -5.8e+04 ... 1.28e+05 1.3e+05
Attributes:
    metadata:  Generated by SplineCV(cv=KFold(n_splits=5, random_state=0, shu...
    damping:   None

xarray.DataArray

'scalars'

northing: 121
easting: 96

-541.2 -544.6 -547.9 -551.3 -554.5 ... -440.8 -440.5 -440.2 -439.7

array([[-541.24413869, -544.57181187, -547.92293689, ..., -360.00006254,
        -357.06767408, -354.19957767],
       [-543.34402687, -546.81675803, -550.35256332, ..., -362.90253226,
        -359.96874159, -357.11431886],
       [-545.05533622, -548.66036837, -552.37518162, ..., -365.66137905,
        -362.73269531, -359.90052825],
       ...,
       [-591.95335283, -595.51882199, -599.06869705, ..., -440.89315875,
        -440.6944619 , -440.40553782],
       [-590.53134833, -594.09076637, -597.64079287, ..., -440.69158328,
        -440.42525249, -440.07197234],
       [-589.1663267 , -592.73504777, -596.30209679, ..., -440.51760947,
        -440.1713932 , -439.74434038]], shape=(121, 96))

Coordinates: (2)

northing

(northing)

float64

-1.62e+06 -1.618e+06 ... -1.38e+06

array([-1620000., -1618000., -1616000., -1614000., -1612000., -1610000.,
       -1608000., -1606000., -1604000., -1602000., -1600000., -1598000.,
       -1596000., -1594000., -1592000., -1590000., -1588000., -1586000.,
       -1584000., -1582000., -1580000., -1578000., -1576000., -1574000.,
       -1572000., -1570000., -1568000., -1566000., -1564000., -1562000.,
       -1560000., -1558000., -1556000., -1554000., -1552000., -1550000.,
       -1548000., -1546000., -1544000., -1542000., -1540000., -1538000.,
       -1536000., -1534000., -1532000., -1530000., -1528000., -1526000.,
       -1524000., -1522000., -1520000., -1518000., -1516000., -1514000.,
       -1512000., -1510000., -1508000., -1506000., -1504000., -1502000.,
       -1500000., -1498000., -1496000., -1494000., -1492000., -1490000.,
       -1488000., -1486000., -1484000., -1482000., -1480000., -1478000.,
       -1476000., -1474000., -1472000., -1470000., -1468000., -1466000.,
       -1464000., -1462000., -1460000., -1458000., -1456000., -1454000.,
       -1452000., -1450000., -1448000., -1446000., -1444000., -1442000.,
       -1440000., -1438000., -1436000., -1434000., -1432000., -1430000.,
       -1428000., -1426000., -1424000., -1422000., -1420000., -1418000.,
       -1416000., -1414000., -1412000., -1410000., -1408000., -1406000.,
       -1404000., -1402000., -1400000., -1398000., -1396000., -1394000.,
       -1392000., -1390000., -1388000., -1386000., -1384000., -1382000.,
       -1380000.])

easting

(easting)

float64

-6e+04 -5.8e+04 ... 1.3e+05

array([-60000., -58000., -56000., -54000., -52000., -50000., -48000., -46000.,
       -44000., -42000., -40000., -38000., -36000., -34000., -32000., -30000.,
       -28000., -26000., -24000., -22000., -20000., -18000., -16000., -14000.,
       -12000., -10000.,  -8000.,  -6000.,  -4000.,  -2000.,      0.,   2000.,
         4000.,   6000.,   8000.,  10000.,  12000.,  14000.,  16000.,  18000.,
        20000.,  22000.,  24000.,  26000.,  28000.,  30000.,  32000.,  34000.,
        36000.,  38000.,  40000.,  42000.,  44000.,  46000.,  48000.,  50000.,
        52000.,  54000.,  56000.,  58000.,  60000.,  62000.,  64000.,  66000.,
        68000.,  70000.,  72000.,  74000.,  76000.,  78000.,  80000.,  82000.,
        84000.,  86000.,  88000.,  90000.,  92000.,  94000.,  96000.,  98000.,
       100000., 102000., 104000., 106000., 108000., 110000., 112000., 114000.,
       116000., 118000., 120000., 122000., 124000., 126000., 128000., 130000.])

Attributes: (2)
metadata :
Generated by SplineCV(cv=KFold(n_splits=5, random_state=0, shuffle=True), dampings=[None], mindists=[0])
damping :
None

[10]:

inner_starting_bathymetry = starting_bathymetry.sel(
    easting=slice(inversion_region[0], inversion_region[1]),
    northing=slice(inversion_region[2], inversion_region[3]),
)

starting_bathymetry_rmse = utils.rmse(bathymetry - inner_starting_bathymetry)
starting_bathymetry_rmse

[10]:

np.float64(25.97734942135522)

[9]:

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

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

RMSE: 0.03 m

[10]:

# compare starting and actual bathymetry grids
grids = polar_utils.grd_compare(
    bathymetry,
    starting_bathymetry,
    fig_height=10,
    plot=True,
    cmap="rain",
    reverse_cpt=True,
    diff_cmap="balance+h0",
    grid1_name="True bathymetry",
    grid2_name="Starting bathymetry",
    title="Difference",
    title_font="18p,Helvetica-Bold,black",
    cbar_unit="m",
    cbar_label="elevation",
    RMSE_decimals=0,
    region=inversion_region,
    inset=False,
    hist=True,
    label_font="16p,Helvetica,black",
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
    points_style="x.2c",
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_11_0.png

[11]:

# delete files if already exist
for p in pathlib.Path().glob(f"{ensemble_path}_starting_bathymetry.nc"):
    p.unlink(missing_ok=True)

# save to file
starting_bathymetry.drop_attrs().to_netcdf(
    f"{ensemble_path}_starting_bathymetry.nc",
)
constraint_points.to_csv(f"{ensemble_path}_constraint_points.csv", index=False)

Set line spacings and noise levels¶

[11]:

num = 10

# Define number of flights lines on log scale
grav_line_numbers = np.unique(np.round(np.geomspace(3, 30, num)))
grav_line_numbers = [int(i) for i in grav_line_numbers]
assert len(grav_line_numbers) == num

# Define noise levels for grav data
grav_noise_levels = [float(round(x, 2)) for x in np.linspace(0, 5, num)]

print("number of grav lines:", grav_line_numbers)
print("grav noise levels:", grav_noise_levels)

number of grav lines: [3, 4, 5, 6, 8, 11, 14, 18, 23, 30]
grav noise levels: [0.0, 0.56, 1.11, 1.67, 2.22, 2.78, 3.33, 3.89, 4.44, 5.0]

[12]:

# turn into dataframe
sampled_params_df = pd.DataFrame(
    itertools.product(
        grav_line_numbers,
        grav_noise_levels,
    ),
    columns=[
        "grav_line_numbers",
        "grav_noise_levels",
    ],
)

sampled_params_dict = dict(
    grav_line_numbers=dict(sampled_values=sampled_params_df.grav_line_numbers),
    grav_noise_levels=dict(sampled_values=sampled_params_df.grav_noise_levels),
)

plotting.plot_latin_hypercube(
    sampled_params_dict,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_15_0.png

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_15_1.png

[13]:

sampled_params_df

[13]:

	grav_line_numbers	grav_noise_levels
0	3	0.00
1	3	0.56
2	3	1.11
3	3	1.67
4	3	2.22
...	...	...
95	30	2.78
96	30	3.33
97	30	3.89
98	30	4.44
99	30	5.00

100 rows × 2 columns

Contaminate gravity with noise¶

[14]:

# fnames to save files for each ensemble
sampled_params_df["grav_df_fname"] = pd.Series()

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # set file names, add to dataframe
    sampled_params_df.loc[i, "grav_df_fname"] = f"{ensemble_path}_grav_df_{i}.csv"

[15]:

grav_grid = original_grav_df.set_index(["northing", "easting"]).to_xarray()

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    grav_df = original_grav_df.copy()

    # contaminated with long-wavelength noise
    grav_df["gravity_anomaly_full_res"] = grav_df.gravity_anomaly_full_res_no_noise
    with utils._log_level(logging.ERROR):
        contaminated = inv_synthetic.contaminate_with_long_wavelength_noise(
            grav_grid.gravity_anomaly_full_res_no_noise,
            coarsen_factor=None,
            spacing=2e3,
            noise_as_percent=False,
            noise=row.grav_noise_levels,
        )
    contaminated_df = vd.grid_to_table(
        contaminated.rename("gravity_anomaly_full_res")
    ).reset_index(drop=True)

    grav_df = pd.merge(  # noqa: PD015
        grav_df.drop(columns=["gravity_anomaly_full_res"], errors="ignore"),
        contaminated_df,
        on=["easting", "northing"],
    )

    # # short-wavelength noise
    # with utils._log_level(logging.ERROR):
    #     contaminated = synthetic.contaminate_with_long_wavelength_noise(
    #         grav_df.set_index(["northing", "easting"])
    #         .to_xarray()
    #         .gravity_anomaly_full_res,
    #         coarsen_factor=None,
    #         spacing=spacing * 2,
    #         noise_as_percent=False,
    #         noise=row.grav_noise_levels,
    #     )
    # contaminated_df = vd.grid_to_table(
    #     contaminated.rename("gravity_anomaly_full_res")
    # ).reset_index(drop=True)

    # grav_df = pd.merge(
    #     grav_df.drop(columns=["gravity_anomaly_full_res"], errors="ignore"),
    #     contaminated_df,
    #     on=["easting", "northing"],
    # )

    # add noise level to uncert column
    grav_df["uncert"] = row.grav_noise_levels

    # save to files
    grav_df.to_csv(row.grav_df_fname, index=False)

sampled_params_df.head()

[15]:

	grav_line_numbers	grav_noise_levels	grav_df_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[16]:

sampled_params_df.to_csv(ensemble_fname, index=False)

[17]:

sampled_params_df = pd.read_csv(ensemble_fname)
sampled_params_df.head()

[17]:

	grav_line_numbers	grav_noise_levels	grav_df_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[18]:

# subplots showing grav data loss from sampling
grids = []
titles = []
cbar_labels = []
for i, row in sampled_params_df.iterrows():
    if i < num:
        grav_df = pd.read_csv(row.grav_df_fname)
        grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()
        dif = (
            grav_grid.gravity_anomaly_full_res_no_noise
            - grav_grid.gravity_anomaly_full_res
        )
        # add to lists
        grids.append(dif)
        titles.append(f"noise level: {round(row.grav_noise_levels, 1)} mGal")
        cbar_labels.append(f"RMSE: {round(utils.rmse(dif), 2)} (mGal)")

fig = maps.subplots(
    grids,
    fig_title="Gravity noise",
    titles=titles,
    cbar_labels=cbar_labels,
    cmap="balance+h0",
    hist=True,
    cpt_lims=polar_utils.get_combined_min_max(grids, robust=True),
    cbar_font="18p,Helvetica,black",
    yshift_amount=-1.2,
)

fig.show()

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_22_0.png

Sample along flightlines¶

[19]:

# fnames to save files for each ensemble
sampled_params_df["grav_survey_df_fname"] = pd.Series()

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # set file names, add to dataframe
    sampled_params_df.loc[i, "grav_survey_df_fname"] = (
        f"{ensemble_path}_grav_survey_df_{i}.csv"
    )

[20]:

logging.getLogger().setLevel(logging.WARNING)

# add empty columns
# average of the east/west and north/south line spacings
sampled_params_df["grav_line_spacing"] = pd.Series()

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_df = pd.read_csv(row.grav_df_fname)

    # create flight lines
    grav_survey_df = synthetic.airborne_survey(
        along_line_spacing=500,
        grav_observation_height=1e3,
        ns_line_number=row.grav_line_numbers,
        ew_line_number=row.grav_line_numbers,
        region=inversion_region,
        grav_grid=grav_df.set_index(["northing", "easting"])
        .to_xarray()
        .gravity_anomaly_full_res,
        plot=False,
    )
    x_spacing = (inversion_region[1] - inversion_region[0]) / row.grav_line_numbers
    y_spacing = (inversion_region[3] - inversion_region[2]) / row.grav_line_numbers
    grav_line_spacing = ((x_spacing + y_spacing) / 2) / 1e3

    # sample no-noise grid onto survey lines
    grav_survey_df = utils.sample_grids(
        grav_survey_df,
        grav_df.set_index(
            ["northing", "easting"]
        ).gravity_anomaly_full_res_no_noise.to_xarray(),
        sampled_name="gravity_anomaly_no_noise",
    )

    # add to dataframe
    sampled_params_df.loc[i, "grav_line_spacing"] = grav_line_spacing

    # save to files
    grav_survey_df.to_csv(row.grav_survey_df_fname, index=False)


sampled_params_df

[20]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333
...	...	...	...	...	...
95	30	2.78	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333
96	30	3.33	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333
97	30	3.89	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333
98	30	4.44	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333
99	30	5.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333

100 rows × 5 columns

[21]:

for i, row in sampled_params_df.iterrows():
    # load data
    grav_survey_df = pd.read_csv(row.grav_survey_df_fname)

    # calc min distance
    coords = vd.grid_coordinates(
        region=inversion_region,
        spacing=100,
    )
    grid = vd.make_xarray_grid(coords, np.ones_like(coords[0]), data_names="z").z
    min_dist = utils.dist_nearest_points(
        grav_survey_df,
        grid,
    ).min_dist
    grav_proximity = min_dist.median().to_numpy() / 1e3
    sampled_params_df.loc[i, "grav_proximity"] = grav_proximity
    sampled_params_df.loc[i, "grav_proximity_skew"] = scipy.stats.skew(
        min_dist.to_numpy().ravel(), nan_policy="omit"
    )

(
    sampled_params_df.grav_proximity.describe(),
    sampled_params_df.grav_proximity_skew.describe(),
)

[21]:

(count    100.000000
 mean       3.430303
 std        2.368128
 min        0.824621
 25%        1.400000
 50%        2.704055
 75%        5.001000
 max        8.333333
 Name: grav_proximity, dtype: float64,
 count    100.000000
 mean       0.504258
 std        0.015236
 min        0.489544
 25%        0.491777
 50%        0.498486
 75%        0.513887
 max        0.537052
 Name: grav_proximity_skew, dtype: float64)

[22]:

w = np.abs(inversion_region[1] - inversion_region[0]) / 1e3
h = np.abs(inversion_region[3] - inversion_region[2]) / 1e3
w, h

[22]:

(np.float64(150.0), np.float64(200.0))

[23]:

sampled_params_df["flight_kms"] = pd.Series()
for i, row in sampled_params_df.iterrows():
    survey_df = pd.read_csv(row.grav_survey_df_fname)
    # flight km is length of lines in km X number of lines
    sampled_params_df.loc[i, "flight_kms"] = (
        w * row.grav_line_numbers + h * row.grav_line_numbers
    )
sampled_params_df.flight_kms.unique()

[23]:

array([np.float64(1050.0), np.float64(1400.0), np.float64(1750.0),
       np.float64(2100.0), np.float64(2800.0), np.float64(3850.0),
       np.float64(4900.0), np.float64(6300.0), np.float64(8050.0),
       np.float64(10500.0)], dtype=object)

[24]:

inversion_area = (
    (inversion_region[1] - inversion_region[0])
    / 1e3
    * (inversion_region[3] - inversion_region[2])
    / 1e3
)
sampled_params_df["flight_kms_per_10000sq_km"] = (
    (sampled_params_df.flight_kms / inversion_area) * 10e3
).astype(float)
sampled_params_df.flight_kms_per_10000sq_km.unique()

[24]:

array([ 350.        ,  466.66666667,  583.33333333,  700.        ,
        933.33333333, 1283.33333333, 1633.33333333, 2100.        ,
       2683.33333333, 3500.        ])

[25]:

df = sampled_params_df
y = "grav_proximity"
x = "grav_proximity_skew"

plt.plot(df[x], df[y], marker="o")
ax = plt.gca()
ax.set_xlabel(x)
ax.set_ylabel(y)

[25]:

Text(0, 0.5, 'grav_proximity')

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_30_1.png

[26]:

df = sampled_params_df
y = "grav_proximity"
x = "flight_kms"

plt.plot(df[x], df[y], marker="o")
ax = plt.gca()
ax.set_xlabel(x)
ax.set_ylabel(y)

[26]:

Text(0, 0.5, 'grav_proximity')

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_31_1.png

[27]:

df = sampled_params_df
y = "grav_proximity_skew"
x = "flight_kms"

plt.plot(df[x], df[y], marker="o")
ax = plt.gca()
ax.set_xlabel(x)
ax.set_ylabel(y)

[27]:

Text(0, 0.5, 'grav_proximity_skew')

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_32_1.png

Filter flight line data¶

[28]:

# add empty columns
sampled_params_df["filter_width_trials_fname"] = np.nan

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # set file names, add to dataframe
    sampled_params_df.loc[i, "filter_width_trials_fname"] = (
        f"{ensemble_path}_filter_width_trials_{i}.csv"
    )

[29]:

logging.getLogger().setLevel(logging.WARNING)

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_df = pd.read_csv(row.grav_df_fname)
    grav_survey_df = pd.read_csv(row.grav_survey_df_fname)

    # trial a range of 1D low-pass filters and see which performs best
    dfs = []
    filter_widths = np.arange(0, 80e3, 4e3)
    for f in filter_widths:
        survey_df = grav_survey_df.copy()

        if f == 0:
            survey_df["grav_anomaly_filt"] = survey_df.gravity_anomaly
        else:
            survey_df["grav_anomaly_filt"] = synthetic.filter_flight_lines(
                survey_df,
                data_column="gravity_anomaly",
                filt_type=f"g{f}",
            )

        # compared filtered line data with no-noise line data
        survey_df["dif"] = (
            survey_df.gravity_anomaly_no_noise - survey_df.grav_anomaly_filt
        )

        dfs.append(
            pd.DataFrame(
                [
                    {
                        "filt_width": f,
                        "rmse": utils.rmse(survey_df.dif),
                    }
                ]
            )
        )

    filter_width_trials = pd.concat(dfs, ignore_index=True)

    # save to files
    filter_width_trials.to_csv(row.filter_width_trials_fname, index=False)


sampled_params_df

[29]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
...	...	...	...	...	...	...	...	...	...	...
95	30	2.78	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
96	30	3.33	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
97	30	3.89	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
98	30	4.44	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
99	30	5.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

100 rows × 10 columns

[30]:

# add empty columns
sampled_params_df["best_filter_width"] = np.nan

norm = plt.Normalize(
    vmin=sampled_params_df.grav_noise_levels.to_numpy().min(),
    vmax=sampled_params_df.grav_noise_levels.to_numpy().max(),
)
for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    filter_width_trials = pd.read_csv(row.filter_width_trials_fname)

    best_ind = filter_width_trials.rmse.idxmin()
    best_filter_width = filter_width_trials.iloc[best_ind].filt_width

    if i == 0:
        fig, ax = plt.subplots()

    ax.plot(
        filter_width_trials.filt_width,
        filter_width_trials.rmse,
        color=plt.cm.viridis(norm(row.grav_noise_levels)),  # pylint: disable=no-member
    )
    ax.set_xlabel("Filter width (m)")
    ax.set_ylabel("RMSE (mGal)")
    ax.scatter(
        x=filter_width_trials.filt_width.iloc[best_ind],
        y=filter_width_trials.rmse.iloc[best_ind],
        marker="*",
        edgecolor="black",
        linewidth=0.5,
        color=plt.cm.viridis(norm(row.grav_noise_levels)),  # pylint: disable=no-member
        s=60,
        zorder=20,
    )

    sampled_params_df.loc[i, "best_filter_width"] = best_filter_width

sm = plt.cm.ScalarMappable(cmap="viridis", norm=norm)
cax = fig.add_axes([0.93, 0.1, 0.05, 0.8])
cbar = plt.colorbar(sm, cax=cax)
cbar.set_label("noise level")

sampled_params_df

[30]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0
...	...	...	...	...	...	...	...	...	...	...	...
95	30	2.78	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0
96	30	3.33	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	40000.0
97	30	3.89	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	44000.0
98	30	4.44	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48000.0
99	30	5.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	52000.0

100 rows × 11 columns

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_36_2.png

[31]:

# variable = grav_noise_levels
# title = "Gravity noise (mGal)"

variable = "grav_line_spacing"
title = "Line spacing (m)"

norm = plt.Normalize(
    vmin=sampled_params_df[variable].to_numpy().min(),
    vmax=sampled_params_df[variable].to_numpy().max(),
)
for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    filter_width_trials = pd.read_csv(row.filter_width_trials_fname)

    best_ind = filter_width_trials.rmse.idxmin()
    best_filter_width = filter_width_trials.iloc[best_ind].filt_width

    if i == 0:
        fig, ax = plt.subplots()

    ax.plot(
        filter_width_trials.filt_width,
        filter_width_trials.rmse,
        color=plt.cm.viridis(norm(row[variable])),  # pylint: disable=no-member
    )
    ax.set_xlabel("Filter width (m)")
    ax.set_ylabel("RMSE (mGal)")
    ax.scatter(
        x=filter_width_trials.filt_width.iloc[best_ind],
        y=filter_width_trials.rmse.iloc[best_ind],
        marker="*",
        edgecolor="black",
        linewidth=0.5,
        color=plt.cm.viridis(norm(row[variable])),  # pylint: disable=no-member
        s=60,
        zorder=20,
    )

sm = plt.cm.ScalarMappable(cmap="viridis", norm=norm)
cax = fig.add_axes([0.93, 0.1, 0.05, 0.8])
cbar = plt.colorbar(sm, cax=cax)
cbar.set_label(title)

sampled_params_df

[31]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0
...	...	...	...	...	...	...	...	...	...	...	...
95	30	2.78	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0
96	30	3.33	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	40000.0
97	30	3.89	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	44000.0
98	30	4.44	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48000.0
99	30	5.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	52000.0

100 rows × 11 columns

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_37_2.png

[32]:

fig = synth_plotting.plot_2var_ensemble(
    sampled_params_df[sampled_params_df.grav_noise_levels > 0],
    x="grav_line_spacing",
    y="grav_noise_levels",
    x_title="Line spacing (km)",
    y_title="Gravity noise (mGal)",
    background="best_filter_width",
    background_title="Optimal filter width (m)",
    background_robust=True,
    constrained_layout=False,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_38_0.png

[33]:

sampled_params_df.best_filter_width.plot.hist()

[33]:

<Axes: ylabel='Frequency'>

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_39_1.png

[34]:

logging.getLogger().setLevel(logging.WARNING)

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_survey_df = pd.read_csv(row.grav_survey_df_fname)

    # filter each line in 1D with a Gaussian filter to remove some noise
    if row.best_filter_width == 0:
        grav_survey_df["gravity_anomaly"] = grav_survey_df.gravity_anomaly
    else:
        grav_survey_df["gravity_anomaly"] = synthetic.filter_flight_lines(
            grav_survey_df,
            data_column="gravity_anomaly",
            filt_type=f"g{row.best_filter_width}",
        )

    # save to files
    grav_survey_df.to_csv(row.grav_survey_df_fname, index=False)

Interpolate flight line data¶

[35]:

logging.getLogger().setLevel(logging.WARNING)

# add empty columns
sampled_params_df["grav_data_loss_rmse"] = np.nan
sampled_params_df["grav_data_loss_mae"] = np.nan

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_df = pd.read_csv(row.grav_df_fname)
    grav_survey_df = pd.read_csv(row.grav_survey_df_fname)

    grav_grid = pygmt.surface(
        data=grav_survey_df[["easting", "northing", "gravity_anomaly"]],
        region=inversion_region,
        spacing=spacing,
        tension=0.25,
        verbose="q",
    )
    grav_df = utils.sample_grids(
        grav_df,
        grav_grid,
        sampled_name="gravity_anomaly",
    )

    # add to dataframe
    dif = grav_df.gravity_anomaly_full_res_no_noise - grav_df.gravity_anomaly
    sampled_params_df.loc[i, "grav_data_loss_rmse"] = float(utils.rmse(dif))
    sampled_params_df.loc[i, "grav_data_loss_mae"] = float(np.abs(dif).mean())

    # save to files
    grav_survey_df.to_csv(row.grav_survey_df_fname, index=False)
    grav_df.to_csv(row.grav_df_fname, index=False)

sampled_params_df

[35]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866
...	...	...	...	...	...	...	...	...	...	...	...	...	...
95	30	2.78	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	0.777220	0.597805
96	30	3.33	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	40000.0	0.884233	0.678400
97	30	3.89	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	44000.0	0.988483	0.756277
98	30	4.44	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48000.0	1.087007	0.829231
99	30	5.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	5.833333	0.824621	0.537052	10500.0	3500.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	52000.0	1.184134	0.901034

100 rows × 13 columns

[36]:

sampled_params_df.to_csv(ensemble_fname, index=False)

[37]:

sampled_params_df = pd.read_csv(ensemble_fname)
sampled_params_df.head()

[37]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866

[38]:

fig = synth_plotting.plot_2var_ensemble(
    sampled_params_df,
    x="grav_line_spacing",
    y="grav_noise_levels",
    x_title="Line spacing (km)",
    y_title="Gravity noise (mGal)",
    background="grav_data_loss_rmse",
    background_title="Gravity data RMSE (m)",
    plot_title="Gravity spacing vs noise; no regional field",
    constrained_layout=False,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_45_0.png

[39]:

x = "grav_noise_levels"
fig = synth_plotting.plot_ensemble_as_lines(
    sampled_params_df,
    figsize=(4, 2),
    x=x,
    x_label="Gravity noise (mGal)",
    y="grav_data_loss_rmse",
    y_label="Gravity data RMSE (m)",
    groupby_col="grav_line_spacing",
    cbar_label="Line spacing (km)",
    # logx=True,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_46_0.png

[40]:

x = "grav_line_spacing"
fig = synth_plotting.plot_ensemble_as_lines(
    sampled_params_df,
    figsize=(4, 2),
    x=x,
    x_label="Line spacing (km)",
    y="grav_data_loss_rmse",
    y_label="Gravity data RMSE (m)",
    groupby_col="grav_noise_levels",
    cbar_label="Gravity noise (mGal)",
    # logx=True,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_47_0.png

[41]:

x = "grav_line_spacing"
fig = synth_plotting.plot_ensemble_as_lines(
    sampled_params_df,
    figsize=(4, 2),
    x=x,
    x_label="Line spacing (km)",
    y="grav_data_loss_rmse",
    y_label="Gravity data RMSE (m)",
    groupby_col="grav_noise_levels",
    cbar_label="Gravity noise (mGal)",
    logx=True,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_48_0.png

[42]:

# subplots showing grav data loss from sampling and noise
grids = []
cbar_labels = []
row_titles = []
column_titles = []
surveys = []

# iterate over the ensemble starting with high noise, low line numbers
# row per noise level and column per line spacing
for i, row in sampled_params_df.sort_values(
    by=["grav_noise_levels", "grav_line_spacing"],
    ascending=False,
).iterrows():
    grav_df = pd.read_csv(row.grav_df_fname)
    grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()
    dif = grav_grid.gravity_anomaly_full_res_no_noise - grav_grid.gravity_anomaly
    # add to lists
    grids.append(dif)
    cbar_labels.append(f"RMSE:{round(utils.rmse(dif), 2)}")
    if i % num == 0:
        column_titles.append(round(row.grav_line_spacing))
    if i < num:
        row_titles.append(round(row.grav_noise_levels, 1))
    # get survey points
    grav_survey_df = pd.read_csv(row.grav_survey_df_fname)
    surveys.append(grav_survey_df)

fig = maps.subplots(
    grids,
    fig_height=8,
    fig_title="Data loss from airborne surveying and noise",
    fig_title_font="100p,Helvetica-Bold,black",
    fig_x_axis_title="Line spacing (km)",
    fig_y_axis_title="Gravity noise (mGal)",
    fig_axis_title_font="80p,Helvetica-Bold,black",
    fig_title_y_offset="8c",
    # fig_x_axis_title_y_offset="-2.5c",
    fig_y_axis_title_x_offset="3c",
    cbar_labels=cbar_labels,
    cmap="balance+h0",
    # hist=True,
    cpt_lims=polar_utils.get_combined_min_max(grids, robust=True),
    cbar_font="18p,Helvetica,black",
    row_titles=row_titles,
    column_titles=column_titles,
    point_sets=surveys,
    points_style="p.02c",
    yshift_amount=-1.2,
)
fig.show()

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_49_0.png

[43]:

# subplots showing true gravity and sampled-noisy gravity
grids = []
cbar_labels = []
row_titles = []
column_titles = []

df = sampled_params_df
# iterate over the ensemble starting with high noise, small line spacing
# row per noise level and column per line spacing
for i, row in df.sort_values(
    by=["grav_noise_levels", "grav_line_spacing"], ascending=[False, True]
).iterrows():
    grav_df = pd.read_csv(row.grav_df_fname)
    grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

    # dif = grav_grid.gravity_anomaly_full_res_no_noise - grav_grid.gravity_anomaly

    dif = grav_grid.gravity_anomaly_full_res_no_noise - grav_grid.gravity_anomaly

    # add to lists
    grids.append(grav_grid.gravity_anomaly_full_res_no_noise)
    grids.append(grav_grid.gravity_anomaly)
    cbar_labels.append(f"RMSE:{round(utils.rmse(dif), 2)}")
    cbar_labels.append(" ")
    if i % np.sqrt(len(df)) == 0:
        column_titles.append(round(row.grav_line_spacing))
        column_titles.append(" ")
    if i < np.sqrt(len(df)):
        row_titles.append(round(row.grav_noise_levels, 1))
        # row_titles.append(" ")

fig = maps.subplots(
    grids,
    dims=(np.sqrt(len(df)), 2 * np.sqrt(len(df))),
    fig_height=8,
    fig_title="No regional",
    fig_x_axis_title="Flight line spacing (km)",
    fig_y_axis_title="Gravity noise (mGal)",
    fig_title_font="100p,Helvetica-Bold,black",
    fig_axis_title_font="80p,Helvetica-Bold,black",
    fig_title_y_offset="8c",
    fig_x_axis_title_y_offset="2.5c",
    fig_y_axis_title_x_offset="3c",
    cbar_labels=cbar_labels,
    cmaps=[x for xs in [["viridis"] * 2 for x in range(len(df))] for x in xs],
    reverse_cpt=True,
    # colorbar=False,
    # hist=True,
    # cpt_lims=polar_utils.get_combined_min_max(grids, robust=True),
    cpt_limits=[
        x
        for xs in [
            [polar_utils.get_combined_min_max(grids, robust=True)] * 2
            for x in range(len(df))
        ]
        for x in xs
    ],
    cbar_font="45p,Helvetica,black",
    row_titles=row_titles,
    column_titles=column_titles,
    row_titles_font="70p,Helvetica,black",
    column_titles_font="70p,Helvetica,black",
    yshift_amount=-1.3,
)

fig.show()

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_50_0.png

Density optimization¶

[44]:

for i, row in sampled_params_df.iterrows():
    # add results file name
    sampled_params_df.loc[i, "results_fname"] = f"{ensemble_path}_results_{i}"
    sampled_params_df.loc[i, "inverted_bathymetry_fname"] = (
        f"{ensemble_path}_inverted_bathymetry_{i}.nc"
    )
sampled_params_df.head()

[44]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[47]:

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_df = pd.read_csv(row.grav_df_fname)

    if row.grav_noise_levels == 0:
        l2_norm_tolerance = 0.2**0.5
    else:
        l2_norm_tolerance = row.grav_noise_levels**0.5

    # set kwargs to pass to the inversion
    kwargs = {
        # set stopping criteria
        "max_iterations": 200,
        "l2_norm_tolerance": l2_norm_tolerance,  # square root of the gravity noise
        "delta_l2_norm_tolerance": 1.008,
        "solver_damping": 0.025,
    }

    # run an optimization to find the optimal density contrast
    # automatically reruns inversion with optimal density contrast
    # and saves the results to <fname>_results.pickle
    _, _ = optimization.optimize_inversion_zref_density_contrast(
        grav_df=grav_df,
        constraints_df=constraint_points,
        density_contrast_limits=(1400, 3300),
        zref=0,
        n_trials=8,
        starting_topography=starting_bathymetry,
        regional_grav_kwargs=dict(
            method="constant",
            constant=0,
        ),
        fname=row.results_fname,
        plot_cv=False,
        progressbar=False,
        **kwargs,
    )

sampled_params_df.head()

Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.
Best density_contrast value (1400) is at the limit of provided values (3300, 1400) and thus is likely not a global minimum, expand the range of values tested to ensure the best parameter value is found.

[47]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[48]:

for i, row in sampled_params_df.iterrows():
    # load study
    with pathlib.Path(f"{row.results_fname}_study.pickle").open("rb") as f:
        study = pickle.load(f)

    # add best density to dataframe
    sampled_params_df.loc[i, "best_density_contrast"] = study.best_params[
        "density_contrast"
    ]

    # load saved inversion results
    with pathlib.Path(f"{row.results_fname}_results.pickle").open("rb") as f:
        topo_results, grav_results, parameters, elapsed_time = pickle.load(f)

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

    # sample the inverted topography at the constraint points
    constraint_points = utils.sample_grids(
        constraint_points,
        final_bathymetry,
        f"inverted_bathymetry_{i}",
        coord_names=("easting", "northing"),
    )
    constraints_rmse = utils.rmse(
        constraint_points.true_upward - constraint_points[f"inverted_bathymetry_{i}"]
    )

    # clip to inversion region
    final_bathymetry = final_bathymetry.sel(
        easting=slice(inversion_region[0], inversion_region[1]),
        northing=slice(inversion_region[2], inversion_region[3]),
    )

    inversion_rmse = utils.rmse(bathymetry - final_bathymetry)

    # save final topography to file
    final_bathymetry.to_netcdf(row.inverted_bathymetry_fname)

    # add to dataframe
    sampled_params_df.loc[i, "constraints_rmse"] = constraints_rmse
    sampled_params_df.loc[i, "inversion_rmse"] = inversion_rmse

sampled_params_df.head()

[48]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402

[49]:

sampled_params_df.to_csv(ensemble_fname, index=False)

[50]:

sampled_params_df = pd.read_csv(ensemble_fname)
sampled_params_df.head()

[50]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402

[51]:

fig, ax = plt.subplots()

for i, row in sampled_params_df.iterrows():
    # load study
    with pathlib.Path(f"{row.results_fname}_study.pickle").open("rb") as f:
        study = pickle.load(f)

    if i == 0:
        label = "Best"
        label_points = "Trials"
    else:
        label = None
        label_points = None

    study_df = study.trials_dataframe()
    # print(study_df)
    df = study_df[["value", "params_density_contrast"]].sort_values(
        by="params_density_contrast"
    )
    best_score = df.value.argmin()

    ax.plot(
        df.params_density_contrast.iloc[best_score],
        df.value.iloc[best_score],
        "s",
        markersize=6,
        color="r",
        label=label,
        zorder=10,
        markeredgecolor="black",
        linewidth=0.5,
    )
    ax.plot(
        df.params_density_contrast,
        df.value,
        # marker="o",
        color="b",
        # color=plt.cm.viridis(norm(row.grav_proximity)), # pylint: disable=no-member
    )
    ax.scatter(
        df.params_density_contrast,
        df.value,
        s=10,
        marker=".",
        color="black",
        edgecolors="black",
        zorder=10,
        label=label_points,
    )

ax.vlines(
    true_density_contrast,
    ax.get_ylim()[0],
    ax.get_ylim()[1],
    color="black",
    linestyle="--",
    label="True",
)
ax.legend(loc="best")
ax.set_title("Density contrast optimization")
ax.set_xlabel("Density contrast (kg/m$^3$)")
ax.set_ylabel("RMSE (m)")

# ax.set_xlim(true_density_contrast-100, true_density_contrast+100)
# ax.set_ylim(0, 100)

[51]:

Text(0, 0.5, 'RMSE (m)')

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_57_1.png

[52]:

fig = synth_plotting.plot_2var_ensemble(
    sampled_params_df,
    figsize=(6, 4),
    x="grav_line_spacing",
    y="grav_noise_levels",
    x_title="Line spacing (km)",
    y_title="Gravity noise (mGal)",
    background="best_density_contrast",
    background_title="Optimal density contrast",
    # background_robust=True,
    constrained_layout=False,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_58_0.png

Redo with optimal density¶

[53]:

for i, row in sampled_params_df.iterrows():
    # add results file name
    sampled_params_df.loc[i, "final_inversion_fname"] = (
        f"{ensemble_path}_final_inversion_{i}"
    )
    sampled_params_df.loc[i, "final_inverted_bathymetry_fname"] = (
        f"{ensemble_path}_final_inverted_bathymetry_{i}.nc"
    )

sampled_params_df.head()

[53]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse	final_inversion_fname	final_inverted_bathymetry_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[54]:

regional_grav_kwargs = dict(
    method="constant",
    constraints_df=constraint_points,
)

[55]:

logging.getLogger().setLevel(logging.WARNING)

for i, row in tqdm(sampled_params_df.iterrows(), total=len(sampled_params_df)):
    # load data
    grav_df = pd.read_csv(row.grav_df_fname)

    if row.grav_noise_levels == 0:
        l2_norm_tolerance = 0.2**0.5
    else:
        l2_norm_tolerance = row.grav_noise_levels**0.5

    # set kwargs to pass to the inversion
    kwargs = {
        # set stopping criteria
        "max_iterations": 200,
        "l2_norm_tolerance": l2_norm_tolerance,  # square root of the gravity noise
        "delta_l2_norm_tolerance": 1.008,
        "solver_damping": 0.025,
    }

    # run the inversion workflow using the optimal damping and density values
    _ = inversion.run_inversion_workflow(
        grav_df=grav_df,
        density_contrast=row.best_density_contrast,
        zref=0,
        starting_topography=starting_bathymetry,
        regional_grav_kwargs=regional_grav_kwargs,
        fname=row.final_inversion_fname,
        create_starting_prisms=True,
        calculate_starting_gravity=True,
        calculate_regional_misfit=True,
        **kwargs,
    )

sampled_params_df.head()

[55]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	best_filter_width	grav_data_loss_rmse	grav_data_loss_mae	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse	final_inversion_fname	final_inverted_bathymetry_fname
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	0.0	2.103457	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	16000.0	2.117755	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	28000.0	2.153094	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	32000.0	2.181365	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	36000.0	2.219025	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...

[ ]:

for i, row in sampled_params_df.iterrows():
    # load saved inversion results
    with pathlib.Path(f"{row.final_inversion_fname}_results.pickle").open("rb") as f:
        topo_results, grav_results, parameters, elapsed_time = pickle.load(f)

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

    # sample the inverted topography at the constraint points
    constraint_points = utils.sample_grids(
        constraint_points,
        final_bathymetry,
        f"final_inverted_bathymetry_{i}",
        coord_names=("easting", "northing"),
    )
    constraints_rmse = utils.rmse(
        constraint_points.true_upward
        - constraint_points[f"final_inverted_bathymetry_{i}"]
    )

    # clip to inversion region
    final_bathymetry = final_bathymetry.sel(
        easting=slice(inversion_region[0], inversion_region[1]),
        northing=slice(inversion_region[2], inversion_region[3]),
    )

    inversion_rmse = utils.rmse(bathymetry - final_bathymetry)

    # save final topography to file
    final_bathymetry.to_netcdf(row.final_inverted_bathymetry_fname)

    # add to dataframe
    sampled_params_df.loc[i, "final_inversion_constraints_rmse"] = constraints_rmse
    sampled_params_df.loc[i, "final_inversion_rmse"] = inversion_rmse

sampled_params_df["final_inversion_topo_improvement_rmse"] = (
    starting_bathymetry_rmse - sampled_params_df.final_inversion_rmse
)

sampled_params_df.head()

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	...	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse	final_inversion_fname	final_inverted_bathymetry_fname	final_inversion_constraints_rmse	final_inversion_rmse	final_inversion_topo_improvement_rmse
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48.115453	40.583114	-14.605765
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48.836411	41.348915	-15.371565
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	50.682037	42.094517	-16.117168
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	52.445568	42.719881	-16.742532
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	54.430651	43.495182	-17.517832

5 rows × 23 columns

[57]:

w = np.abs(inversion_region[1] - inversion_region[0]) / 1e3
h = np.abs(inversion_region[3] - inversion_region[2]) / 1e3
w, h

[57]:

(np.float64(150.0), np.float64(200.0))

[58]:

sampled_params_df["flight_kms"] = pd.Series()
for i, row in sampled_params_df.iterrows():
    survey_df = pd.read_csv(row.grav_survey_df_fname)
    # flight km is length of lines in km X number of lines
    sampled_params_df.loc[i, "flight_kms"] = (
        w * row.grav_line_numbers + h * row.grav_line_numbers
    )
sampled_params_df.flight_kms.unique()

[58]:

array([np.float64(1050.0), np.float64(1400.0), np.float64(1750.0),
       np.float64(2100.0), np.float64(2800.0), np.float64(3850.0),
       np.float64(4900.0), np.float64(6300.0), np.float64(8050.0),
       np.float64(10500.0)], dtype=object)

[59]:

inversion_area = (
    (inversion_region[1] - inversion_region[0])
    / 1e3
    * (inversion_region[3] - inversion_region[2])
    / 1e3
)
sampled_params_df["flight_kms_per_10000sq_km"] = (
    (sampled_params_df.flight_kms / inversion_area) * 10e3
).astype(float)
sampled_params_df.flight_kms_per_10000sq_km.unique()

[59]:

array([ 350.        ,  466.66666667,  583.33333333,  700.        ,
        933.33333333, 1283.33333333, 1633.33333333, 2100.        ,
       2683.33333333, 3500.        ])

[13]:

sampled_params_df.to_csv(ensemble_fname, index=False)

[11]:

sampled_params_df = pd.read_csv(ensemble_fname)
sampled_params_df.head()

[11]:

	grav_line_numbers	grav_noise_levels	grav_df_fname	grav_survey_df_fname	grav_line_spacing	grav_proximity	grav_proximity_skew	flight_kms	flight_kms_per_10000sq_km	filter_width_trials_fname	...	grav_data_loss_mae	results_fname	inverted_bathymetry_fname	best_density_contrast	constraints_rmse	inversion_rmse	final_inversion_fname	final_inverted_bathymetry_fname	final_inversion_constraints_rmse	final_inversion_rmse
0	3	0.00	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	1.163586	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1400.0	56.844169	50.830768	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48.115453	40.583114
1	3	0.56	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	1.203714	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1474.0	48.649062	41.422841	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	48.836411	41.348915
2	3	1.11	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	1.262860	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1478.0	50.601559	42.122256	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	50.682037	42.094517
3	3	1.67	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	1.313212	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1482.0	52.395224	42.720566	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	52.445568	42.719881
4	3	2.22	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	58.333333	8.333333	0.489544	1050.0	350.0	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	...	1.367866	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	1486.0	54.442471	43.497402	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	../results/Ross_Sea/ensembles/Ross_Sea_ensembl...	54.430651	43.495182

5 rows × 22 columns

Examine results¶

[62]:

# subplots showing inverted bathymetry error
grids = []
cbar_labels = []
row_titles = []
column_titles = []

# iterate over the ensemble starting with high noise, low line numbers
# row per noise level and column per line spacing
for i, row in sampled_params_df.sort_values(
    by=["grav_noise_levels", "grav_line_spacing"],
    ascending=False,
).iterrows():
    dif = bathymetry - xr.open_dataarray(row.final_inverted_bathymetry_fname)
    # add to lists
    grids.append(dif)
    cbar_labels.append(f"RMSE: {round(utils.rmse(dif), 2)} (mGal)")
    if i % num == 0:
        column_titles.append(f"Spacing: {round(row.grav_line_spacing)} km")
    if i < num:
        row_titles.append(f"Noise: {round(row.grav_noise_levels, 1)} mGal")

fig = maps.subplots(
    grids,
    fig_height=8,
    fig_title="Inverted bathymetry error",
    fig_title_font="75p,Helvetica-Bold,black",
    fig_title_y_offset="5c",
    cbar_labels=cbar_labels,
    cmap="balance+h0",
    # hist=True,
    cpt_lims=polar_utils.get_combined_min_max(
        grids,
        robust=True,
        absolute=True,
    ),
    cbar_font="18p,Helvetica,black",
    row_titles=row_titles,
    column_titles=column_titles,
    row_titles_font="25p,Helvetica,black",
    column_titles_font="25p,Helvetica,black",
    yshift_amount=-1.2,
)
fig.show()

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_70_0.png

[63]:

fig = synth_plotting.plot_2var_ensemble(
    sampled_params_df,
    x="grav_line_spacing",
    y="grav_noise_levels",
    x_title="Line spacing (km)",
    y_title="Gravity noise (mGal)",
    background="final_inversion_rmse",
    background_title="Bathymetry RMSE (m)",
    background_robust=True,
    plot_contours=[starting_bathymetry_rmse],
    plot_title="with optimal density",
    constrained_layout=False,
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_72_0.png

[64]:

_ = synth_plotting.plot_ensemble_as_lines(
    sampled_params_df,
    y="final_inversion_rmse",
    x="grav_noise_levels",
    groupby_col="grav_line_spacing",
    x_label="Gravity noise (mGal)",
    cbar_label="Line spacing (km)",
    trend_line=True,
    horizontal_line=starting_bathymetry_rmse,
    horizontal_line_label="Starting RMSE",
    horizontal_line_label_loc="lower right",
    plot_title="No regional field",
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_73_0.png

[65]:

_ = synth_plotting.plot_ensemble_as_lines(
    sampled_params_df,
    y="final_inversion_rmse",
    x="grav_line_spacing",
    groupby_col="grav_noise_levels",
    x_label="Line spacing (km)",
    cbar_label="Gravity noise (mGal)",
    trend_line=True,
    horizontal_line=starting_bathymetry_rmse,
    horizontal_line_label="Starting RMSE",
    horizontal_line_label_loc="lower right",
    plot_title="No regional field",
)

_images/ensemble_02_grav_spacing_vs_noise_no_regional_density_estimation_74_0.png