1

我正在尝试找出输油管道缺少的纬度和经度点，以便可以为我们的项目进行 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.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

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

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

4

0 回答 0