Estimating a regional field#

Here we will present the 4 methods we provide for estimating the regional component of gravity.

Import packages#

[1]:

from __future__ import annotations

%load_ext autoreload
%autoreload 2

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

from invert4geom import plotting, regional, synthetic, utils

Create observed gravity data#

Create 2 synthetic surfaces, an upper one with shorter wavelength topographic features, and a lower once with longer wavelength features. The upper surface will produce the gravity component of interest, the residual, and the lower surface will produce the gravity component we aim to estimate and remove, the regional. We will create prism layers for both of these layers, forward model them onto a set of gravity observation points. The combined forward gravity effect of these layes, with some added noise, will represent the observed gravity.

[2]:

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

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

# create lower synthetic topography data
lower_topography = synthetic.synthetic_topography_regional(
    spacing,
    region,
    scale=3,
    yoffset=-1000,
)

[3]:

# create constraint points by sampling upper surface at random locations
num_constraints = 30
coords = vd.scatter_points(region=region, size=num_constraints, random_state=7)
constraint_points = pd.DataFrame(data={"easting": coords[0], "northing": coords[1]})

# sample true topography at these points
constraint_points = utils.sample_grids(
    constraint_points, upper_topography, "upward", coord_names=("easting", "northing")
)

Prism layers#

[14]:

# the density contrast is between sediment (~1800 kg/m3) and air (~1 kg/m3)
density_contrast = 1800

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

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

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

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

# to add an offset, prisms are created between -2km the height of the topography
zref = -2e3  # lower_topography.values.mean()

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

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

plotting.show_prism_layers(
    [upper_prisms, lower_prisms],
    color_by="thickness",
    log_scale=False,
    zscale=20,
    backend="static",
)

../_images/user_guide_estimating_regional_field_7_0.png

Forward gravity of prism layers#

[15]:

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

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

grav_df = vd.grid_to_table(observations)

# forward gravity of upper surface prisms
grav_df["upper_surface_grav"] = upper_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=True,
)
# forward gravity of lower surface prisms
grav_df["lower_surface_grav"] = lower_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=True,
)
grav_df["grav"] = grav_df.upper_surface_grav + grav_df.lower_surface_grav
grav_df

[15]:

	northing	easting	upward	upper_surface_grav	lower_surface_grav	grav
0	0.0	0.0	1000.0	4.702821	6.858162	11.560984
1	0.0	500.0	1000.0	5.761512	7.729221	13.490732
2	0.0	1000.0	1000.0	5.940439	8.493371	14.433810
3	0.0	1500.0	1000.0	5.886537	9.126971	15.013508
4	0.0	2000.0	1000.0	5.748903	9.635873	15.384776
...	...	...	...	...	...	...
4936	30000.0	38000.0	1000.0	1.911976	12.415943	14.327919
4937	30000.0	38500.0	1000.0	1.900899	11.744157	13.645055
4938	30000.0	39000.0	1000.0	1.866645	10.909692	12.776336
4939	30000.0	39500.0	1000.0	1.760582	9.901674	11.662256
4940	30000.0	40000.0	1000.0	1.391119	8.749302	10.140421

4941 rows × 6 columns

[16]:

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

# grid the contaminated gravity
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

grav_grid.observed_grav.plot()
grav_df

[16]:

	northing	easting	upward	upper_surface_grav	lower_surface_grav	grav	observed_grav
0	0.0	0.0	1000.0	4.702821	6.858162	11.560984	11.626785
1	0.0	500.0	1000.0	5.761512	7.729221	13.490732	13.427617
2	0.0	1000.0	1000.0	5.940439	8.493371	14.433810	14.756958
3	0.0	1500.0	1000.0	5.886537	9.126971	15.013508	15.068894
4	0.0	2000.0	1000.0	5.748903	9.635873	15.384776	15.119878
...	...	...	...	...	...	...	...
4936	30000.0	38000.0	1000.0	1.911976	12.415943	14.327919	15.208171
4937	30000.0	38500.0	1000.0	1.900899	11.744157	13.645055	13.156471
4938	30000.0	39000.0	1000.0	1.866645	10.909692	12.776336	13.497419
4939	30000.0	39500.0	1000.0	1.760582	9.901674	11.662256	10.902134
4940	30000.0	40000.0	1000.0	1.391119	8.749302	10.140421	9.877312

4941 rows × 7 columns

../_images/user_guide_estimating_regional_field_10_1.png

Regional estimation methods#

DC-Shift#

[24]:

# estimate regional with a 30km low pass filter
grav_df = regional.regional_dc_shift(
    grav_grid=grav_grid.observed_grav,
    grav_df=grav_df,
    constraint_points=constraint_points,
    regional_col_name="dc_shift_reg",
)

# grid the results
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

# compare with true regional
_ = polar_utils.grd_compare(
    grav_grid.lower_surface_grav,
    grav_grid.dc_shift_reg,
    plot_type="xarray",
    robust=True,
    plot=True,
    grid1_name="True regional field",
    grid2_name="Regional field from method: DC-Shift",
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
)

