Optimal buffer zone width

Optimal buffer zone width#

To limit gravity edge effects, many inversions use a buffer zone to extend the model beyond the extent of the gravity data. This limits the edge effects to the portion of the model outside the area of interest. A large buffer zone results is smaller errors in the area of interest, but the larger model can lead to significantly slower runtime of the code. Here we show a method to choose the optimal buffer zone width to limit edge effects, but not significantly slow down the inversion.

The edge effects are due to the contrast between the edge of the model and the void space beyond. The larger and/or denser the prisms near the edge, the large the contrast with the void space will be and thus the larger the edge effects will be.

Import packages#

[1]:

# set EPSG for plotting functions
import os

import polartoolkit as ptk
import verde as vd
import xrft

import invert4geom

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

/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/UQpy/__init__.py:6: UserWarning:

pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.

Define model parameters#

[2]:

spacing = 1000
region = (0, 40000, 0, 30000)
zref = 0
density_contrast = 2670
gravity_obs_height = 1e3

Quantify edge effects#

Here we show a technique for trying to quantify edge effects. We create a prism layer which represents a simplistic version of our true prism layer which will be used in the inversion. This simple layer of prisms will have flat tops and bottoms, and constant density values. These values should represent the mean values of the true prism model. We then calculate the forward gravity effect of these prisms across an inner_region which represents our area of interest. A buffer zone can be set as a percentage of region width. Since the prisms are flat, the maximum gravity value will be in the center of the prism layer and the gravity values will decay towards the edges of the region of interest as the edge effect increases. The gravity decay can be estimated as the percentile difference between the maximum and minimum gravity values in the region of interest.

To represent our prism layer we want to use in the inversion (from above cells), we use the same zref, the mean topography for the top values, and the same density contrast values.

[7]:

max_decay, _, _, _ = invert4geom.gravity_decay_buffer(
    buffer_perc=0,
    obs_height=gravity_obs_height,
    top=topography.values.mean(),
    zref=zref,
    spacing=spacing,
    inner_region=region,
    density=density_contrast,
)

../_images/how_to_determining_buffer_width_11_0.png

../_images/how_to_determining_buffer_width_11_1.png

../_images/how_to_determining_buffer_width_11_2.png

[8]:

print(f"Max decay: {round(max_decay, 1)}%")

Max decay: 49.3%

As you can see above, the technique for estimating the gravity decay (edge effects) works pretty well. For this scenario, with no buffer zone, the maximum amount of decay is 50%. This would introduce significant errors into the inversion. Below we demonstrate how to choose the optimal buffer zone width.

Find optimal buffer zone width#

We can run an optimization on a range of buffer zone widths and find the value which gets us closest to the desired amount of gravitational decay. This set amount of decay is chosen with target. High values will result in a small buffer zone (inversions will run faster), but will include noticable edge effects at the model edge. Conversely, small values will results in large buffer zones (slowing down the code), but will limit edge effects. Here we use a target of 3% decay. The optimization will test 20 values of buffer zone widths between 1% and 100% of the region width.

Do estimate the gravity decay, we create a flat layer of prisms and calculate it’s forward gravity. The max value is then compared to the edge values to determine the amount of decay. Use the parameters inner_region, spacing, top, zref and density to mimic your inversion model as closely as possible.

[9]:

study, results = invert4geom.optimal_buffer(
    target=3,
    buffer_perc_limits=(1, 100),
    n_trials=20,
    obs_height=gravity_obs_height,
    top=topography.values.mean(),
    zref=zref,
    spacing=spacing,
    inner_region=region,
    density=density_contrast,
)
_, buffer_width, buffer_cells, _ = results

../_images/how_to_determining_buffer_width_15_1.png

../_images/how_to_determining_buffer_width_15_2.png

../_images/how_to_determining_buffer_width_15_3.png

[10]:

print(f"Optimal buffer width: {buffer_width / 1e3} km")
print(f"Optimal number of buffer cells: {buffer_cells}")

Optimal buffer width: 7.0 km
Optimal number of buffer cells: 7

Use the optimal buffer width#

We will recalculate the forward gravity of the topography, but adding a buffer zone with the width determined above. If you’re using real topography data, you can just use a grid which is larger by the amount of the buffer zone instead of padding the grid.

[11]:

