问题描述
我正在尝试计算/绘制一条线上的一个点，最靠近另一条线。我可以用笛卡尔坐标和大圆距离等计算该线上的点。我遇到问题的地方是准确计算以英里为单位的距离并绘制它。
我有几行由 lats/lons 定义。我试图将其中一条线的终点作为点 p3。 p1 和 p2 是相邻直线上的两个点。
当我计算交叉轨道距离时，它似乎不正确。 （下图）
附带说明：我可以使用任意数字在笛卡尔坐标系中准确计算。使用我下面的代码：
import numpy as np
import matplotlib.pyplot as plt
x1=1; y1=5
x2=6; y2=7
x0=0; y0=0
p1 = x1,y1
p2 = x2,y2
p3 = x0,y0
d = abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / np.sqrt(np.square(x2-x1) + np.square(y2-y1))
print('d = '+str(d))
A = y1-y2
B = x2-x1
C = (x1*y2)-(x2*y1)
d2 = abs( A*x0 + B*y0 + C ) / np.sqrt(np.square(B) + np.square(A))
print('d2 = '+str(d2))
# point along the line at distance d
# https://en.wikipedia.org/wiki/distance_from_a_point_to_a_line#:~:text=In%20Euclidean%20geometry%2C%20the%20distance,nearest%20point%20on%20the%20line.
x4 = ( B*(B*x0-A*y0) - A*C ) / (A*A + B*B)
y4 = ( A*(A*y0-B*x0) - B*C ) / (A*A + B*B)
P4 = x4;y4
%matplotlib notebook
fig,ax = plt.subplots(figsize=(5,5))
plt.scatter(x1,y1,color="red",marker=".",label="p1")
plt.scatter(x2,y2,color="blue",label="p2")
plt.scatter(x0,y0,color="black",label="p3")
plt.scatter(x4,y4,color="green",label="p4")
# plt.axline((x1,y1),(x2,y2),linestyle='--',color='black',zorder=0)
plt.plot((x0,x4),(y0,y4),linestyle=':',zorder=0)
plt.axis('equal')
plt.legend()
这是剧情：
所以，我可以使方程起作用，我只是在将它们转换为/从笛卡尔/球面/椭圆和获得准确距离时遇到困难。
我知道我正在将这些纬度/经度绘制为笛卡尔图，但它仍然没有在我脑海中加起来。
这是我用来计算/绘图的代码：
在这个网站上找到： https://gis.stackexchange.com/questions/209540/projecting-cross-track-distance-on-great-circle
from math import radians,cos,sin,asin,sqrt,degrees,atan2
def validate_point(p):
lat,lon = p
assert -90 <= lat <= 90,"bad latitude"
assert -180 <= lon <= 180,"bad longitude"
# original formula from http://www.movable-type.co.uk/scripts/latlong.html
def distance_haversine(p1,p2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
haversine
formula:
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
_ ____
c = 2 ⋅ atan2( √a,√(1−a) )
d = R ⋅ c
where φ is latitude,λ is longitude,R is earth’s radius (mean radius = 6,371km);
note that angles need to be in radians to pass to trig functions!
"""
lat1,lon1 = p1
lat2,lon2 = p2
for p in [p1,p2]:
validate_point(p)
R = 6371 # km - earths's radius
# convert decimal degrees to radians
lat1,lon1,lat2,lon2 = map(radians,[lat1,lon2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a)) # 2 * atan2(sqrt(a),sqrt(1-a))
d = R * c # convert to km
d = d * 0.621371 # convert km to miles
return d
def spherical2Cart(lat,lon):
clat=(90-lat)*np.pi/180.
lon=lon*np.pi/180.
x=np.cos(lon)*np.sin(clat)
y=np.sin(lon)*np.sin(clat)
z=np.cos(clat)
return np.array([x,y,z])
def cart2Spherical(x,z):
r=np.sqrt(x**2+y**2+z**2)
clat=np.arccos(z/r)/np.pi*180
lat=90.-clat
lon=np.arctan2(y,x)/np.pi*180
lon=(lon+360)%360
return np.array([lat,lon,np.ones(lat.shape)])
def greatCircle(lat1,lon2,r=None,verbose=True):
'''Compute the great circle distance on a sphere
<lat1>,<lat2>: scalar float or nd-array,latitudes in degree for
location 1 and 2.
<lon1>,<lon2>: scalar float or nd-array,longitudes in degree for
location 1 and 2.
<r>: scalar float,spherical radius.
Return <arc>: great circle distance on sphere.
'''
if r is None:
r=6371 # km
d2r=lambda x:x*np.pi/180
lat1,lon2=map(d2r,lon2])
dlon=abs(lon1-lon2)
numerator=(cos(lat2)*sin(dlon))**2 + \
(cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dlon))**2
numerator=np.sqrt(numerator)
denominator=sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(dlon)
dsigma=np.arctan2(numerator,denominator)
arc=r*dsigma
return arc
def getCrosstrackPoint(lat1,lat3,lon3):
'''Get the closest point on great circle path to the 3rd point
<lat1>,<lon1>: scalar float or nd-array,latitudes and longitudes in
degree,start point of the great circle.
<lat2>,end point of the great circle.
<lat3>,<lon3>: scalar float or nd-array,a point away from the great circle.
Return <latp>,<lonp>: latitude and longitude of point P on the great
circle that connects P1,P2,and is closest
to point P3.
'''
x1,z1=spherical2Cart(lat1,lon1)
x2,z2=spherical2Cart(lat2,lon2)
x3,y3,z3=spherical2Cart(lat3,lon3)
D,E,F=np.cross([x1,z1],[x2,z2])
a=E*z3-F*y3
b=F*x3-D*z3
c=D*y3-E*x3
f=c*E-b*F
g=a*F-c*D
h=b*D-a*E
tt=np.sqrt(f**2+g**2+h**2)
xp=f/tt
yp=g/tt
zp=h/tt
result1=cart2Spherical(xp,yp,zp)
result2=cart2Spherical(-xp,-yp,-zp)
d1=greatCircle(result1[0],result1[1],lon3,r=1)
d2=greatCircle(result2[0],result2[1],r=1)
if d1>d2:
return result2[0],result2[1]
else:
return result1[0],result1[1]
def getCrosstrackdistance(lat1,r=None):
'''Compute cross-track distance
<lat1>,a point away from the great circle.
Return <dxt>: great cicle distance between point P3 to the closest point
on great circle that connects P1 and P2.
NOTE that the sign of dxt tells which side of the 3rd point
P3 is on.
'''
if r is None:
r=CONS.EARTH_RADIUS
# get angular distance between P1 and P3
delta13=greatCircle(lat1,r=1.)
# bearing between P1,P3
theta13=getbearing(lat1,lon3)*np.pi/180
# bearing between P1,P2
theta12=getbearing(lat1,lon2)*np.pi/180
dtheta=np.arcsin(sin(delta13)*sin(theta13-theta12))
dxt=r*dtheta
return dxt
def getAlongTrackdistance(lat1,r=None):
'''Compute the distance from the start point to closest point to 3rd point
<lat1>,a point away from the great circle.
Return <dat>: distance from the start point to the closest point on
the great circle connecting P1 and P2.
See also getCrosstrackdistance(),getCrosstrackPoint().
'''
if r is None:
r=CONS.EARTH_RADIUS
# angular distance from P1 to P3
delta13=greatCircle(lat1,r=1.)
# angular distance from Pcloset to P3
dxt=getCrosstrackdistance(lat1,r=1)
dat=r*np.arccos(cos(delta13)/cos(dxt/r))
return dat
def getbearing(lat1,long1,long2):
import math
dLon = (long2 - long1)
y = math.sin(dLon) * math.cos(lat2)
x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(dLon)
brng = math.atan2(y,x)
brng = np.rad2deg(brng)
return brng
x1 = -75.9671285; y1 = 41.67304238
x2 = -75.96262168; y2 = 41.66425696
x3 = -75.96017288; y3 = 41.67049662
p1=[x1,y1]
p2=[x2,y2]
p3=[x3,y3]
pp=getCrosstrackPoint(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],)
print('pp',pp)
dxt=getCrosstrackdistance(p1[0],r=1)
print('dxt',dxt)
dat=getAlongTrackdistance(p1[0],r=1)
print('dat',dat)
dxt2=greatCircle(pp[0],pp[1],r=1)
print('dxt2',dxt2)
dat2=greatCircle(pp[0],p1[0],r=1)
print('dat2',dat2)
%matplotlib notebook
fig,ax = plt.subplots(figsize=(8,8))
plt.plot(p1[0],marker="o",label="Great Circle Point p1",markersize=12)
plt.plot(p2[0],label="Great Circle Point p2",markersize=12)
plt.plot(p3[0],label="Point p3",markersize=12)
plt.plot(pp[0],color="Black",marker="*",label="Crosstrack Point",markersize=12)
plt.plot([pp[0],p3[0]],[pp[1],p3[1]],label='Cross track')
print('Cross track dist = ',distance_haversine(pp,p3))
for q in range(df.shape[0]):
toe_lat = df['wgs84_toe_lat'].iloc[q]
toe_lon = df['wgs84_toe_lon'].iloc[q]
heel_lat = df['wgs84_heel_lat'].iloc[q]
heel_lon = df['wgs84_heel_lon'].iloc[q]
plt.plot([toe_lon,heel_lon],[toe_lat,heel_lat],'green')
if q==0:
plt.plot([toe_lon,'green',label="Well Lateral")
plt.ylabel('Latitude')
plt.xlabel('Longitude')
plt.grid(True)
plt.axis('equal')
plt.legend()
当我计算从 p1 到 p3 的半正弦距离时，它计算出 0.48 英里，但 GIS 软件显示为 0.4 英里。这几乎没有我需要的那么准确。
这是否意味着我正在评估的线/点非常接近笛卡尔坐标会更准确？如果是这样，我如何使用纬度/经度以英里为单位准确计算/绘制它们？
当我尝试将纬度/经度转换为弧度时
具有类似的功能：
def conv_latlon2xy(lat=None,lon=None):
lat,lon = np.deg2rad(lat),np.deg2rad(lon)
# R = 6371 # radius of the earth
x = R * np.cos(lat) * np.cos(lon)
y = R * np.cos(lat) * np.sin(lon)
z = R *np.sin(lat)
return x,y
我没有得到准确的结果。
在这种情况下，使用纬度/经度计算从点到线的最近距离并能够以英里为单位准确测量/绘制该距离（距离小于 1 英里）的最佳方法是什么？我可以在笛卡尔和球形/椭圆形中计算它，但我无法在两者之间跳转。
解决方法
我发现解决这个问题的最好方法是使用 pyproj 和类似下面的代码：
transformer = pyproj.Transformer.from_crs(4326,2271) # WSG84 -> NAD83PAnorth ft
df['xx_mid'],df['yy_mid'] = transformer.transform(df['WGS84_mid_lat'].copy(),df['WGS84_mid_lon'].copy())
这会将纬度/经度转换为英尺或米。