../_images/user_guide_estimating_regional_field_13_0.png

Filter#

[18]:

# estimate regional with a 30km low pass filter
grav_df = regional.regional_filter(
    filter_width=30e3,
    grav_grid=grav_grid.observed_grav,
    grav_df=grav_df,
    regional_col_name="filter_reg",
)

# grid the results
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

# compare with true regional
_ = polar_utils.grd_compare(
    grav_grid.lower_surface_grav,
    grav_grid.filter_reg,
    plot_type="xarray",
    robust=True,
    plot=True,
    grid1_name="True regional field",
    grid2_name="Regional field from method: Filter",
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
)

../_images/user_guide_estimating_regional_field_15_0.png

Trend#

[19]:

# estimate regional with fitting a 4th order trend
grav_df = regional.regional_trend(
    trend=4,
    grav_grid=grav_grid.observed_grav,
    grav_df=grav_df,
    fill_method="verde",
    regional_col_name="trend_reg",
)

# grid the results
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

# compare with true regional
_ = polar_utils.grd_compare(
    grav_grid.lower_surface_grav,
    grav_grid.trend_reg,
    plot_type="xarray",
    robust=True,
    plot=True,
    grid1_name="True regional field",
    grid2_name="Regional field from method: Trend",
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
)

../_images/user_guide_estimating_regional_field_17_0.png

Equivalent Sources#

[20]:

# estimate regional with fitting deep equivalent sources
grav_df = regional.regional_eq_sources(
    source_depth=-100e3,
    grav_df=grav_df,
    input_grav_name="observed_grav",
    eq_damping=None,
    block_size=2e3,
    regional_col_name="eq_sources_reg",
)

# grid the results
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

# compare with true regional
_ = polar_utils.grd_compare(
    grav_grid.lower_surface_grav,
    grav_grid.eq_sources_reg,
    plot_type="xarray",
    robust=True,
    plot=True,
    grid1_name="True regional field",
    grid2_name="Regional field from method: Equivalent sources",
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
)

../_images/user_guide_estimating_regional_field_19_0.png

Constraint point minimization#

[21]:

# estimate regional with the constraints method
grav_df = regional.regional_constraints(
    constraint_points=constraint_points,
    grav_grid=grav_grid.observed_grav,
    grav_df=grav_df,
    region=region,
    spacing=spacing,
    regional_col_name="constraints_reg",
    grid_method="pygmt",
    tension_factor=1,
)

# grid the results
grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()

# compare with true regional
_ = polar_utils.grd_compare(
    grav_grid.lower_surface_grav,
    grav_grid.constraints_reg,
    plot_type="xarray",
    robust=True,
    plot=True,
    grid1_name="True regional field",
    grid2_name="Regional field from method: Constraints",
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    grounding_line=False,
    points=constraint_points.rename(columns={"easting": "x", "northing": "y"}),
)

../_images/user_guide_estimating_regional_field_21_0.png

[23]:

names = ["dc_shift_reg", "filter_reg", "trend_reg", "eq_sources_reg", "constraints_reg"]

fig = maps.subplots(
    [grav_grid[n] for n in names],
    region,
    margins="2c",
    fig_width=20,
    # grd2cpt=True,
    fig_title="Regional components",
    subplot_titles=names,
    cbar_units=["m"] * len(names),
    autolabel=True,
)

fig.show()

makecpt [ERROR]: Option T: min >= max
ERROR:root:Module 'makecpt' failed with status code 72:
makecpt [ERROR]: Option T: min >= max
Traceback (most recent call last):
  File "/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/polartoolkit/maps.py", line 432, in plot_grd
    pygmt.makecpt(
  File "/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/pygmt/helpers/decorators.py", line 603, in new_module
    return module_func(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/pygmt/helpers/decorators.py", line 776, in new_module
    return module_func(*bound.args, **bound.kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/pygmt/src/makecpt.py", line 165, in makecpt
    lib.call_module(module="makecpt", args=arg_str)
  File "/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/pygmt/clib/session.py", line 624, in call_module
    raise GMTCLibError(
pygmt.exceptions.GMTCLibError: Module 'makecpt' failed with status code 72:
makecpt [ERROR]: Option T: min >= max
WARNING:root:grid region can't be extracted.

../_images/user_guide_estimating_regional_field_22_1.png

The above plot and code show four techniques to estimate the regional field from gravity data. Each of these techniques requires the user input to pick the value of a certain parameter, each of which has a large effect on the outcome. Below these hyperparameters for each technique are explained:

Filter technique: the width of the gaussian filter

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

Equivalent Sources technique: the source depth and optionally a damping parameter

Constraint Point Minimization technique: the tension factor if using minimum curvature gridding, or the damping value if using bi-harmonic spline gridding

Above we just picked these hyperparameter values arbitrarily. Below we show a more informed technique for choosing each of these hyperparameter values.

[1]:

# STILL TO COME

[ ]:

Estimating a regional field

Contents