多个GPS点位相对单个参考点的带正负号距离计算问题

我有80个共线点位，这条线穿过一个特定参考点。我需要计算每个台站到该参考点的距离，要求距离值可正可负。

我尝试了如下方法，但它只能返回正的距离值：

stla = st[0].stats.sac['stla'] # 台站纬度
stlo = st[0].stats.sac['stlo'] # 台站经度
# 参考点坐标
crlo = 104.499765143 
crla = 29.4742139251
dist1 = gps2dist_azimuth(stla,stlo, crla,crlo)
dist1=dist1[0]

我得到的输出如下（单位：米）：

1549.9072755468023
 2314.5301492042786
 1978.8805972408672
 604.7927668562546
 195.07039250670383
 2177.129827481811
 212.6443078339846
 863.9573297833116
 1615.5449194502921
 1219.2825679960995
 1564.226848063852
 2600.316245974806
 741.5826151164387
 665.5972988769541
 2335.6051118556975

我期望的输出应该是如下形式（单位为Km，可忽略单位差异）：

0.5983164169
0.637797120571
0.732341381346  
0.776296768282  
0.83927019696   
0.864898580629  
0.927615280754  
0.9785794379    
1.07543888694   
1.09081523648   
1.14869546005   
1.18406238008   
1.25599403922   
1.33394492802   
1.3999813227    
1.44882008293
1.51207083251
1.57852778971

我知道可以通过如下GMT命令实现，指定一个点位作为中心点：

gmt project 1.txt -C104.499765143/29.4742139251 -A87.64 -Fpq -Q 

其中1.txt为台站位置文件，-C参数指定中心点坐标，-A参数代表方位角。请问我要如何实现需求？

方案1：纯Python原生实现 #

你使用的gps2dist_azimuth（ObsPy自带工具函数）返回值的第二个元素就是台站相对于参考点的方位角（0~360度，正北为0顺时针递增），结合已知的测线方位角87.64即可给距离添加符号：

import math
from obspy.geodetics import gps2dist_azimuth

# 固定参数
crlo = 104.499765143 
crla = 29.4742139251
line_az = 87.64 # 测线方位角，和GMT命令的-A参数一致

for tr in st: # 遍历所有台站
    stla = tr.stats.sac['stla']
    stlo = tr.stats.sac['stlo']
    # 返回值依次为：距离(米)、台站相对参考点方位角、反向方位角
    dist_m, az, baz = gps2dist_azimuth(stla, stlo, crla, crlo)
    # 计算方位差余弦判断符号：余弦为正代表台站在测线正方向，为负代表在反方向
    az_diff = math.radians(az - line_az)
    sign = 1 if math.cos(az_diff) > 0 else -1
    # 输出带符号的千米单位距离
    signed_dist_km = sign * dist_m / 1000
    print(signed_dist_km)

方案2：直接调用GMT命令对齐结果 #

如果需要和GMT输出完全一致，可以直接在Python中调用现有GMT命令：

import subprocess
import os

# 写入台站坐标到临时文件
with open("1.txt", "w") as f:
    for tr in st:
        f.write(f"{tr.stats.sac['stlo']} {tr.stats.sac['stla']}\n")

# 执行GMT命令
cmd = "gmt project 1.txt -C104.499765143/29.4742139251 -A87.64 -Fpq -Q"
result = subprocess.check_output(cmd, shell=True, text=True)

# 解析结果，第一列即为带符号的距离（单位为千米）
signed_dists = [float(line.split()[0]) for line in result.strip().split("\n")]
print(signed_dists)

# 清理临时文件
os.remove("1.txt")
if os.path.exists("gmt.history"):
    os.remove("gmt.history")

内容的提问来源于stack exchange，提问作者Aqeel