Cross-over levelling dependence on line length#

[1]:
%load_ext autoreload
%autoreload 2


import logging

import matplotlib.pyplot as plt
import numpy as np
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"

Load data#

This notebook uses the data and intersections dataframes from the previous notebook Cross-over-errors.

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

data_df = data_df.rename(columns={"distance_along_line": "along_line_distance"})

# define flight lines vs tie lines
data_df["line_type"] = np.where(data_df.line >= 142, 1, 0)

# keep a subset of tie lines and just 1 flight line
line = 77
lines_to_keep = [
    150,
    151,
    152,
    153,
    154,
    155,
    156,
    157,
    158,
    159,
    160,
    161,
    162,
    164,
    165,
    166,
    167,
]
lines_to_keep = np.append(lines_to_keep, line)

data_df = data_df[(data_df.line.isin(lines_to_keep))]

data_df.head()
[2]:
easting northing height line unixtime along_line_distance mag line_type
155245 613796.184341 108259.149034 3793.1 77 1.230970e+09 72.986528 35.54 0
155246 614014.105815 108288.966298 3791.7 77 1.230970e+09 292.943124 35.18 0
155247 614232.757811 108315.269624 3791.9 77 1.230970e+09 513.172378 34.43 0
155248 614414.514181 108337.253486 3793.5 77 1.230970e+09 696.255616 33.60 0
155249 614594.969358 108363.481597 3794.7 77 1.230970e+09 878.614978 32.61 0
[3]:
vmin, vmax = vd.minmax(data_df.mag, min_percentile=5, max_percentile=95)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag",
    s=1,
    cmap="viridis",
    vmin=vmin,
    vmax=vmax,
    title="Unlevelled data",
)
ax.set_aspect("equal")
_images/levelling_06_crossover_levelling_dependence_on_line_length_4_0.png
[ ]:
# Create survey with structural bindings
survey = airbornegeo.Survey(
    data_df,
    line_column="line",
    line_type_column="line_type",
    distance_column="along_line_distance",
)

# Calculate theoretical intersection points
survey.create_intersection_table(method="groups")

# Interpolate mag data at intersections
survey.interpolate_intersections(
    to_interp="mag",
    window_width=500,
    method="cubic",
    extrapolate=False,
)

# Calculate cross-over errors
survey.calculate_crossover_errors(data_col="mag")

# Capture baseline intersections for later resets (forked branching pattern)
inters = survey.intersections

# Update data_df reference to match survey
data_df = survey.data
inters.head()
[ ]:
survey.plot_line_and_crosses(
    y=["mag"], x="along_line_distance", line=line, y_axes=[1], plot_inters=True
)
[6]:
line_df = data_df[data_df.line == line]
ints = inters[(inters.line1 == line) | (inters.line2 == line)]
crossover_error_col = "crossover_error_0"

for ind, row in line_df[line_df.is_intersection].iterrows():
    # search intersections for mistie values
    crossover_error_row = ints[
        ((ints.line1 == line) & (ints.line2 == row.intersecting_line))
        | ((ints.line1 == row.intersecting_line) & (ints.line2 == line))
    ]
    # add misties to line dataframe
    line_df.loc[ind, crossover_error_col] = crossover_error_row[
        crossover_error_col
    ].to_numpy()


# line_df = airbornegeo.interpolate_missing_pointwise(
#     line_df,
#     to_interp=crossover_error_col,
#     interp_on="distance_along_line",
#     method="linear",
#     extrapolate=True,
# )

Level with fitting a trend to the cross-overs#

[ ]:
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend0",
    degree=0,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend0"] = data_df.mag - data_df.mag_levelled_trend0
[ ]:
# Reset intersections to baseline before second branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend1",
    degree=1,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend1"] = data_df.mag - data_df.mag_levelled_trend1
[ ]:
# Reset intersections to baseline before third branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend2",
    degree=2,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend2"] = data_df.mag - data_df.mag_levelled_trend2

fig, axs = plt.subplots(1, 1, figsize=(10, 5))

