Cross-over errors

3. Cross-over errors#

[1]:

%load_ext autoreload
%autoreload 2


import logging

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

import airbornegeo

logging.getLogger("airbornegeo").setLevel("INFO")
logging.basicConfig()
pio.renderers.default = "notebook"

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.

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

# 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()

[2]:
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 
[3]:

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

[3]:

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

[4]:

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

[4]:

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

[5]:

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

[6]:

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

3.2. Find intersections#

[7]:

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

[8]:

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

INFO:airbornegeo:found 662 intersections

[8]:
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/cross_over_errors_10_5.png

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 columns we want to be interpolated.

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
[ ]: windows
- - - - - - -          x         - - - - - -
                    [     ] # first window
                 [           ] # doubled once
            [                     ] # doubled twice

[9]:

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

[10]:

inters.head()

[10]:
line tie geometry max_dist easting northing dist_along_flight_line dist_along_flight_tie flight_interpolation_type tie_interpolation_type flight_height tie_height
0 1 143 POINT (1158153 254083) 16.441255 1158153.0 254083.0 545768.072695 138071.419338 interpolated interpolated 4168.038983 4008.397691
1 1 144 POINT (1190820 259865) 24.124887 1190820.0 259865.0 578952.394781 131585.125284 interpolated interpolated 4174.500813 4010.840823
2 1 145 POINT (1223524 265689) 17.269650 1223524.0 265689.0 612181.972720 147253.496664 interpolated interpolated 3544.149925 3995.961216
3 1 146 POINT (1256193 271497) 31.608179 1256193.0 271497.0 645372.035087 136211.666619 interpolated interpolated 3564.579994 3578.174277
4 1 147 POINT (1288901 277249) 50.174377 1288901.0 277249.0 678598.170740 155115.224990 interpolated interpolated 3536.160719 3545.947314 
[11]:

data_df[data_df.is_intersection].head()

[11]:
easting northing height line unixtime distance_along_line mag tie geometry is_intersection intersecting_line interpolation_type
7455 637640.0 161826.0 4122.332685 1 NaN 16811.103180 -53.846042 NaN POINT (637640 161826) True 150.0 interpolated
7456 670007.0 167673.0 4132.634362 1 NaN 49821.820318 -80.387990 NaN POINT (670007 167673) True 151.0 interpolated
7457 702402.0 173414.0 4156.031337 1 NaN 82726.336899 -44.550155 NaN POINT (702402 173414) True 152.0 interpolated
7458 702430.0 173419.0 4156.517464 1 NaN 82767.019948 -44.908585 NaN POINT (702430 173419) True 153.0 interpolated
7459 734809.0 179103.0 4121.060700 1 NaN 115638.011785 -120.195968 NaN POINT (734809 179103) True 154.0 interpolated 
[12]:

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

[12]:

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

[13]:

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

[14]:

# 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

[15]:

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

inters.head()

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