Tie line levelling

5. Tie line levelling#

[1]:

# %load_ext autoreload
# %autoreload 2


import logging

import cmocean
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import plotly.io as pio
import verde as vd

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

5.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.

[2]:

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",
    ]
]

# to run faster, 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

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

data_df.head()

[2]:
easting northing height line unixtime distance_along_line mag tie geometry
0 621099.093988 159056.748193 4110.45 1 1.229500e+09 27.181565 -33.385 False POINT (621099.094 159056.748)
1 621206.517953 159071.301512 4114.50 1 1.229500e+09 135.587530 -35.120 False POINT (621206.518 159071.302)
2 621313.794221 159085.113599 4117.90 1 1.229500e+09 243.749367 -36.875 False POINT (621313.794 159085.114)
3 621420.722903 159099.184618 4120.80 1 1.229500e+09 351.600337 -38.585 False POINT (621420.723 159099.185)
4 621527.158177 159114.303863 4123.15 1 1.229500e+09 459.104857 -40.205 False POINT (621527.158 159114.304) 
[3]:

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

5.2. Find intersections#

[4]:

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

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

[4]:
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

_images/tie_line_levelling_6_3.png

5.3. Add intersections as rows to the dataframe#

[5]:

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

Segments: 100%|██████████| 121/121 [00:09<00:00, 13.43it/s]

[6]:

# see which lines dont have intersections
airbornegeo.lines_without_intersections(data_df, inters)

[6]:

[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)]

5.4. Calculate initial cross-over errors#

[7]:

inters = airbornegeo.calculate_crossover_errors(
    data_df,
    inters,
    data_col="mag",
    plot_map=True,
)
_images/tie_line_levelling_11_1.png 
[8]:

inters.head()

[8]:
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

5.5. Level lines to ties#

Now we will level only the lines, holding the ties constant. We will just use a trend degree of 0, which allows only a DC shift of the lines. We will save the levelled data to a new column mag_levelled_lines_trend0. If you don’t want to keep track on new columns, you can just use the same name as the data column.

We assign the returned inters dataframe to a new variable name so later on in the notebook we can access the original again.

[9]:

data_df, inters_lines_to_ties = airbornegeo.line_levelling(
    data_df,
    inters,
    lines_to_level=data_df[data_df.tie == False].line.unique(),
    data_col="mag",
    levelled_col="mag_levelled_lines_to_ties_trend0",
    degree=0,
)
inters_lines_to_ties.head()

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

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

[11]:

airbornegeo.plot_line_and_crosses(
    data_df,
    line=5,
    x="distance_along_line",
    y=["mag", "mag_levelled_lines_to_ties_trend0"],
    y_axes=[1, 1],
    plot_inters=[True, False],
)

[12]:

inters_lines_to_ties = airbornegeo.calculate_crossover_errors(
    data_df,
    inters_lines_to_ties,
    data_col="mag_levelled_lines_to_ties_trend0",
    plot_map=True,
)
_images/tie_line_levelling_17_1.png 
[13]:

airbornegeo.plot_levelling_convergence(inters_lines_to_ties)
_images/tie_line_levelling_18_0.png

From the above profile, we can see line 77 has been shifted down to try and minimize the cross-over errors. The map, histogram and levelling convergence figures show we have reduced the cross-over errors, bringing the RMSE from ~43nT to ~24nT.

[14]:

fig, axs = plt.subplots(1, 3, figsize=(15, 40))

max_abs = vd.maxabs(data_df.mag, percentile=95)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag",
    s=1,
    ax=axs[0],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Unlevelled",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

data_df["levelling_correction"] = (
    data_df.mag - data_df.mag_levelled_lines_to_ties_trend0
)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="levelling_correction",
    s=1,
    ax=axs[1],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Levelling correction",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag_levelled_lines_to_ties_trend0",
    s=1,
    ax=axs[2],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Lines levelled to ties, trend 0",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

plt.tight_layout()
plt.show()
_images/tie_line_levelling_20_0.png

5.6. Level ties to lines#

We can also level the tie lines to the flight lines. We do this be supplying the lines_to_level parameter with the names of the tie lines. We will give this levelled data a new name mag_levelled_ties_trend0

[15]:

data_df, inters_ties_to_lines = airbornegeo.line_levelling(
    data_df,
    inters,
    lines_to_level=data_df[data_df.tie == True].line.unique(),
    data_col="mag",
    levelled_col="mag_levelled_ties_to_lines_trend0",
    degree=0,
)
inters_ties_to_lines.head()

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

inters_ties_to_lines = airbornegeo.calculate_crossover_errors(
    data_df,
    inters_ties_to_lines,
    data_col="mag_levelled_ties_to_lines_trend0",
    plot_map=True,
)
_images/tie_line_levelling_23_1.png 
[17]:

airbornegeo.plot_levelling_convergence(inters_ties_to_lines)
_images/tie_line_levelling_24_0.png 
[18]:

fig, axs = plt.subplots(1, 3, figsize=(15, 40))

max_abs = vd.maxabs(data_df.mag, percentile=95)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag",
    s=1,
    ax=axs[0],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Unlevelled",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

data_df["levelling_correction"] = (
    data_df.mag - data_df.mag_levelled_ties_to_lines_trend0
)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="levelling_correction",
    s=1,
    ax=axs[1],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Levelling correction",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag_levelled_ties_to_lines_trend0",
    s=1,
    ax=axs[2],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Ties levelled to lines, trend 0",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

plt.tight_layout()
plt.show()
_images/tie_line_levelling_25_0.png

We can repeat these steps, increasing the trend order if we need more levelling, and alternating between levelling lines or ties. In some case, we want to iteratively repeat the levelling. This is shown below.

5.7. Iteratively level lines to ties#

This will perform up to 5 iterations of levelling the lines to the ties. Once the cross-over errors stop changing, the iterations will stop. Most of the time it will stop after 1 iteration, but as you’ll see in further notebooks, ones we performed weighted levelling, more iterations can be useful.

[19]:

data_df, inters_iterative = airbornegeo.iterative_line_levelling(
    data_df,
    inters,
    lines_to_level=data_df[data_df.tie == False].line.unique(),
    data_col="mag",
    levelled_col="mag_levelled_lines_to_ties_trend0_iterative",
    degree=0,
    iterations=5,
)
inters_iterative.head()

Iteration:  20%|██        | 1/5 [00:18<01:14, 18.60s/it]

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

airbornegeo.plot_levelling_convergence(inters_iterative)
_images/tie_line_levelling_29_0.png

Instead of just iteratively repeating levelling lines to ties, for each iteration we can alternative between levelling lines to ties, and ties to lines, as shown below.

5.8. Alternativing iterative levelling#

This will perform up to 5 iterations where each iteration first levels the lines to the ties, then the ties to the lines. Since each iteration performs two instances of leveling, there will be up to 10 new misties columns in the intersection table. If the misties values don’t decrease, or begin increasing, the iterations will stop. We will also increase the degree trend order to 1, to allow the lines to be shifted as well as tilted.

[21]:

data_df, inters_alternating_iterative = (
    airbornegeo.alternating_iterative_line_levelling(
        data_df,
        inters,
        data_col="mag",
        levelled_col="mag_levelled_alternating_trend1",
        degree=1,
        iterations=5,
    )
)
inters_alternating_iterative.head()

Iteration: 100%|██████████| 5/5 [01:19<00:00, 15.97s/it]

[21]:
line tie geometry max_dist easting northing dist_along_flight_line dist_along_flight_tie flight_interpolation_type tie_interpolation_type ... mistie_1 mistie_2 mistie_3 mistie_4 mistie_5 mistie_6 mistie_7 mistie_8 mistie_9 mistie_10
0 1 143 POINT (1158153 254083) 16.441255 1158153.0 254083.0 545768.072695 138071.419338 interpolated interpolated ... 85.390637 53.830704 50.428968 47.406472 45.290888 43.475684 42.384815 41.138719 40.572418 39.700654
1 1 144 POINT (1190820 259865) 24.124887 1190820.0 259865.0 578952.394781 131585.125284 interpolated interpolated ... 14.857246 4.930539 1.394638 -1.040580 -3.329380 -4.819601 -6.004392 -7.036234 -7.651205 -8.375117
2 1 145 POINT (1223524 265689) 17.269650 1223524.0 265689.0 612181.972720 147253.496664 interpolated interpolated ... 53.413067 37.042316 33.372066 31.524115 29.061863 27.896855 26.618016 25.800627 25.136921 24.561012
3 1 146 POINT (1256193 271497) 31.608179 1256193.0 271497.0 645372.035087 136211.666619 interpolated interpolated ... 62.375033 48.273736 44.469298 43.210468 40.574970 39.736369 38.363591 37.761503 37.049120 36.621833
4 1 147 POINT (1288901 277249) 50.174377 1288901.0 277249.0 678598.170740 155115.224990 interpolated interpolated ... -31.710499 -10.507246 -14.446019 -15.118106 -17.927039 -18.440708 -19.907525 -20.295314 -21.056429 -21.335793

5 rows × 21 columns

[22]:

airbornegeo.plot_levelling_convergence(inters_alternating_iterative)
_images/tie_line_levelling_32_0.png 
[23]:

fig, axs = plt.subplots(1, 3, figsize=(15, 40))

max_abs = vd.maxabs(data_df.mag, percentile=95)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag",
    s=1,
    ax=axs[0],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Unlevelled",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

data_df["levelling_correction"] = data_df.mag - data_df.mag_levelled_alternating_trend1
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="levelling_correction",
    s=1,
    ax=axs[1],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Levelling correction",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag_levelled_alternating_trend1",
    s=1,
    ax=axs[2],
    cmap=cmocean.cm.balance,
    vmin=-max_abs,
    vmax=max_abs,
    colorbar=False,
    title="Alternativing iterative levelling, trend 1",
)
ax.set_yticks([])
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.1)

plt.tight_layout()
plt.show()
_images/tie_line_levelling_33_0.png