Python——球面两点距离及两条直线夹角的计算

问题描述：
平常为了得出地理位置上两点的实际距离（譬如北京与杭州之间的实际距离），除了利用经纬度计算出两点的空间距离，还需要考虑地形因素。由于之间考虑地形造成误差较大，因此采用微分的办法来解决，简单来说就是将两点细分为多点间的距离（当然这个多点是有限的）。

图中计算AB之间的距离，可以计算出中间多个点位的距离（如AC），然后计算在AB直线上的投影。在此之前，需要计算出球面上两条直线间的夹角以及两点在球面上的距离。

1）球面上两条直线间的夹角
方法一：简易版，将球面弧线看成是直角坐标系下的直线，采用向量乘积。角度a = arc cos (a1a2/|a1||a2|)
方法二：将球面坐标系转为三维笛卡尔坐标系，然后利用向量乘积求角度。
这里考虑到向量的方向性及准确性，采用方法二会更好点。

def latlong_to_3d(latr,lonr):"""Convert a point given latitude and longitude in radians to3-dimensional space,assuming a sphere radius of one."""return np.array((math.cos(latr) * math.cos(lonr),math.cos(latr) * math.sin(lonr),math.sin(latr)))
def angle_between_vectors_degrees(u,v):"""Return the angle between two vectors in any dimension space,in degrees."""return np.degrees(math.acos(np.dot(u,v) / (np.linalg.norm(u) * np.linalg.norm(v))))
# Convert the points to numpy latitude/longitude radians space
a = np.radians(np.array(A))
b = np.radians(np.array(B))
c = np.radians(np.array(C))# The points in 3D space
a3 = latlong_to_3d(*a)
b3 = latlong_to_3d(*b)
c3 = latlong_to_3d(*c)
# Vectors in 3D space
a3vec = a3 - b3
c3vec = c3 - b3
# Find the angle between the vectors in 2D space
angle3deg = angle_between_vectors_degrees(a3vec,c3vec)

2)两点间的距离
S=R·arc cos[cosβ1cosβ2cos（α1-α2）+sinβ1sinβ2]
具体原理见球面两点距离原理

或是采用google map办法：
地球上任两点，其经度分别为A1、A2（E正，W负），纬度分别为B1、B2（N正，S负)。
令A0=（A1-A2）÷2，B0=（BI-B2）÷2
f=√sinB0×sinB0+cosB1×cosB2×sinA0×sinA0
S = 2Rf

在依次计算逐个站点在直线上的投影，然后考虑两点间的高程数据，进行累加得到AB间实际距离。
附上代码如下：

import pandas as pd
from dbfread import DBF
import os
import numpy as np
from math import radians, cos, sin, asin, sqrt
import mathdef geodistance(lng1,lat1,lng2,lat2):lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度dlon=lng2-lng1dlat=lat2-lat1a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径，6371km#distance=round(distance/1000,3)return distancedef latlong_to_3d(latr,lonr):"""Convert a point given latitude and longitude in radians to3-dimensional space,assuming a sphere radius of one."""return np.array((math.cos(latr) * math.cos(lonr),math.cos(latr) * math.sin(lonr),math.sin(latr)))
def angle_between_vectors_degrees(u,v):"""Return the angle between two vectors in any dimension space,in degrees."""return np.degrees(math.acos(np.dot(u,v) / (np.linalg.norm(u) * np.linalg.norm(v))))path = r'D:\20210208'
filetype = '.dbf'#指定文件类型
total=[]
def get_filename(path,filetype):name =[]final_name = []for root,dirs,files in os.walk(path):for i in files:if filetype in i:name.append(i.replace(filetype,''))#生成不带‘.csv’后缀的文件名组成的列表final_name = [item +'.dbf' for item in name]#生成‘.csv’后缀的文件名组成的列表return final_name#输出由有‘.csv’后缀的文件名组成的列表path_all = get_filename(path,filetype)
for j in path_all:table = DBF(path + '/' + j) #, encoding='GBK')df = pd.DataFrame(iter(table))total_distance = 0.0total_distance_plain = 0.0A = (df['POINT_Y'][0],df['POINT_X'][0])B = (df['POINT_Y'][len(df)-1],df['POINT_X'][len(df)-1])a = np.radians(np.array(A))b = np.radians(np.array(B))for n in range(len(df)-1):C = (df['POINT_Y'][n],df['POINT_X'][n])D = (df['POINT_Y'][n+1],df['POINT_X'][n+1])c = np.radians(np.array(C))d = np.radians(np.array(D))print(j,A,B,C,D)# The points in 3D spacea3 = latlong_to_3d(*a)b3 = latlong_to_3d(*b)c3 = latlong_to_3d(*c)d3 = latlong_to_3d(*d)# Vectors in 3D spaceb3vec = b3 - a3d3vec = d3 - a3if D == B:angle3deg_a = 0else:# Find the angle between the vectors in 2D spaceangle3deg_a = angle_between_vectors_degrees(b3vec,d3vec)# Vectors in 3D spacea3vec = a3 - d3c3vec = c3 - d3if C == A:angle3deg_d = 0else:# Find the angle between the vectors in 2D spaceangle3deg_d = angle_between_vectors_degrees(a3vec,c3vec)if angle3deg_d < 90:angle3deg = angle3deg_a + angle3deg_delse:angle3deg = angle3deg_a + 180 - angle3deg_dx = geodistance(df['POINT_X'][n],df['POINT_Y'][n],df['POINT_X'][n+1],df['POINT_Y'][n+1])h = df['GRID_CODE'][n]-df['GRID_CODE'][n+1]x_t = x*cos(angle3deg)l = (x_t*x_t+h*h)**0.5total_distance = total_distance+l#地表距离total_distance_plain = total_distance_plain+x#弧面距离total.append((j,total_distance))
print(total)

Python——球面两点距离及两条直线夹角的计算相关推荐

使用向量点乘的方式计算两条直线夹角的代码案例（平面内，也可拓展至三维空间）
案例代码: bool JudgeDoubleIsZero(double f){const double EPSINON = 0.000001F;if ((f >= -EPSINON) & ...
opencv求解两条直线的交点
假设现在有一个点集,需要拟合出最能够表达点集轮廓的几条直线,并求直线之间的交点. 从点集中拟合直线可以采用的方法:随机抽样一致性(RANSAC),霍夫变换(though transform) 思路1 ...
求空间两条直线之间的距离
1. 前言最近老板让写一段空间点匹配的代码, 其中涉及到求空间两直线之间的距离,写起来满费劲的, 这里做一个记录. 2. 处理思路空间两直线之间的位置关系主要可以分为: 重合, 平行, 相交, 异 ...
学习OpenCV3：判断两条直线平行，并计算平行距离
一.问题已知两条直线 l 1 ( x 1 , y 1 , x 2 , y 2 ) l_1(x_1,y_1,x_2,y_2) l1(x1,y1,x2,y2)和 l 2 ( x 3 , y ...
算法提高两条直线（二分法 + 枚举）
试题算法提高两条直线资源限制时间限制:1.0s 内存限制:256.0MB 问题描述给定平面上n个点. 求两条直线,这两条直线互相垂直,而且它们与x轴的夹角为45度,并且n个点中离这两条直线的 ...
两条直线求交点c语言,C§ 3.3.1两条直线的交点坐标（5页）-原创力文档
§ 3.1两条直线的交点坐标学习目标 1．掌握判断两直线相交的方法:会求两直线交点坐标: 2.体会判断两直线相交中的数形结合思想. 学习过程一.课前准备: (预习教材P112~ P114,找出疑 ...
空间两条直线段的最短距离及最近点计算
如果这两条直线段不共线,假设直线段l0的两端点为:P0.P1:直线段l1的两端点为Q0.Q1,:求两直线段的最短距离? 直线段l0我们可以用方程表示为: (1) 直线段l1我们也可以用方程表示 ...
空间两条直线的最短距离及最近点计算
直线的信息可以以两个端点的形式给出,也可以以一个直线上的点和直线的方向向量给出.本文中假设这两条直线不共线,即这两条直线既不重合也不相交. 1.如果这两条直线是以两个端点的形式给出,那么假设直线l0的 ...
判断两条直线是否相交c语言,计算几何-两条线段是否相交（三种算法）
原标题:计算几何-两条线段是否相交(三种算法) 计算几何中,判断线段是否相交是最基本的题目. 所谓几何, 最基本的当然就是坐标, 从坐标中我们可以知道位置和方向,比如:一个点就是一个位置,两点确定一条 ...
使用C++面向对象思想计算两条直线交点
使用C++面向对象思想计算两条直线交点以下是使用C++面向对象思想计算两条直线交点的示例代码: #include <iostream>using namespace std;class ...

Python——球面两点距离及两条直线夹角的计算

Python——球面两点距离及两条直线夹角的计算相关推荐

最新文章

热门文章