Levelling to a grid

1. Levelling to a grid#

[1]:

%load_ext autoreload
%autoreload 2

import boule
import numpy as np
import pandas as pd
import plotly.io as pio
import polartoolkit as ptk
import verde as vd

import airbornegeo

pio.renderers.default = "notebook"

/home/airbornegeo/airbornegeo/.pixi/envs/default/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

1.1. Load survey data#

This is a subset of the BAS AGAP survey over Antarctica’s Gamburtsev Subglacial Mountains. The file is download and subset in the notebook AGAP_magnetic_survey, and the BAS processing steps are repeated in the notebook processing_AGAP_magnetic_survey.

[2]:

data_df = pd.read_csv("data/AGAP_gravity_survey_processed.csv")
data_df = data_df[
    [
        "easting",
        "northing",
        "height",
        "line",
        "unixtime",
        "distance_along_line",
        "grav_disturbance_filt",
    ]
]
data_df.head()

[2]:
easting northing height line unixtime distance_along_line grav_disturbance_filt
0 1.000024e+06 226237.330771 4156.1 1 1.229507e+09 0.000000 49.38
1 1.000083e+06 226246.631269 4156.0 1 1.229507e+09 59.842447 49.45
2 1.000142e+06 226255.809132 4156.1 1 1.229507e+09 119.693401 49.52
3 1.000201e+06 226264.969079 4156.4 1 1.229507e+09 179.545645 49.58
4 1.000260e+06 226274.156809 4156.6 1 1.229507e+09 239.285174 49.65 
[3]:

# plot the data
airbornegeo.plotly_points(
    data_df[::20],
    color_col="grav_disturbance_filt",
    robust=True,
    size=2,
)

1.2. Get a grid of gravity disturbance at 10km#

[4]:

# Download satellite gravity data from EIGEN-6C4 model
region = vd.get_region((data_df.easting, data_df.northing))
eigen = ptk.fetch.gravity(
    version="eigen",
    spacing=5e3,
    region=region,
)
eigen_df = vd.grid_to_table(eigen)
eigen_df = eigen_df.rename(columns={"x": "easting", "y": "northing"})

# reproject from EPSG 3031 to lat lon
eigen_df["lon"], eigen_df["lat"] = airbornegeo.reproject(
    eigen_df.easting,
    eigen_df.northing,
    input_crs="EPSG:3031",
    output_crs="EPSG:4326",
)

# calculated normal gravity at all EIGEN observation locations
eigen_df["normal_gravity"] = boule.WGS84.normal_gravity(
    latitude=eigen_df.lat,
    height=eigen_df.ellipsoidal_height,
)

# calculate gravity disturbance
eigen_df["disturbance"] = eigen_df.gravity - eigen_df.normal_gravity

# convert to a dataset
eigen_ds = eigen_df.set_index(["northing", "easting"]).to_xarray()

eigen_ds

[4]:

<xarray.Dataset> Size: 479kB
Dimensions:             (northing: 94, easting: 127)
Coordinates:
  * northing            (northing) float64 752B 1.2e+05 1.25e+05 ... 5.85e+05
  * easting             (easting) float64 1kB 1e+06 1.005e+06 ... 1.63e+06
Data variables:
    ellipsoidal_height  (northing, easting) float32 48kB 1e+04 1e+04 ... 1e+04
    gravity             (northing, easting) float32 48kB 9.8e+05 ... 9.797e+05
    lon                 (northing, easting) float64 96kB 83.16 83.19 ... 70.26
    lat                 (northing, easting) float64 96kB -80.75 -80.7 ... -74.16
    normal_gravity      (northing, easting) float64 96kB 9.8e+05 ... 9.798e+05
    disturbance         (northing, easting) float64 96kB 38.05 38.12 ... -20.55
[5]:

eigen_ds.disturbance.plot()

[5]:

<matplotlib.collections.QuadMesh at 0x7fc9e078cad0>
_images/levelling_to_a_grid_7_1.png

1.3. Upward continue survey data to same height as gravity grid#

In order to accurately compare the survey gravity data to the grid, the gravity data should be upward continued so it’s at the same altitude.

[6]:

# plot the data
airbornegeo.plotly_points(
    data_df[::20],
    color_col="height",
    robust=True,
    size=2,
)

[7]:

blocked_survey = airbornegeo.block_reduce(
    data_df,
    np.median,
    spacing=1000,
    reduce_by="distance_along_line",
    groupby_column="line",
)
blocked_survey

Segments: 100%|██████████| 100/100 [00:04<00:00, 24.14it/s]

