4. Using equivalent sources for cross-over errors

For some data types, such as gravity and magnetics, the value of the data depends on the the elevation it was collected at. For an intersection, if the two lines were flown at different altitudes, the intersection point is only the intersection in 2D space, not in 3D. To make it a true 3D intersection, we can perform upward continuation on the data so that the values of both lines represent the values that would have been observed if the flights were flown at the same altitude. We don’t actually need to upward continue the entire lines, just the single points which make up the intersection.

%load_ext autoreload
%autoreload 2


import logging

import geopandas as gpd
import numpy as np
import pandas as pd
import plotly.io as pio

import airbornegeo

logging.getLogger("airbornegeo").setLevel("INFO")
logging.basicConfig()
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

4.1. Load 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.

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

# for testing limit number of lines
data_df = data_df[~data_df.line.between(133, 142)]
data_df = data_df[~data_df.line.between(168, 176)]
data_df = data_df[
    (data_df.line.isin(data_df.line.unique()[::2])) | (data_df.line >= 143)
]

# drop rows with any nans
data_df = data_df.dropna(how="any")

# define flight lines vs tie lines with column 'tie' which is True or False
data_df["tie"] = False
data_df.loc[data_df.line >= 142, "tie"] = True

data_df.head()

	easting	northing	height	line	unixtime	distance_along_line	mag	tie
0	621099.093988	159056.748193	4110.45	1	1.229500e+09	27.181565	-33.385	False
1	621206.517953	159071.301512	4114.50	1	1.229500e+09	135.587530	-35.120	False
2	621313.794221	159085.113599	4117.90	1	1.229500e+09	243.749367	-36.875	False
3	621420.722903	159099.184618	4120.80	1	1.229500e+09	351.600337	-38.585	False
4	621527.158177	159114.303863	4123.15	1	1.229500e+09	459.104857	-40.205	False

4.2. Find intersections

# convert dataframe into geodataframe
data_df = gpd.GeoDataFrame(
    data_df,
    geometry=gpd.points_from_xy(data_df.easting, data_df.northing),
    crs="EPSG:3031",
)

# calculate theoretical intersection points
inters = airbornegeo.create_intersection_table(data_df)

Line/tie combinations: 100%|██████████| 3630/3630 [00:00<00:00, 7052.12it/s]
Potential intersections: 100%|██████████| 670/670 [00:06<00:00, 106.66it/s]
INFO:airbornegeo:found 662 intersections

4.3. Add intersections as rows to the dataframe

This just performs a standard interpolation, we haven’t used equivalent sources yet.

data_df, inters = airbornegeo.interpolate_intersections(
    data_df,
    inters,
    to_interp=["mag", "height"],
    window_width=500,
    method="cubic",
    extrapolate=False,
)

Line 206: 100%|██████████| 121/121 [00:08<00:00, 14.21it/s]

lines_without_inters = airbornegeo.lines_without_intersections(data_df, inters)
lines_without_inters

[np.int64(107),
 np.int64(109),
 np.int64(129),
 np.int64(131),
 np.int64(188),
 np.int64(189),
 np.int64(190),
 np.int64(192),
 np.int64(193),
 np.int64(194),
 np.int64(203)]

# drop lines without intersections
data_df = data_df[~data_df.line.isin(lines_without_inters)]

4.4. Calculate intersection cross-over errors

Now we have interpolated magnetic values for each intersecting line at the intersection points, we can calculate the cross-over errors for each intersection

inters = airbornegeo.calculate_crossover_errors(
    data_df,
    inters,
    data_col="mag",
    plot_map=True,
)

4.5. Inspect flight altitudes

airbornegeo.plotly_points(
    data_df[::50],
    color_col="height",
    hover_cols=["line"],
)

4.6. Inspect flight altitude difference at intersection points

inters["inter_height_diff"] = np.abs(inters["flight_height"] - inters["tie_height"])
airbornegeo.plotly_points(
    inters,
    color_col="inter_height_diff",
    hover_cols=["line", "tie"],
    robust=False,
    size=6,
)

airbornegeo.plot_line_and_crosses(
    data_df,
    line=148,
    x="distance_along_line",
    y=["height", "mag"],
    y_axes=[1, 2],
    plot_inters=True,
)

From the above plot showing line 148, we can see many of the cross lines were flown several hundred meters higher or lower.

4.7. Fit equivalent sources to each line

eqs = airbornegeo.eq_sources_1d(
    data_df,
    data_column="mag",
    depth="default",
    damping=None,
    block_size=500,
    groupby_column="line",
)

Groups: 100%|██████████| 110/110 [01:36<00:00,  1.14it/s]

4.8. Upward continue the intersection points to common altitudes

For each intersection point, we find the highest of the two flight elevations, and upward continue the lower data point to this height. This only changes the individual intersection point, not any of the rest of the line. You can see in the below profile of line 148 that the line data is the same between mag and the new column mag_eqs, but if you look at the intersection points (diamonds), you will see some of the have shifted slightly.

data_df["mag_eqs"] = airbornegeo.update_intersections_with_eq_sources(
    data_df,
    fitted_equivalent_sources=eqs,
    data_column="mag",
)

Segments: 100%|██████████| 110/110 [00:03<00:00, 31.51it/s]

airbornegeo.plot_line_and_crosses(
    data_df,
    line=148,
    y=["height", "mag", "mag_eqs"],
    y_axes=[1, 2, 2],
    plot_inters=[True, True, True],
)

4.9. Recalculate cross-over errors with the updated intersection values

inters = airbornegeo.calculate_crossover_errors(
    data_df,
    inters,
    data_col="mag_eqs",
    plot_map=True,
)

airbornegeo.plot_levelling_convergence(inters)

We can see that by just upward continuing the intersection points, we have slightly lowered the cross-over errors.