1

2 つの大円の間の交点 (緯度/経度、度) を計算しようとしています (そのため、大円ごとに 1 つずつ、始点と終点を表す 2 つの緯度/経度のペアを入力し、コードは想定されています緯度/経度の交点を返します)。thisやthisなど、他の投稿の手順に従いました。以下に私のコードを投稿します: 定義により、2 つの異なる大円は地球上の 2 つの場所で交差するため、大円のペアごとに 2 つの交点が返されます。残念ながら、間違った結果を返しています (Google Earth のすべてのポイントをプロットして数回テストしただけです)。

import numpy as np
import math

# Define points in great circle 1
p1_lat1 = 32.498520
p1_long1 = -106.816846
p1_lat2 = 38.199999
p1_long2 = -102.371389

# Define points in great circle 2
p2_lat1 = 34.086771
p2_long1 = -107.313379
p2_lat2 = 34.910553
p2_long2 = -98.711786

# Convert points in great circle 1, degrees to radians
p1_lat1_rad = ((math.pi * p1_lat1) / 180.0)
p1_long1_rad = ((math.pi * p1_long1) / 180.0)
p1_lat2_rad = ((math.pi * p1_lat2) / 180.0)
p1_long2_rad = ((math.pi * p1_long2) / 180.0)

# Convert points in great circle 2, degrees to radians
p2_lat1_rad = ((math.pi * p2_lat1) / 180.0)
p2_long1_rad = ((math.pi * p2_long1) / 180.0)
p2_lat2_rad = ((math.pi * p2_lat2) / 180.0)
p2_long2_rad = ((math.pi * p2_long2) / 180.0)

# Put in polar coordinates
x1 = math.cos(p1_lat1_rad) * math.sin(p1_long1_rad)
y1 = math.cos(p1_lat1_rad) * math.cos(p1_long1_rad)
z1 = math.sin(p1_lat1_rad)
x2 = math.cos(p1_lat2_rad) * math.sin(p1_long2_rad)
y2 = math.cos(p1_lat2_rad) * math.cos(p1_long2_rad)
z2 = math.sin(p1_lat2_rad)
cx1 = math.cos(p2_lat1_rad) * math.sin(p2_long1_rad)
cy1 = math.cos(p2_lat1_rad) * math.cos(p2_long1_rad)
cz1 = math.sin(p2_lat1_rad)
cx2 = math.cos(p2_lat2_rad) * math.sin(p2_long2_rad)
cy2 = math.cos(p2_lat2_rad) * math.cos(p2_long2_rad)
cz2 = math.sin(p2_lat2_rad)

# Get normal to planes containing great circles
# np.cross product of vector to each point from the origin
N1 = np.cross([x1, y1, z1], [x2, y2, z2])
N2 = np.cross([cx1, cy1, cz1], [cx2, cy2, cz2])

# Find line of intersection between two planes
L = np.cross(N1, N2)

# Find two intersection points
X1 = L / np.sqrt(L[0]**2 + L[1]**2 + L[2]**2)
X2 = -X1
i_lat1 = math.asin(X1[2]) * 180./np.pi
i_long1 = math.atan2(X1[1], X1[0]) * 180./np.pi
i_lat2 = math.asin(X2[2]) * 180./np.pi
i_long2 = math.atan2(X2[1], X2[0]) * 180./np.pi

# Print results
print i_lat1, i_long1, i_lat2, i_long2

これは を返しますが、交点の 1 つについて(34.314378256910636, -164.51906362183226, -34.314378256910636, 15.480936378167755)値を返す必要があります。(34.280552, -105.549375)

私が間違っていることについてのアイデアはありますか？ご助力ありがとうございます！

編集：将来の参考のために（誰かが同じ問題に遭遇した場合に備えて）、Ruei-Min Linのアドバイスに従って修正したため、44行目と45行目を次のように置き換えました。

N1 = np.cross([y1, x1, z1], [y2, x2, z2])
N2 = np.cross([cy1, cx1, cz1], [cy2, cx2, cz2])
4

1 に答える 1

2

実際、次の設定から開始できます。

(p1_lat1, p1_long1) = (0, 0)
(p1_lat2, p1_long2) = (90, 0)
(p2_lat1, p2_long1) = (0, 0)
(p2_lat2, p2_long2) = (0, 90)

何が起こるかを見るために。結果をトレースすることにより、x1 と y1 を交換する必要があります。

import numpy as np
import math

# Define points in great circle 1
p1_lat1 = 32.498520
p1_long1 = -106.816846
p1_lat2 = 38.199999
p1_long2 = -102.371389

# Define points in great circle 2
p2_lat1 = 34.086771
p2_long1 = -107.313379
p2_lat2 = 34.910553
p2_long2 = -98.711786

# Convert points in great circle 1, degrees to radians
p1_lat1_rad = ((math.pi * p1_lat1) / 180.0)
p1_long1_rad = ((math.pi * p1_long1) / 180.0)
p1_lat2_rad = ((math.pi * p1_lat2) / 180.0)
p1_long2_rad = ((math.pi * p1_long2) / 180.0)

# Convert points in great circle 2, degrees to radians
p2_lat1_rad = ((math.pi * p2_lat1) / 180.0)
p2_long1_rad = ((math.pi * p2_long1) / 180.0)
p2_lat2_rad = ((math.pi * p2_lat2) / 180.0)
p2_long2_rad = ((math.pi * p2_long2) / 180.0)

# Put in polar coordinates
x1 = math.cos(p1_lat1_rad) * math.cos(p1_long1_rad)
y1 = math.cos(p1_lat1_rad) * math.sin(p1_long1_rad)
z1 = math.sin(p1_lat1_rad)
x2 = math.cos(p1_lat2_rad) * math.cos(p1_long2_rad)
y2 = math.cos(p1_lat2_rad) * math.sin(p1_long2_rad)
z2 = math.sin(p1_lat2_rad)
cx1 = math.cos(p2_lat1_rad) * math.cos(p2_long1_rad)
cy1 = math.cos(p2_lat1_rad) * math.sin(p2_long1_rad)
cz1 = math.sin(p2_lat1_rad)
cx2 = math.cos(p2_lat2_rad) * math.cos(p2_long2_rad)
cy2 = math.cos(p2_lat2_rad) * math.sin(p2_long2_rad)
cz2 = math.sin(p2_lat2_rad)

# Get normal to planes containing great circles
# np.cross product of vector to each point from the origin
N1 = np.cross([x1, y1, z1], [x2, y2, z2])
N2 = np.cross([cx1, cy1, cz1], [cx2, cy2, cz2])

# Find line of intersection between two planes
L = np.cross(N1, N2)

# Find two intersection points
X1 = L / np.sqrt(L[0]**2 + L[1]**2 + L[2]**2)
X2 = -X1
i_lat1 = math.asin(X1[2]) * 180./np.pi
i_long1 = math.atan2(X1[1], X1[0]) * 180./np.pi
i_lat2 = math.asin(X2[2]) * 180./np.pi
i_long2 = math.atan2(X2[1], X2[0]) * 180./np.pi

# Print results
print i_lat1, i_long1, i_lat2, i_long2
于 2015-04-06T08:16:17.623 に答える