3. Cross-over errors

%load_ext autoreload
%autoreload 2


import logging

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

import airbornegeo

# setup logging to get some additional info from the airbornegeo functions
logging.getLogger("airbornegeo").setLevel("INFO")
logging.basicConfig()
pio.renderers.default = "notebook"

/home/sungw937/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

3.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)
]

# 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

data_df[data_df.tie].line.sort_values().unique()

array([143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,
       156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 177,
       178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190,
       191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203,
       204, 205, 206])

data_df[~data_df.tie].line.sort_values().unique()

array([  1,   3,   5,   7,   9,  11,  13,  15,  17,  19,  21,  23,  25,
        27,  29,  31,  33,  35,  37,  39,  41,  43,  45,  47,  49,  51,
        53,  55,  57,  59,  61,  63,  65,  67,  69,  71,  73,  75,  77,
        79,  81,  83,  85,  87,  89,  91,  93,  95,  97,  99, 101, 103,
       105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 125, 127, 129,
       131])

airbornegeo.plotly_points(
    data_df,
    color_col="line",
    hover_cols=["distance_along_line"],
)

3.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)
inters

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

	line	tie	geometry	max_dist	easting	northing
0	1	143	POINT (1158153 254083)	16.441255	1158153.0	254083.0
1	1	144	POINT (1190820 259865)	24.124887	1190820.0	259865.0
2	1	145	POINT (1223524 265689)	17.269650	1223524.0	265689.0
3	1	146	POINT (1256193 271497)	31.608179	1256193.0	271497.0
4	1	147	POINT (1288901 277249)	50.174377	1288901.0	277249.0
...	...	...	...	...	...	...
657	125	199	POINT (1318165 302398)	33.689493	1318165.0	302398.0
658	125	200	POINT (1350929 308239)	36.343660	1350929.0	308239.0
659	127	148	POINT (1319894 292699)	22.093103	1319894.0	292699.0
660	127	199	POINT (1319922 292702)	23.297236	1319922.0	292702.0
661	127	200	POINT (1352642 298507)	36.743873	1352642.0	298507.0

662 rows × 6 columns

3.3. Add intersections as rows to the dataframe

We can now add a row to each line of the dataframe which represents the intersection point. We don’t know the data values for the intersection point, so we will also need to interpolate the data from the line data on either side of the intersection. We can specify the interpolation type (‘linear’, ‘nearest’, ‘quadratic’, ‘cubic’ etc.) and a window width of data to use for the interpolation in meters. We also specify which column we want to be interpolated (the magnetic data).

If the window width is not wide enough to capture data on either side, it will be doubled twice to try and capture enough data. If it still too narrow to capture data on either side, the intersection point will either be dropped, or if extrapolate is set to true, the values will be filled with extrapolation, instead of interpolation. This extrapolation is also necessary if the intersections lie beyond the end of a line, which occurs if you used the buffer_dist argument for create_intersection_table.

 - : data
 x : intersection point
[ ]: data windows for interpolation
- - - - - - -          x         - - - - - -
                    [     ] # first window
                 [           ] # doubled once
            [                     ] # doubled twice

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

Segments: 100%|██████████| 121/121 [00:08<00:00, 14.46it/s]

inters.head()

	line	tie	geometry	max_dist	easting	northing	dist_along_flight_line	dist_along_flight_tie	flight_interpolation_type	tie_interpolation_type
0	1	143	POINT (1158153 254083)	16.441255	1158153.0	254083.0	545768.072695	138071.419338	interpolated	interpolated
1	1	144	POINT (1190820 259865)	24.124887	1190820.0	259865.0	578952.394781	131585.125284	interpolated	interpolated
2	1	145	POINT (1223524 265689)	17.269650	1223524.0	265689.0	612181.972720	147253.496664	interpolated	interpolated
3	1	146	POINT (1256193 271497)	31.608179	1256193.0	271497.0	645372.035087	136211.666619	interpolated	interpolated
4	1	147	POINT (1288901 277249)	50.174377	1288901.0	277249.0	678598.170740	155115.224990	interpolated	interpolated

data_df[data_df.is_intersection].head()

	easting	northing	height	line	unixtime	distance_along_line	mag	tie	geometry	is_intersection	intersecting_line	mag_interpolation_type
7453	637640.0	161826.0	NaN	1	NaN	16811.103180	-53.846042	NaN	POINT (637640 161826)	True	150.0	interpolated
7454	670007.0	167673.0	NaN	1	NaN	49821.820318	-80.387990	NaN	POINT (670007 167673)	True	151.0	interpolated
7455	702402.0	173414.0	NaN	1	NaN	82726.336899	-44.550155	NaN	POINT (702402 173414)	True	152.0	interpolated
7456	702430.0	173419.0	NaN	1	NaN	82767.019948	-44.908585	NaN	POINT (702430 173419)	True	153.0	interpolated
7457	734809.0	179103.0	NaN	1	NaN	115638.011785	-120.195968	NaN	POINT (734809 179103)	True	154.0	interpolated

# see which lines dont have intersections
airbornegeo.lines_without_intersections(data_df, 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)]

3.4. Inspect the intersections

The below two cells allow you to look at a line in profile view and see it’s intersections. We pass the column ‘mag’ to be plotted on the y axis, and distance along line on the x axis. The small dots, which form lines since they are close together, show the magnetic data of the line. The diamonds show the intersection points. The diamonds are plotted with the y axis values of the intersecting line. If the diamond if vertically offset from the points, that means the intersecting line’s magnetic value at the intersection is quite different. This shows the cross-over errors. For example, in the below plot, the lines 85 and 77 have quite large cross-over erros with the plotted line, 161.

# # if you want to individually inspect each line and its intersections
# airbornegeo.inspect_intersections(
#     data_df,
#     plot_variable=["mag"],
#     # plot_all=True,
# )

# if you want to just look at 1 line
airbornegeo.plot_line_and_crosses(
    data_df,
    line=161,
    x="distance_along_line",
    y=["mag"],
    y_axes=[1],
    plot_inters=True,
)

3.5. 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,
)

inters.head()

	line	tie	geometry	max_dist	easting	northing	dist_along_flight_line	dist_along_flight_tie	flight_interpolation_type	tie_interpolation_type	mistie_0
0	1	143	POINT (1158153 254083)	16.441255	1158153.0	254083.0	545768.072695	138071.419338	interpolated	interpolated	100.458661
1	1	144	POINT (1190820 259865)	24.124887	1190820.0	259865.0	578952.394781	131585.125284	interpolated	interpolated	32.151483
2	1	145	POINT (1223524 265689)	17.269650	1223524.0	265689.0	612181.972720	147253.496664	interpolated	interpolated	72.936553
3	1	146	POINT (1256193 271497)	31.608179	1256193.0	271497.0	645372.035087	136211.666619	interpolated	interpolated	84.125118
4	1	147	POINT (1288901 277249)	50.174377	1288901.0	277249.0	678598.170740	155115.224990	interpolated	interpolated	-7.731396