使用方位角、倾角和测量深度 (MD) 查找井眼/管道的纬度/经度坐标

如何解决使用方位角、倾角和测量深度 (MD) 查找井眼/管道的纬度/经度坐标

我正在尝试找出输油管道缺少的经纬度点，以便为我们的项目绘制 3D 绘图。我正在使用第三方 API，它为我提供了起始经纬度管道，它还提供了管道全长的方位角、倾角和MD。我使用下面的代码来计算缺失的纬度/经度点：

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.2,0.8,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.3,8.7,6.8,5.3,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.4,90.6,91.1,89.9,90.9,92,92.1,92.4,92.2,89.6,89.5,91.5,90.4,90.2,89.2,88.8,88.7,88.2,88.6,91.9,90.8,92.3,92.5,91.3,89.4,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,299,299.2,298.8,303.5,300.8,299.5,299.1,302.7,302.5,300.7,302.4,303,304.5,303.1,299.7,298.5,298.1,296.8,298.6,301.2,301,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,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],'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

surface_longitude=-99.550764 和surface_latitude=28.443821对应于管道的起始纬度/经度。

代码首先将起始纬度和经度值转换为英尺（使用 geopandas），将方位角和倾角值转换为弧度，MD 已经以英尺为单位。然后，该公式用于计算以英尺为单位的缺失纬度/经度值。

之后，下面的代码将其转换回度数。

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

使用这种方法后，我无法获得正确的 lat 和 lon 值。我不确定这是公式问题还是我的实现有问题。请求有关如何获得正确值的任何建议/建议。