[7]:
distance_along_line easting northing height unixtime grav_disturbance_filt line
0 477.521406 1.000496e+06 226310.158049 4159.1 1.229507e+09 49.89 1
1 1477.462561 1.001484e+06 226460.426302 4160.3 1.229507e+09 50.73 1
2 2504.019094 1.002497e+06 226630.494762 4156.0 1.229507e+09 51.09 1
3 3518.522095 1.003499e+06 226786.791629 4159.6 1.229507e+09 51.07 1
4 4528.427996 1.004498e+06 226936.380167 4165.7 1.229507e+09 51.43 1
... ... ... ... ... ... ... ...
21373 68170.289011 1.586103e+06 507079.737473 2049.2 1.230382e+09 -1.66 100
21374 69150.087140 1.586259e+06 506112.838506 2051.8 1.230382e+09 -2.58 100
21375 70159.126161 1.586438e+06 505120.357453 2053.1 1.230382e+09 -4.03 100
21376 71178.859535 1.586603e+06 504114.078250 2058.5 1.230382e+09 -5.99 100
21377 72178.780154 1.586775e+06 503129.075505 2072.0 1.230382e+09 -8.16 100

21378 rows × 7 columns

[10]:

# fit a set of equivalent sources to each line individually
eqs = airbornegeo.eq_sources_1d(
    blocked_survey,
    data_column="grav_disturbance_filt",
    depth="default",
    damping=None,
    block_size=1000,  # for speed, block reduce sources
    groupby_column="line",
)
eqs

Groups: 100%|██████████| 100/100 [00:04<00:00, 24.49it/s]

[10]:

{1: EquivalentSources(block_size=1000),
 2: EquivalentSources(block_size=1000),
 3: EquivalentSources(block_size=1000),
 4: EquivalentSources(block_size=1000),
 5: EquivalentSources(block_size=1000),
 6: EquivalentSources(block_size=1000),
 7: EquivalentSources(block_size=1000),
 8: EquivalentSources(block_size=1000),
 9: EquivalentSources(block_size=1000),
 10: EquivalentSources(block_size=1000),
 11: EquivalentSources(block_size=1000),
 12: EquivalentSources(block_size=1000),
 13: EquivalentSources(block_size=1000),
 14: EquivalentSources(block_size=1000),
 15: EquivalentSources(block_size=1000),
 16: EquivalentSources(block_size=1000),
 17: EquivalentSources(block_size=1000),
 18: EquivalentSources(block_size=1000),
 19: EquivalentSources(block_size=1000),
 20: EquivalentSources(block_size=1000),
 21: EquivalentSources(block_size=1000),
 22: EquivalentSources(block_size=1000),
 23: EquivalentSources(block_size=1000),
 24: EquivalentSources(block_size=1000),
 25: EquivalentSources(block_size=1000),
 26: EquivalentSources(block_size=1000),
 27: EquivalentSources(block_size=1000),
 28: EquivalentSources(block_size=1000),
 29: EquivalentSources(block_size=1000),
 30: EquivalentSources(block_size=1000),
 31: EquivalentSources(block_size=1000),
 32: EquivalentSources(block_size=1000),
 33: EquivalentSources(block_size=1000),
 34: EquivalentSources(block_size=1000),
 35: EquivalentSources(block_size=1000),
 36: EquivalentSources(block_size=1000),
 37: EquivalentSources(block_size=1000),
 38: EquivalentSources(block_size=1000),
 39: EquivalentSources(block_size=1000),
 40: EquivalentSources(block_size=1000),
 41: EquivalentSources(block_size=1000),
 42: EquivalentSources(block_size=1000),
 43: EquivalentSources(block_size=1000),
 44: EquivalentSources(block_size=1000),
 45: EquivalentSources(block_size=1000),
 46: EquivalentSources(block_size=1000),
 47: EquivalentSources(block_size=1000),
 48: EquivalentSources(block_size=1000),
 49: EquivalentSources(block_size=1000),
 50: EquivalentSources(block_size=1000),
 51: EquivalentSources(block_size=1000),
 52: EquivalentSources(block_size=1000),
 53: EquivalentSources(block_size=1000),
 54: EquivalentSources(block_size=1000),
 55: EquivalentSources(block_size=1000),
 56: EquivalentSources(block_size=1000),
 57: EquivalentSources(block_size=1000),
 58: EquivalentSources(block_size=1000),
 59: EquivalentSources(block_size=1000),
 60: EquivalentSources(block_size=1000),
 61: EquivalentSources(block_size=1000),
 62: EquivalentSources(block_size=1000),
 63: EquivalentSources(block_size=1000),
 64: EquivalentSources(block_size=1000),
 65: EquivalentSources(block_size=1000),
 66: EquivalentSources(block_size=1000),
 67: EquivalentSources(block_size=1000),
 68: EquivalentSources(block_size=1000),
 69: EquivalentSources(block_size=1000),
 70: EquivalentSources(block_size=1000),
 71: EquivalentSources(block_size=1000),
 72: EquivalentSources(block_size=1000),
 73: EquivalentSources(block_size=1000),
 74: EquivalentSources(block_size=1000),
 75: EquivalentSources(block_size=1000),
 76: EquivalentSources(block_size=1000),
 77: EquivalentSources(block_size=1000),
 78: EquivalentSources(block_size=1000),
 79: EquivalentSources(block_size=1000),
 80: EquivalentSources(block_size=1000),
 81: EquivalentSources(block_size=1000),
 82: EquivalentSources(block_size=1000),
 83: EquivalentSources(block_size=1000),
 84: EquivalentSources(block_size=1000),
 85: EquivalentSources(block_size=1000),
 86: EquivalentSources(block_size=1000),
 87: EquivalentSources(block_size=1000),
 88: EquivalentSources(block_size=1000),
 89: EquivalentSources(block_size=1000),
 90: EquivalentSources(block_size=1000),
 91: EquivalentSources(block_size=1000),
 92: EquivalentSources(block_size=1000),
 93: EquivalentSources(block_size=1000),
 94: EquivalentSources(block_size=1000),
 95: EquivalentSources(block_size=1000),
 96: EquivalentSources(block_size=1000),
 97: EquivalentSources(block_size=1000),
 98: EquivalentSources(block_size=1000),
 99: EquivalentSources(block_size=1000),
 100: EquivalentSources(block_size=1000)}

