2. Eotvos Correction#
The Eotvos correction is a necessary correction to apply for any gravity surveys conducted on a moving platform, such as a ship or airplane. The correction accounts for the relative motion between the vehicle and the Earth’s surface, which generates an additional centrifugal force known as the Eotvos effect.
[1]:
%load_ext autoreload
%autoreload 2
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import plotly.io as pio
import seaborn as sns
import airbornegeo
pio.renderers.default = "notebook"
2.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_survey. It has a pre-compute Eotvos correction, which we will compare our computed values to.
[2]:
data_df_full = pd.read_csv("data/AGAP_gravity_survey.csv")
print(data_df_full.columns)
data_df_full["index"] = data_df_full.index
data_df_full.head()
Index(['Lon', 'Lat', 'Height_WGS1984', 'Date', 'Time', 'ST', 'CC', 'RB',
'XACC', 'LACC', 'Still', 'Base', 'ST_real', 'Beam_vel', 'rec_grav',
'Abs_grav', 'VaccCor', 'EotvosCor', 'LatCor', 'FaCor', 'HaccCor',
'Free_air', 'FAA_filt', 'FAA_clip', 'Level_cor', 'FAA_level',
'Fa_4600m', 'easting', 'northing', 'line_name', 'line', 'unixtime'],
dtype='str')
[2]:
| Lon | Lat | Height_WGS1984 | Date | Time | ST | CC | RB | XACC | LACC | ... | FAA_clip | Level_cor | FAA_level | Fa_4600m | easting | northing | line_name | line | unixtime | index | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 77.252450 | -80.583923 | 4156.1 | 2008-12-17 | 0 days 09:42:48 | 11934.47 | 2.61 | -659.0 | -49.0 | 273.0 | ... | 49.38 | 7.03 | 42.4 | 40.8 | 1.000024e+06 | 226237.330771 | 11_DA500 | 1 | 1.229507e+09 | 0 |
| 1 | 77.252672 | -80.583377 | 4156.0 | 2008-12-17 | 0 days 09:42:49 | 11934.47 | 2.72 | -368.6 | -321.0 | 230.0 | ... | 49.45 | 7.04 | 42.4 | 40.8 | 1.000083e+06 | 226246.631269 | 11_DA500 | 1 | 1.229507e+09 | 1 |
| 2 | 77.252901 | -80.582831 | 4156.1 | 2008-12-17 | 0 days 09:42:50 | 11888.95 | -2.08 | 703.1 | 433.0 | 146.0 | ... | 49.52 | 7.04 | 42.5 | 40.9 | 1.000142e+06 | 226255.809132 | 11_DA500 | 1 | 1.229507e+09 | 2 |
| 3 | 77.253131 | -80.582285 | 4156.4 | 2008-12-17 | 0 days 09:42:51 | 11888.95 | 0.50 | 625.1 | 566.0 | 223.0 | ... | 49.58 | 7.03 | 42.5 | 40.9 | 1.000201e+06 | 226264.969079 | 11_DA500 | 1 | 1.229507e+09 | 3 |
| 4 | 77.253358 | -80.581740 | 4156.6 | 2008-12-17 | 0 days 09:42:52 | 11888.95 | -1.73 | 575.1 | 108.0 | 205.0 | ... | 49.65 | 7.04 | 42.6 | 41.0 | 1.000260e+06 | 226274.156809 | 11_DA500 | 1 | 1.229507e+09 | 4 |
5 rows × 33 columns
[3]:
# turn to geopandas geodataframe
data_df_full = gpd.GeoDataFrame(
data_df_full,
geometry=gpd.points_from_xy(x=data_df_full.easting, y=data_df_full.northing),
crs="EPSG:3031",
)
[4]:
# get only the raw columns
# we will perform the corrections ourselves and compare to their values
data_df = data_df_full[
[
"EotvosCor",
"Lon",
"Lat",
"Height_WGS1984",
"easting",
"northing",
"unixtime",
"line",
"geometry",
"index",
]
]
data_df.head()
[4]:
| EotvosCor | Lon | Lat | Height_WGS1984 | easting | northing | unixtime | line | geometry | index | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 68.15 | 77.252450 | -80.583923 | 4156.1 | 1.000024e+06 | 226237.330771 | 1.229507e+09 | 1 | POINT (1000023.79 226237.331) | 0 |
| 1 | 68.37 | 77.252672 | -80.583377 | 4156.0 | 1.000083e+06 | 226246.631269 | 1.229507e+09 | 1 | POINT (1000082.905 226246.631) | 1 |
| 2 | 68.37 | 77.252901 | -80.582831 | 4156.1 | 1.000142e+06 | 226255.809132 | 1.229507e+09 | 1 | POINT (1000142.048 226255.809) | 2 |
| 3 | 68.16 | 77.253131 | -80.582285 | 4156.4 | 1.000201e+06 | 226264.969079 | 1.229507e+09 | 1 | POINT (1000201.195 226264.969) | 3 |
| 4 | 68.12 | 77.253358 | -80.581740 | 4156.6 | 1.000260e+06 | 226274.156809 | 1.229507e+09 | 1 | POINT (1000260.224 226274.157) | 4 |
2.2. Aircraft velocity#
One component of the Eotvos correction calculation is the velocity of the aircraft. Below we calculated the relative distance between successive rows in the dataframe. This assumes the dataframe is sorted by line, and then by time. From the relative distance, we can then calculate the ground speed of the aircraft.
[5]:
data_df["relative_distance"] = airbornegeo.relative_distance(data_df)
data_df["ground_speed"] = airbornegeo.ground_speed(data_df)
airbornegeo.plotly_points(
data_df[::50],
color_col="ground_speed",
)
data_df.head()
[5]:
| EotvosCor | Lon | Lat | Height_WGS1984 | easting | northing | unixtime | line | geometry | index | relative_distance | ground_speed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 68.15 | 77.252450 | -80.583923 | 4156.1 | 1.000024e+06 | 226237.330771 | 1.229507e+09 | 1 | POINT (1000023.79 226237.331) | 0 | 0.000000 | NaN |
| 1 | 68.37 | 77.252672 | -80.583377 | 4156.0 | 1.000083e+06 | 226246.631269 | 1.229507e+09 | 1 | POINT (1000082.905 226246.631) | 1 | 59.842447 | 59.842447 |
| 2 | 68.37 | 77.252901 | -80.582831 | 4156.1 | 1.000142e+06 | 226255.809132 | 1.229507e+09 | 1 | POINT (1000142.048 226255.809) | 2 | 59.850955 | 59.850955 |
| 3 | 68.16 | 77.253131 | -80.582285 | 4156.4 | 1.000201e+06 | 226264.969079 | 1.229507e+09 | 1 | POINT (1000201.195 226264.969) | 3 | 59.852244 | 59.852244 |
| 4 | 68.12 | 77.253358 | -80.581740 | 4156.6 | 1.000260e+06 | 226274.156809 | 1.229507e+09 | 1 | POINT (1000260.224 226274.157) | 4 | 59.739529 | 59.739529 |
2.3. Aircraft track#
Another component of the Eotvos correction calculation is the track of the aircraft. This is similar to the aircraft’s heading. The heading is the angle the nose of the plane points, which due to wind, can differ from the angle the path of the plane makes. This is in degrees clockwise from geographic north.
By setting ellipsoid=True, we calculate the track using the WGS84 ellipsoid, instead of using a simplified spheric model of the Earth. It is more accurate, but a bit slower to compute.
[6]:
data_df["track"] = airbornegeo.track(
data_df,
latitude_column="Lat",
longitude_column="Lon",
groupby_column="line",
ellipsoid=False, # this is faster
)
airbornegeo.plotly_points(data_df[::50], color_col="track", hover_cols=["unixtime"])
data_df.head()
[6]:
| EotvosCor | Lon | Lat | Height_WGS1984 | easting | northing | unixtime | line | geometry | index | relative_distance | ground_speed | track | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 68.15 | 77.252450 | -80.583923 | 4156.1 | 1.000024e+06 | 226237.330771 | 1.229507e+09 | 1 | POINT (1000023.79 226237.331) | 0 | 0.000000 | NaN | NaN |
| 1 | 68.37 | 77.252672 | -80.583377 | 4156.0 | 1.000083e+06 | 226246.631269 | 1.229507e+09 | 1 | POINT (1000082.905 226246.631) | 1 | 59.842447 | 59.842447 | 3.805916 |
| 2 | 68.37 | 77.252901 | -80.582831 | 4156.1 | 1.000142e+06 | 226255.809132 | 1.229507e+09 | 1 | POINT (1000142.048 226255.809) | 2 | 59.850955 | 59.850955 | 3.925778 |
| 3 | 68.16 | 77.253131 | -80.582285 | 4156.4 | 1.000201e+06 | 226264.969079 | 1.229507e+09 | 1 | POINT (1000201.195 226264.969) | 3 | 59.852244 | 59.852244 | 3.943093 |
| 4 | 68.12 | 77.253358 | -80.581740 | 4156.6 | 1.000260e+06 | 226274.156809 | 1.229507e+09 | 1 | POINT (1000260.224 226274.157) | 4 | 59.739529 | 59.739529 | 3.899163 |
2.4. Plot the Eotvos correction included in the dataset#
Below we show the already computed Eotvos correction from the AGAP dataset.
[7]:
airbornegeo.plotly_points(
data_df_full[::50],
color_col="EotvosCor",
)
2.5. Eotvos correction following Glicken 1962#
We offer several methods of computing the Eotvos correction. The most simple is the Glicken 1962 formular, which depends on only the latitudes, the track angle, and the ground speed.
[8]:
data_df["eotvos_correction"] = airbornegeo.eotvos_correction_glicken(
data_df.Lat.values,
data_df.track.values,
data_df.ground_speed.values,
)
airbornegeo.plotly_points(
data_df[::50],
color_col="eotvos_correction",
)
data_df.head()
[8]:
| EotvosCor | Lon | Lat | Height_WGS1984 | easting | northing | unixtime | line | geometry | index | relative_distance | ground_speed | track | eotvos_correction | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 68.15 | 77.252450 | -80.583923 | 4156.1 | 1.000024e+06 | 226237.330771 | 1.229507e+09 | 1 | POINT (1000023.79 226237.331) | 0 | 0.000000 | NaN | NaN | NaN |
| 1 | 68.37 | 77.252672 | -80.583377 | 4156.0 | 1.000083e+06 | 226246.631269 | 1.229507e+09 | 1 | POINT (1000082.905 226246.631) | 1 | 59.842447 | 59.842447 | 3.805916 | 65.687502 |
| 2 | 68.37 | 77.252901 | -80.582831 | 4156.1 | 1.000142e+06 | 226255.809132 | 1.229507e+09 | 1 | POINT (1000142.048 226255.809) | 2 | 59.850955 | 59.850955 | 3.925778 | 66.003487 |
| 3 | 68.16 | 77.253131 | -80.582285 | 4156.4 | 1.000201e+06 | 226264.969079 | 1.229507e+09 | 1 | POINT (1000201.195 226264.969) | 3 | 59.852244 | 59.852244 | 3.943093 | 66.049747 |
| 4 | 68.12 | 77.253358 | -80.581740 | 4156.6 | 1.000260e+06 | 226274.156809 | 1.229507e+09 | 1 | POINT (1000260.224 226274.157) | 4 | 59.739529 | 59.739529 | 3.899163 | 65.711171 |
2.6. Compare the results#
[11]:
data_df_full = data_df_full.drop(columns=["eotvos_correction"], errors="ignore")
data_df_full = data_df_full.merge(
data_df[["index", "eotvos_correction"]], on="index", how="left"
)
data_df_full["eotvos_diff"] = (
data_df_full["EotvosCor"] - data_df_full["eotvos_correction"]
)
sns.histplot((data_df_full["eotvos_diff"]), kde=True)
plt.title("Eotvos correction difference with published values")
plt.show()
[12]:
airbornegeo.plotly_points(
data_df_full[::50],
color_col="eotvos_diff",
)
From above we can see we have fairly-accurately recreated the Eotvos correction values from the AGAP survey. Below we show a few other methods for computing the correction.
2.7. Eotvos correction following Harlan 1968#
Another common method to compute the correction is from Harlan 1968. The paper offers several equation, but we follow equation 16, which uses ground speed instead of air speed.
[15]:
data_df["eotvos_correction"] = airbornegeo.eotvos_correction_harlan(
data_df.Lat.values,
data_df.track.values,
data_df.ground_speed.values,
data_df.Height_WGS1984.values,
)
airbornegeo.plotly_points(
data_df[::50],
color_col="eotvos_correction",
)
data_df.head()
[15]:
| EotvosCor | Lon | Lat | Height_WGS1984 | easting | northing | unixtime | line | geometry | index | relative_distance | ground_speed | track | eotvos_correction | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 68.15 | 77.252450 | -80.583923 | 4156.1 | 1.000024e+06 | 226237.330771 | 1.229507e+09 | 1 | POINT (1000023.79 226237.331) | 0 | 0.000000 | NaN | NaN | NaN |
| 1 | 68.37 | 77.252672 | -80.583377 | 4156.0 | 1.000083e+06 | 226246.631269 | 1.229507e+09 | 1 | POINT (1000082.905 226246.631) | 1 | 59.842447 | 59.842447 | 3.805916 | 18.438218 |
| 2 | 68.37 | 77.252901 | -80.582831 | 4156.1 | 1.000142e+06 | 226255.809132 | 1.229507e+09 | 1 | POINT (1000142.048 226255.809) | 2 | 59.850955 | 59.850955 | 3.925778 | 19.021738 |
| 3 | 68.16 | 77.253131 | -80.582285 | 4156.4 | 1.000201e+06 | 226264.969079 | 1.229507e+09 | 1 | POINT (1000201.195 226264.969) | 3 | 59.852244 | 59.852244 | 3.943093 | 19.107004 |
| 4 | 68.12 | 77.253358 | -80.581740 | 4156.6 | 1.000260e+06 | 226274.156809 | 1.229507e+09 | 1 | POINT (1000260.224 226274.157) | 4 | 59.739529 | 59.739529 | 3.899163 | 18.859979 |
[ ]: