问题描述：
平常为了得出地理位置上两点的实际距离（譬如北京与杭州之间的实际距离），除了利用经纬度计算出两点的空间距离，还需要考虑地形因素。由于之间考虑地形造成误差较大，因此采用微分的办法来解决，简单来说就是将两点细分为多点间的距离（当然这个多点是有限的）。
图中计算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 to 3-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 spacea = np.radians(np.array(A))b = np.radians(np.array(B))c = np.radians(np.array(C))# The points in 3D spacea3 = latlong_to_3d(*a)b3 = latlong_to_3d(*b)c3 = latlong_to_3d(*c)# Vectors in 3D spacea3vec = a3 - b3c3vec = c3 - b3# Find the angle between the vectors in 2D spaceangle3deg = 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 pdfrom dbfread import DBFimport osimport numpy as npfrom math import radians, cos, sin, asin, sqrtimport mathdef geodistance(lng1,lat1,lng2,lat2): lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度 dlon=lng2-lng1 dlat=lat2-lat1 a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 distance=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 to 3-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.0 total_distance_plain = 0.0 A = (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 space a3 = latlong_to_3d(*a) b3 = latlong_to_3d(*b) c3 = latlong_to_3d(*c) d3 = latlong_to_3d(*d) # Vectors in 3D space b3vec = b3 - a3 d3vec = d3 - a3 if D == B: angle3deg_a = 0 else: # Find the angle between the vectors in 2D space angle3deg_a = angle_between_vectors_degrees(b3vec,d3vec) # Vectors in 3D space a3vec = a3 - d3 c3vec = c3 - d3 if C == A: angle3deg_d = 0 else: # Find the angle between the vectors in 2D space angle3deg_d = angle_between_vectors_degrees(a3vec,c3vec) if angle3deg_d < 90: angle3deg = angle3deg_a + angle3deg_d else: angle3deg = angle3deg_a + 180 - angle3deg_d x = 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.5 total_distance = total_distance+l#地表距离 total_distance_plain = total_distance_plain+x#弧面距离 total.append((j,total_distance)) print(total)