[11]:

# upward continue each line to 10 km
blocked_survey["upward_continued_10km"] = airbornegeo.upward_continue_by_line(
    blocked_survey,
    eqs,
    height=10e3,
)
blocked_survey.head()

Groups:   0%|          | 0/100 [00:00<?, ?it/s]Groups: 100%|██████████| 100/100 [00:00<00:00, 337.91it/s]

[11]:
distance_along_line easting northing height unixtime grav_disturbance_filt line upward_continued_10km
0 477.521406 1.000496e+06 226310.158049 4159.1 1.229507e+09 49.89 1 40.232278
1 1477.462561 1.001484e+06 226460.426302 4160.3 1.229507e+09 50.73 1 40.972657
2 2504.019094 1.002497e+06 226630.494762 4156.0 1.229507e+09 51.09 1 41.713109
3 3518.522095 1.003499e+06 226786.791629 4159.6 1.229507e+09 51.07 1 42.431563
4 4528.427996 1.004498e+06 226936.380167 4165.7 1.229507e+09 51.43 1 43.141623 
[12]:

# plot the upward continued data
airbornegeo.plotly_points(
    blocked_survey,
    color_col="upward_continued_10km",
    robust=True,
    size=2,
)

1.3.1. Sample the satellite gravity grid values into the dataframe#

[13]:

# sample grid along lines
blocked_survey["sampled_grid_values"] = airbornegeo.sample_grid(
    blocked_survey,
    eigen_ds.disturbance,
    coord_names=("easting", "northing"),
)
blocked_survey.head()

grdtrack [WARNING]: Some input points were outside the grid domain(s).

[13]:
distance_along_line easting northing height unixtime grav_disturbance_filt line upward_continued_10km sampled_grid_values
0 477.521406 1.000496e+06 226310.158049 4159.1 1.229507e+09 49.89 1 40.232278 48.675663
1 1477.462561 1.001484e+06 226460.426302 4160.3 1.229507e+09 50.73 1 40.972657 48.836891
2 2504.019094 1.002497e+06 226630.494762 4156.0 1.229507e+09 51.09 1 41.713109 49.005095
3 3518.522095 1.003499e+06 226786.791629 4159.6 1.229507e+09 51.07 1 42.431563 49.172386
4 4528.427996 1.004498e+06 226936.380167 4165.7 1.229507e+09 51.43 1 43.141623 49.342525 
[14]:

# plot the sampled grid values
airbornegeo.plotly_points(
    blocked_survey,
    color_col="sampled_grid_values",
    robust=True,
    size=2,
)

[15]:

# plot the difference
blocked_survey["data_to_grid_diff"] = (
    blocked_survey.upward_continued_10km - blocked_survey.sampled_grid_values
)
airbornegeo.plotly_points(
    blocked_survey,
    color_col="data_to_grid_diff",
    robust=True,
    size=2,
)

1.4. Level the lines to the grid#

[16]:

blocked_survey["levelled_trend_0"] = airbornegeo.level_to_grid(
    blocked_survey,
    degree=0,  # DC shift
    data_column="upward_continued_10km",
    grid_column="sampled_grid_values",
    groupby_column="line",
)
blocked_survey["levelled_trend_1"] = airbornegeo.level_to_grid(
    blocked_survey,
    degree=1,  # DC shift + tilt
    data_column="upward_continued_10km",
    grid_column="sampled_grid_values",
    groupby_column="line",
)
blocked_survey.head()

