I'm trying to find out the missing lat and lon points for a oil pipeline so that it can be 3D-plotted for our project.I am using a third-party api which provides me with the starting lat and long of the pipeline ,it also provides the azimuth,inclination and MD for the whole lenght of the pipeline. I have used the below code to calculate the missing lat/lon points:
import geopandas as gpd
import pandas as pd
import numpy as np
import math
tempdf=pd.DataFrame({'drift_angle':[0, 0.1, 0.4, 0.4, 0.2, 0.2, 0.1, 0.2, 0.8, 0.4, 0.1, 0.3, 3.9, 5, 8.8, 9.9, 10.8, 10.7,9.7, 8.6, 8.1, 6.2, 7.6, 8.2, 9, 8.4, 6.5, 8.1, 8.3, 8.3, 8.8, 8.7, 6.8, 5.3, 3.9, 0.4, 0.8, 0.3, 0.3,8.2, 17.4, 25.7, 33.4, 39.5, 46.9, 55.8, 64.9, 71.5, 80.9, 91.8, 89.8, 90.5, 91.2, 90.3, 91, 91.2, 91.4,90.6, 91.1, 89.9, 90.9, 92, 92.1, 92.4, 92.2, 89.6, 89.5, 89.5, 91.5, 92, 92.4, 90.4, 90.2, 90.5, 89.2,88.8, 88.7, 88.2, 88.6, 90.9, 91.9, 90.9, 90.8, 91.1, 90.9, 90.5, 90.8, 90.2, 91.5, 91.8, 92.3, 92.5,91.1, 90.9, 91.3, 91.5, 89.4, 88.2, 88.3, 88.3, 88.3, 88.3],
'drift_direction':[0, 247.4, 68.5, 48, 42.8, 77.1, 105.1, 143.9, 180.9, 161.7, 160.2, 236.5, 260.3, 254.3,255.7, 260.4, 260.9, 251, 253.9, 258.6, 249.3, 248.1, 241.8, 243, 247, 250.1, 250.6, 253, 253.3, 254.2,255.1, 245.9, 241.9, 237.7, 249, 314.5, 41.5, 42.4, 322.4, 304.7, 309.4, 299.3, 296.4, 297.4, 300.9,300.4, 302.1, 300.6, 300, 300.1, 298.4, 298.2, 297.9, 296.6, 299.4, 300.4, 299, 296.4, 299.2, 298.8,300.9, 303.5, 300.9, 300.8, 300.6, 299.5, 299.1, 299, 302.7, 302.5, 302.5, 300.6, 300.7, 302.4, 299.5,299.3, 299.2, 298.8, 299.5, 303, 304.5, 303.1, 302.1, 299.7, 298.5, 297.9, 298.1, 296.8, 298.2, 299.5,299.1, 298.6, 299.2, 299.7, 299.1, 298.6, 301.2, 301, 300.5, 300.5, 300.2, 300.2]})
tempdf['Inclination'] = np.radians(tempdf['drift_angle'])
tempdf['Azimuth'] = np.radians(tempdf['drift_direction'])
tdf =tempdf.copy()
#print(tempdf)
#print(tdf)
tdf['surface_longitude']=-99.550764
tdf['surface_latitude']=28.443821
tdf['measured_depth']=[0, 229, 410, 593, 776, 959, 1143, 1327, 1511, 1695, 1786, 1963, 2146, 2329, 2513, 2696, 2879, 3063, 3246, 3430, 3613, 3796, 3979, 4163, 4347, 4530, 4714, 4897, 5080, 5264, 5447, 5630, 5814, 5997, 6181, 6365, 6549, 6653, 6745, 6837, 6929, 7021, 7113, 7205, 7297, 7389, 7480, 7572, 7664, 7756, 7848, 7939, 8031, 8123, 8215, 8307, 8398, 8490, 8582, 8674, 8766, 8858, 8949, 9041, 9133, 9225, 9316, 9408, 9500, 9592, 9684, 9776, 9867, 9959, 10051, 10143, 10235, 10326, 10417, 10509, 10601, 10693, 10785, 10877, 10969, 11060, 11152, 11244, 11336, 11428, 11520, 11611, 11703, 11794, 11886, 11978, 12070, 12161, 12253, 12345, 12399, 12466]
#tdf.insert(9,'measured_depth',[0, 229, 410, 593, 776, 959, 1143, 1327, 1511, 1695, 1786, 1963, 2146, 2329, 2513, 2696, 2879, 3063, 3246, 3430, 3613, 3796, 3979, 4163, 4347, 4530, 4714, 4897, 5080, 5264, 5447, 5630, 5814, 5997, 6181, 6365, 6549, 6653, 6745, 6837, 6929, 7021, 7113, 7205, 7297, 7389, 7480, 7572, 7664, 7756, 7848, 7939, 8031, 8123, 8215, 8307, 8398, 8490, 8582, 8674, 8766, 8858, 8949, 9041, 9133, 9225, 9316, 9408, 9500, 9592, 9684, 9776, 9867, 9959, 10051, 10143, 10235, 10326, 10417, 10509, 10601, 10693, 10785, 10877, 10969, 11060, 11152, 11244, 11336, 11428, 11520, 11611, 11703, 11794, 11886, 11978, 12070, 12161, 12253, 12345, 12399, 12466])
gdf = tdf.set_geometry(gpd.points_from_xy(tdf.surface_longitude, tdf.surface_latitude), crs="EPSG:4326")
projected_df = gdf.to_crs("EPSG:32040")
projected_df['LON_32040'] = projected_df.geometry.x
projected_df['LAT_32040'] = projected_df.geometry.y
projected_df.sort_values(by=['measured_depth'])
#print(projected_df.measured_depth)
l = projected_df.index.values
for i in range(len(l)):
if (i == 0):
X1 = projected_df.loc[l[i], 'LON_32040']
Y1 = projected_df.loc[l[i], 'LAT_32040']
projected_df.loc[l[i], 'Longitude'] = X1
projected_df.loc[l[i], 'Latitude'] = Y1
else:
dMD = projected_df.loc[l[i], 'measured_depth'] - projected_df.loc[l[i-1], 'measured_depth']
B = math.acos(math.cos(projected_df.loc[l [i], 'Inclination'] - projected_df.loc[l[i-1], 'Inclination']) - (math.sin(projected_df.loc[l[i-1], 'Inclination'])*math.sin(projected_df.loc[l[i], 'Inclination'] )*(1-math.cos(projected_df.loc[l[i], 'Azimuth']-projected_df.loc[l[i-1], 'Azimuth']))))
if B == 0.0:
RF = 1
else:
RF = 2 / B * math.tan(B / 2)
dX = dMD/2 * (math.sin(projected_df.loc[l[i-1], 'Inclination'])*math.sin(projected_df.loc[l[i-1], 'Azimuth']) + math.sin(projected_df.loc[l[i], 'Inclination'] )*math.sin(projected_df.loc[l[i], 'Azimuth']))*RF
dY = dMD/2 * (math.sin(projected_df.loc[l[i-1], 'Inclination'])*math.cos(projected_df.loc[l[i-1], 'Azimuth']) + math.sin(projected_df.loc[l[i], 'Inclination'] )*math.cos(projected_df.loc[l[i], 'Azimuth']))*RF
X2 = X1 + dX
Y2 = Y1 + dY
projected_df.loc[l[i], 'Longitude'] = X2
projected_df.loc[l[i], 'Latitude'] = Y2
The surface_longitude=-99.550764 and surface_latitude=28.443821 corresponds to the starting lat/lon of the pipeline.
The code first converts the starting lat and lon values to feet(using geopandas) and the azimuth and inclination values to radians, MD is already in feet. Then after that the formula is used to calculate the missing lat/lon values in feet.
After that below code converts it back to degrees.
gdf = projected_df.set_geometry(gpd.points_from_xy(projected_df.Longitude, projected_df.Latitude), crs="EPSG:32040")
projected_df = gdf.to_crs("EPSG:4326")
projected_df['LON_4326'] = projected_df.geometry.x
projected_df['LAT_4326'] = projected_df.geometry.y
After using this approach I wasn't able to acheieve the correct values of lat and lon. I'm not sure if it is a formula issue or something wrong in my implementation.Requesting any advice/suggestions on how to get to the correct values.