ax = line_df[line_df.is_intersection].plot.scatter(
    "along_line_distance",
    crossover_error_col,
    marker="^",
    s=30,
    c="r",
    ax=axs,
    label="mistie values",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend0",
    marker="^",
    s=1,
    c="g",
    ax=axs,
    label="Levelling correction trend 0",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend1",
    marker="^",
    s=1,
    c="purple",
    ax=axs,
    label="Levelling correction trend 1",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend2",
    marker="^",
    s=1,
    c="orange",
    ax=axs,
    label="Levelling correction trend 2",
)

# zoom in on the end of the line
max_dist_to_shorten_line = 400e3
ax.set_xlim(
    max_dist_to_shorten_line, data_df[data_df.line == line].along_line_distance.max()
)
ax.set_ylim(-20, 140)
plt.show()

Wavelength-based levelling#

In the above cells, we calculate cross-over errors, and fit trends of specified orders to these values, yielding a levelling correction. For trend orders greater than 1, these levelling corrections are dependent on line length or the position of the cross-overs, in this case, all of them being at the end of the line.

An alternative method is to interpolate these cross-over errors along the entire line and low-pass filter the results. This means we can choose the wavelength of our levelling correction, and our levelling corrections are independent of line length.

[ ]:
# Reset intersections to baseline before fourth branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_filter",
    filter_type="g100000",
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_100km_filter"] = data_df.mag - data_df.mag_levelled_filter
[ ]:
# Reset intersections to baseline before fifth branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_filter",
    filter_type="g500000",
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_500km_filter"] = data_df.mag - data_df.mag_levelled_filter

fig, axs = plt.subplots(1, 1, figsize=(10, 5))

ax = line_df[line_df.is_intersection].plot.scatter(
    "along_line_distance",
    crossover_error_col,
    marker="^",
    s=30,
    c="r",
    ax=axs,
    label="mistie values",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_100km_filter",
    marker="^",
    s=1,
    c="g",
    ax=axs,
    label="Levelling correction 100 km filter",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_500km_filter",
    marker="^",
    s=1,
    c="purple",
    ax=axs,
    label="Levelling correction 500 km filter",
)

# zoom in on the end of the line
ax.set_xlim(
    max_dist_to_shorten_line, data_df[data_df.line == line].along_line_distance.max()
)
ax.set_ylim(-20, 140)
plt.show()

Repeat with a shorter version of the line#

[12]:
# shorten the line of interest
to_drop = data_df[
    (data_df.line == line) & (data_df.along_line_distance < max_dist_to_shorten_line)
]
data_df = data_df.drop(to_drop.index)
[13]:
vmin, vmax = vd.minmax(data_df.mag, min_percentile=5, max_percentile=95)
ax = data_df[::10].plot.scatter(
    "easting",
    "northing",
    c="mag",
    s=1,
    cmap="viridis",
    vmin=vmin,
    vmax=vmax,
    title="Unlevelled data",
)
ax.set_aspect("equal")
_images/levelling_06_crossover_levelling_dependence_on_line_length_17_0.png
[ ]:
# Create survey with structural bindings (on shortened dataset)
survey = airbornegeo.Survey(
    data_df,
    line_column="line",
    line_type_column="line_type",
    distance_column="along_line_distance",
)

# Calculate theoretical intersection points
survey.create_intersection_table(method="groups")

# Interpolate mag data at intersections
survey.interpolate_intersections(
    to_interp="mag",
    window_width=500,
    method="cubic",
    extrapolate=False,
)

# Calculate cross-over errors
survey.calculate_crossover_errors(data_col="mag")

# Capture baseline intersections for this section (forked branching pattern)
inters = survey.intersections

# Update data_df reference to match survey
data_df = survey.data
inters.head()
[ ]:
survey.plot_line_and_crosses(
    y=["mag"], x="along_line_distance", line=line, y_axes=[1], plot_inters=True
)
[16]:
line_df = data_df[data_df.line == line]
ints = inters[(inters.line1 == line) | (inters.line2 == line)]
crossover_error_col = "crossover_error_0"

for ind, row in line_df[line_df.is_intersection].iterrows():
    # search intersections for mistie values
    crossover_error_row = ints[
        ((ints.line1 == line) & (ints.line2 == row.intersecting_line))
        | ((ints.line1 == row.intersecting_line) & (ints.line2 == line))
    ]
    # add misties to line dataframe
    line_df.loc[ind, crossover_error_col] = crossover_error_row[
        crossover_error_col
    ].to_numpy()


# line_df = airbornegeo.interpolate_missing_pointwise(
#     line_df,
#     to_interp=crossover_error_col,
#     interp_on="distance_along_line",
#     method="linear",
#     extrapolate=True,
# )

Level with fitting a trend to the cross-overs#

[ ]:
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend0",
    degree=0,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend0"] = data_df.mag - data_df.mag_levelled_trend0
[ ]:
# Reset intersections to baseline before second branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend1",
    degree=1,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend1"] = data_df.mag - data_df.mag_levelled_trend1
[ ]:
# Reset intersections to baseline before third branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_trend2",
    degree=2,
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_trend2"] = data_df.mag - data_df.mag_levelled_trend2

fig, axs = plt.subplots(1, 1, figsize=(10, 5))

ax = line_df[line_df.is_intersection].plot.scatter(
    "along_line_distance",
    crossover_error_col,
    marker="^",
    s=30,
    c="r",
    ax=axs,
    label="mistie values",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend0",
    marker="^",
    s=1,
    c="g",
    ax=axs,
    label="Levelling correction trend 0",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend1",
    marker="^",
    s=1,
    c="purple",
    ax=axs,
    label="Levelling correction trend 1",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_trend2",
    marker="^",
    s=1,
    c="orange",
    ax=axs,
    label="Levelling correction trend 2",
)

ax.set_xlim(
    max_dist_to_shorten_line, data_df[data_df.line == line].along_line_distance.max()
)
ax.set_ylim(-20, 140)

plt.show()

Wavelength-based levelling#

In the above cells, we calculate cross-over errors, and fit trends of specified orders to these values, yielding a levelling correction. For trend orders greater than 1, these levelling corrections are dependent on line length or the position of the cross-overs, in this case, all of them being at the end of the line.

An alternative method is to interpolate these cross-over errors along the entire line and low-pass filter the results. This means we can choose the wavelength of our levelling correction, and our levelling corrections are independent of line length.

[ ]:
# Reset intersections to baseline before fourth branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_filter",
    filter_type="g100000",
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_100km_filter"] = data_df.mag - data_df.mag_levelled_filter
[ ]:
# Reset intersections to baseline before fifth branch
survey.intersections = inters
survey.crossover_pair_levelling(
    lines_to_level=[line],
    data_col="mag",
    levelled_col="mag_levelled_filter",
    filter_type="g500000",
)

# Update data_df reference
data_df = survey.data

# Calculate levelling correction
data_df["levelling_correction_500km_filter"] = data_df.mag - data_df.mag_levelled_filter

fig, axs = plt.subplots(1, 1, figsize=(10, 5))

ax = line_df.plot.scatter(
    "along_line_distance",
    crossover_error_col,
    s=1,
    ax=axs,
    label="interpolated mistie values",
)
ax = line_df[line_df.is_intersection].plot.scatter(
    "along_line_distance",
    crossover_error_col,
    marker="^",
    s=30,
    c="r",
    ax=axs,
    label="mistie values",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_100km_filter",
    marker="^",
    s=1,
    c="g",
    ax=axs,
    label="Levelling correction 100 km filter",
)
ax = data_df[data_df.line == line].plot.scatter(
    "along_line_distance",
    "levelling_correction_500km_filter",
    marker="^",
    s=1,
    c="purple",
    ax=axs,
    label="Levelling correction 500 km filter",
)

ax.set_xlim(
    max_dist_to_shorten_line, data_df[data_df.line == line].along_line_distance.max()
)
ax.set_ylim(-20, 140)
plt.show()
[ ]: