我正在尝试计算/绘制一条线上最靠近另一条线的点。我可以计算那条线上的点，无论是笛卡尔坐标还是大圆距离等。我遇到问题的地方是以英里为单位准确计算距离并绘制它。
我有几行由纬度/经度定义。我试图将其中一条线的终点作为点，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",marker=".",label="p2")
plt.scatter(x0,y0,color="black",marker=".",label="p3")
plt.scatter(x4,y4,color="green",marker=".",label="p4")
# plt.axline((x1, y1), (x2, y2), linestyle='--', color='black', zorder=0)
plt.plot((x0, x4), (y0, y4), linestyle=':', color='black', 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, lon1, lat2, 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,y,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,lon1,lat2,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,lon1,lat2,lon2=map(d2r,[lat1,lon1,lat2,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,lon1,lat2,lon2,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>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, 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,y1,z1=spherical2Cart(lat1,lon1)
x2,y2,z2=spherical2Cart(lat2,lon2)
x3,y3,z3=spherical2Cart(lat3,lon3)
D,E,F=np.cross([x1,y1,z1],[x2,y2,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],lat3,lon3,r=1)
d2=greatCircle(result2[0],result2[1],lat3,lon3,r=1)
if d1>d2:
return result2[0],result2[1]
else:
return result1[0],result1[1]
def getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute cross-track distance
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, 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,lon1,lat3,lon3,r=1.)
# bearing between P1, P3
theta13=getBearing(lat1,lon1,lat3,lon3)*np.pi/180
# bearing between P1, P2
theta12=getBearing(lat1,lon1,lat2,lon2)*np.pi/180
dtheta=np.arcsin(sin(delta13)*sin(theta13-theta12))
dxt=r*dtheta
return dxt
def getAlongTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute the distance from the start point to closest point to 3rd point
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, 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,lon1,lat3,lon3,r=1.)
# angular distance from Pcloset to P3
dxt=getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=1)
dat=r*np.arccos(cos(delta13)/cos(dxt/r))
return dat
def getBearing(lat1, long1, lat2, 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],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print('dxt',dxt)
dat=getAlongTrackDistance(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print('dat',dat)
dxt2=greatCircle(pp[0],pp[1],p3[0],p3[1],r=1)
print('dxt2',dxt2)
dat2=greatCircle(pp[0],pp[1],p1[0],p1[1],r=1)
print('dat2',dat2)
%matplotlib notebook
fig,ax = plt.subplots(figsize=(8,8))
plt.plot(p1[0],p1[1],color="blue",marker="o",label="Great Circle Point p1",markersize=12)
plt.plot(p2[0],p2[1],color="green",marker="o",label="Great Circle Point p2",markersize=12)
plt.plot(p3[0],p3[1],color="red",marker="o",label="Point p3",markersize=12)
plt.plot(pp[0],pp[1],color="Black",marker="*",label="Crosstrack Point",markersize=12)
plt.plot([pp[0],p3[0]],[pp[1],p3[1]],linestyle='--',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,heel_lon],[toe_lat,heel_lat],'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 英里的最佳方法是什么？我可以用笛卡尔和球形/椭圆来计算它，但我不能在两者之间进行跳跃。