基于 geopandas 读取 postgis数据 + 求地球两点的方位角-球面距离-欧氏距离

import numpy as np
import nvector as nv

# 基于 geopandas 取 postgis数据
sql = ‘SELECT * FROM node_table‘
g_out = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="geom")
values = g_out.values
sql = ‘SELECT geom FROM node_table WHERE ST_DWithin(geom, \‘POINT(116.3211279 39.984223)\‘, 0.002)‘
g_out2 = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="geom")
point_ndarray = g_out2.values
point = point_ndarray[0, 0]
lon = point.x
lat = point.y

# 求两点的方位角
wgs84 = nv.FrameE(name=‘WGS84‘)
pointA = wgs84.GeoPoint(latitude=39.984323, longitude=116.3211279, z=200, degrees=True)
pointB = wgs84.GeoPoint(latitude=39.984223, longitude=116.3211279, z=200, degrees=True)
# Step1: Find p_AB_N (delta decomposed in N)
p_AB_N = pointA.delta_to(pointB)
# x, y, z = p_AB_N.pvector.ravel()
# valtxt = ‘{0:8.2f}, {1:8.2f}, {2:8.2f}‘.format(x, y, z)
# print(‘Ex1: delta north, east, down = {}‘.format(valtxt))
# Step2: Also find the direction (azimuth) to B, relative to north
azimuth = p_AB_N.azimuth_deg
print(‘azimuth = {0:4.2f} deg‘.format(azimuth))
# 求两点的球面距离和欧氏距离
wgs84 = nv.FrameE(name=‘WGS84‘)
point1 = wgs84.GeoPoint(latitude=39.984323, longitude=116.3211279, degrees=True)
point2 = wgs84.GeoPoint(latitude=39.984223, longitude=116.3211279, degrees=True)
s_12, _azi1, _azi2 = point1.distance_and_azimuth(point2)
p_12_E = point2.to_ecef_vector() - point1.to_ecef_vector()
d_12 = p_12_E.length
msg = ‘Ellipsoidal and Euclidean distance = {:5.2f} km, {:5.2f} km‘.format(s_12 / 1000, d_12 / 1000)
print(msg)