# pad the topography by with optimal buffer width
pad_width = {
    "northing": buffer_cells,
    "easting": buffer_cells,
}
buffer_topography = xrft.pad(
    topography,
    pad_width,
    mode="linear_ramp",
    constant_values=None,
    end_values=topography.median(),
)

[12]:

fig = ptk.plot_grid(
    buffer_topography,
    show_region=region,
    cbar_label="elevation (m)",
    scalebar=True,
    hist=True,
)
fig.show()

../_images/how_to_determining_buffer_width_19_0.png

[13]:

buffer_model = invert4geom.create_model(
    zref=zref,
    density_contrast=density_contrast,
    topography=buffer_topography.to_dataset(name="upward"),
)

# calculate gravity anomaly
grav_data.inv.forward_gravity(buffer_model, "gravity_anomaly_buffer")

[14]:

_ = ptk.grid_compare(
    grav_data.gravity_anomaly,
    grav_data.gravity_anomaly_buffer,
    title="Forward gravity comparison",
    grid1_name="No buffer",
    grid2_name=f"{buffer_width / 1e3} km buffer",
    fig_height=8,
    hist=True,
    inset=False,
    coast=False,
    robust=True,
)

../_images/how_to_determining_buffer_width_21_0.png

[15]:

# profile start and stop locations
start = (region[0], (region[3] - region[2]) / 2)
stop = (region[1], (region[3] - region[2]) / 2)

cpt_lims = ptk.get_combined_min_max(
    [grav_data.gravity_anomaly, grav_data.gravity_anomaly_buffer]
)

# plot topography
fig = ptk.plot_grid(
    grav_data.gravity_anomaly,
    hist=True,
    cpt_lims=cpt_lims,
    title="Forward gravity no buffer",
    cbar_label="Gravity anomaly (mGal)",
)
# plot profile location, and endpoints on map
fig.plot(
    x=[start[0], stop[0]],
    y=[start[1], stop[1]],
    pen="2p,red",
)
fig.text(
    x=start[0],
    y=start[1],
    text="A",
    fill="white",
    font="12p,Helvetica,black",
    justify="CM",
    clearance="+tO",
    no_clip=True,
)
fig.text(
    x=stop[0],
    y=stop[1],
    text="B",
    fill="white",
    font="12p,Helvetica,black",
    justify="CM",
    clearance="+tO",
    no_clip=True,
)

# plot gravity anomaly
fig = ptk.plot_grid(
    grav_data.gravity_anomaly_buffer,
    fig=fig,
    origin_shift="x",
    cpt_lims=cpt_lims,
    hist=True,
    title="Forward gravity with buffer",
    cbar_label="Gravity anomaly (mGal)",
)
# plot profile location, and endpoints on map
fig.plot(
    x=[start[0], stop[0]],
    y=[start[1], stop[1]],
    pen="2p,red",
)
fig.text(
    x=start[0],
    y=start[1],
    text="A",
    fill="white",
    font="12p,Helvetica,black",
    justify="CM",
    clearance="+tO",
    no_clip=True,
)
fig.text(
    x=stop[0],
    y=stop[1],
    text="B",
    fill="white",
    font="12p,Helvetica,black",
    justify="CM",
    clearance="+tO",
    no_clip=True,
)
fig.show()

../_images/how_to_determining_buffer_width_22_0.png

[16]:

# plot profile across region
data_dict = ptk.make_data_dict(
    ["forward gravity no buffer", "forward gravity with buffer"],
    [grav_data.gravity_anomaly, grav_data.gravity_anomaly_buffer],
    ["red", "blue"],
)

layers_dict = ptk.make_data_dict(
    ["prism tops", "prism bottoms"],
    [buffer_model.top, buffer_model.bottom],
    ["blue", "darkorange"],
)

fig, _, _ = ptk.plot_profile(
    "points",
    add_map=True,
    map_background=grav_data.gravity_anomaly,
    subplot_orientation="vertical",
    coast=False,
    inset=False,
    gridlines=False,
    fig_width=10,
    fig_height=8,
    data_height=6,
    start=start,
    stop=stop,
    layers_dict=layers_dict,
    data_dict=data_dict,
    fill_layers=False,
)
fig.show()

../_images/how_to_determining_buffer_width_23_0.png

[ ]: