feat: add DC-shift regional estimation method

src/invert4geom/regional.py

@@ -13,6 +13,61 @@
 from invert4geom import optimization, utils
 
 
+def regional_dc_shift(
+    grav_df: pd.DataFrame,
+    dc_shift: float | None = None,
+    grav_grid: xr.DataArray | None = None,
+    constraint_points: pd.DataFrame | None = None,
+    coord_names: tuple[str, str] = ("easting", "northing"),
+    regional_col_name: str = "reg",
+) -> pd.DataFrame:
+    """
+    separate the regional field by applying a constant shift (DC-shift) to the gravity
+    data. If constraint points of the layer of interested are supplied, the DC shift
+    will minimize the residual misfit at these constraint points.
+
+    Parameters
+    ----------
+    grav_df : pd.DataFrame
+        gravity data with columns defined by coord_names and input_grav_name.
+    dc_shift : float
+        shift to apply to the data
+    grav_grid : xr.DataArray
+        gridded gravity misfit data
+    constraint_points : pd.DataFrame
+        a dataframe of constraint points with columns X and Y columns defined by the
+        coord_names parameter.
+    coord_names : tuple
+        names of the X and Y column names in constraint points dataframe
+    regional_col_name : str
+        name for the new column in grav_df for the regional field.
+
+    Returns
+    -------
+    pd.DataFrame
+        grav_df with new regional column
+    """
+    if constraint_points is not None:
+        # get the gravity values at the constraint points
+        constraints_df = constraint_points.copy()
+
+        # sample gravity at constraint points
+        constraints_df = utils.sample_grids(
+            df=constraints_df,
+            grid=grav_grid,
+            sampled_name="sampled_grav",
+            coord_names=coord_names,
+        )
+
+        # use RMS of sampled value for DC shift
+        dc_shift = utils.rmse(constraints_df.sampled_grav)
+
+    grav_df[regional_col_name] = dc_shift
+
+    # return the new dataframe
+    return grav_df
+
+
 def regional_filter(
     filter_width: float,
     grav_grid: xr.DataArray,

tests/test_regional.py

@@ -39,6 +39,52 @@ def dummy_df() -> pd.DataFrame:
 
 
 # %%
+def test_regional_dc_shift_constraints():
+    """
+    test the regional_dc_shift function with a supplied constraints
+    """
+    grav_df = dummy_df()
+    region = (0, 200, 200, 400)
+
+    # create 10 random point within the region
+    num_constraints = 10
+    coords = vd.scatter_points(region=region, size=num_constraints, random_state=0)
+    points = pd.DataFrame(data={"easting": coords[0], "northing": coords[1]})
+
+    # print(grav_df.describe())
+
+    df = regional.regional_dc_shift(
+        grav_df=grav_df,
+        grav_grid=dummy_grid().misfit,
+        constraint_points=points,
+        regional_col_name="reg",
+    )
+
+    # print(df.describe())
+
+    # test whether regional field has been removed correctly
+    # by whether the means of the reg and misfit are similar
+    # print(np.mean(df.reg), np.mean(df.misfit))
+    print(np.mean(df.reg), np.mean(df.misfit))
+    assert np.mean(df.reg) == pytest.approx(np.mean(df.misfit), rel=1000)
+
+
+def test_regional_dc_shift():
+    """
+    test the regional_dc_shift function with a supplied DC shift
+    """
+
+    grav_df = dummy_grid().to_dataframe().reset_index()
+
+    df = regional.regional_dc_shift(
+        grav_df=grav_df,
+        dc_shift=-200,
+        regional_col_name="reg",
+    )
+
+    assert df.reg.mean() == -200
+
+
 @pytest.mark.parametrize("fill_method", ["rioxarray", "verde"])
 @pytest.mark.parametrize("trend", [0, 2])
 def test_regional_trend(fill_method, trend):