[16]:
distance_along_line easting northing height unixtime grav_disturbance_filt line upward_continued_10km sampled_grid_values data_to_grid_diff levelled_trend_0 levelled_trend_1
0 477.521406 1.000496e+06 226310.158049 4159.1 1.229507e+09 49.89 1 40.232278 48.675663 -8.443385 42.596650 39.485331
1 1477.462561 1.001484e+06 226460.426302 4160.3 1.229507e+09 50.73 1 40.972657 48.836891 -7.864234 43.337029 40.252951
2 2504.019094 1.002497e+06 226630.494762 4156.0 1.229507e+09 51.09 1 41.713109 49.005095 -7.291987 44.077481 41.021368
3 3518.522095 1.003499e+06 226786.791629 4159.6 1.229507e+09 51.07 1 42.431563 49.172386 -6.740823 44.795935 41.767460
4 4528.427996 1.004498e+06 226936.380167 4165.7 1.229507e+09 51.43 1 43.141623 49.342525 -6.200902 45.505995 42.505032 
[17]:

# plot the levelling correction
blocked_survey["levelling_correction_trend_0"] = (
    blocked_survey.upward_continued_10km - blocked_survey.levelled_trend_0
)
airbornegeo.plotly_points(
    blocked_survey,
    color_col="levelling_correction_trend_0",
    robust=True,
    size=2,
)

[18]:

# plot the levelling correction
blocked_survey["levelling_correction_trend_1"] = (
    blocked_survey.upward_continued_10km - blocked_survey.levelled_trend_1
)
airbornegeo.plotly_points(
    blocked_survey,
    color_col="levelling_correction_trend_1",
    robust=True,
    size=2,
)

[19]:

# look at just 1 of the leveled lines
airbornegeo.plotly_profiles(
    blocked_survey[blocked_survey.line == 2],
    x="distance_along_line",
    y=["sampled_grid_values", "upward_continued_10km", "levelled_trend_0"],
)

[20]:

# look at just 1 of the leveled lines
airbornegeo.plotly_profiles(
    blocked_survey[blocked_survey.line == 2],
    x="distance_along_line",
    y=["sampled_grid_values", "upward_continued_10km", "levelled_trend_1"],
)

1.5. Level the entire survey to the grid#

By not supplying the groupby_column argument, instead of levelling each line in 1D, we can level the entire survey in 2D.

[21]:

blocked_survey["levelled_trend_1"] = airbornegeo.level_to_grid(
    blocked_survey,
    degree=1,
    data_column="upward_continued_10km",
    grid_column="sampled_grid_values",
)
blocked_survey["levelled_trend_2"] = airbornegeo.level_to_grid(
    blocked_survey,
    degree=2,
    data_column="upward_continued_10km",
    grid_column="sampled_grid_values",
)
blocked_survey.head()

[21]:
distance_along_line easting northing height unixtime grav_disturbance_filt line upward_continued_10km sampled_grid_values data_to_grid_diff levelled_trend_0 levelled_trend_1 levelling_correction_trend_0 levelling_correction_trend_1 levelled_trend_2
0 477.521406 1.000496e+06 226310.158049 4159.1 1.229507e+09 49.89 1 40.232278 48.675663 -8.443385 42.596650 48.066543 -2.364372 0.746947 67.674905
1 1477.462561 1.001484e+06 226460.426302 4160.3 1.229507e+09 50.73 1 40.972657 48.836891 -7.864234 43.337029 48.774343 -2.364372 0.719706 68.193819
2 2504.019094 1.002497e+06 226630.494762 4156.0 1.229507e+09 51.09 1 41.713109 49.005095 -7.291987 44.077481 49.481169 -2.364372 0.691741 68.709263
3 3518.522095 1.003499e+06 226786.791629 4159.6 1.229507e+09 51.07 1 42.431563 49.172386 -6.740823 44.795935 50.166525 -2.364372 0.664103 69.204949
4 4528.427996 1.004498e+06 226936.380167 4165.7 1.229507e+09 51.43 1 43.141623 49.342525 -6.200902 45.505995 50.843707 -2.364372 0.636591 69.693441 
[22]:

# plot the levelling correction
blocked_survey["levelling_correction_trend_1"] = (
    blocked_survey.upward_continued_10km - blocked_survey.levelled_trend_1
)
airbornegeo.plotly_points(
    blocked_survey,
    color_col="levelling_correction_trend_1",
    robust=True,
    size=2,
)

[23]:

# plot the levelling correction
blocked_survey["levelling_correction_trend_2"] = (
    blocked_survey.upward_continued_10km - blocked_survey.levelled_trend_2
)
airbornegeo.plotly_points(
    blocked_survey,
    color_col="levelling_correction_trend_2",
    robust=True,
    size=2,
)