import numpy as np from scipy.spatial.transform import Rotation as R import pyproj from pyproj import Proj, transform #0.016938035523210708 0.58455146147157355 -0.48870579156409283 0.64744060819180593 -129342.74756339534 3469822.8668770161 5343696.2087743161 # 四元数 (QW, QX, QY, QZ) qw, qx, qy, qz = 0.016938035523210708, 0.58455146147157355, -0.48870579156409283, 0.64744060819180593 # 替换为实际值 # 平移向量 (TX, TY, TZ) tx, ty, tz = -129342.74756339534, 3469822.8668770161, 5343696.2087743161 # 将四元数转换为旋转矩阵 rotation = R.from_quat([qx, qy, qz, qw]) rotation_matrix = rotation.as_matrix() # 计算相机中心坐标(ECEF坐标系下) camera_center = -np.dot(rotation_matrix.T, np.array([tx, ty, tz])) print("投影中心坐标:(ecef坐标系下)", camera_center) # 定义坐标转换 # 从ECEF转换到WGS84 ecef = Proj(proj='geocent', ellps='WGS84', datum='WGS84') wgs84 = Proj(proj='latlong', ellps='WGS84', datum='WGS84') # 执行转换 lon, lat, alt = transform(ecef, wgs84, camera_center[0],camera_center[1], camera_center[2], radians=False) # 输出WGS84经纬度坐标 print("经度:", lon, "纬度:", lat, "高度:", alt) # 从WGS84转换到UTM Zone 50N utm_zone_50n = Proj(proj="utm", zone=50, ellps='WGS84', datum='WGS84', south=False) utm_x, utm_y = transform(wgs84, utm_zone_50n, lon, lat) # 输出UTM坐标 print("WGS84坐标系下:UTM Zone 50N X:", utm_x, "Y:", utm_y,"Z:",alt) # 将四元数转换为旋转矩阵 rotation = R.from_quat([qx, qy, qz, qw]) rotation_matrix = rotation.as_matrix() # 将旋转矩阵转换为欧拉角 (Omega, Phi, Kappa) # 摄影测量中通常使用 ZYX 旋转顺序 omega, phi, kappa = rotation.as_euler('ZYX', degrees=True) print("旋转矩阵:\n", rotation_matrix) print("Omega:", omega, "Phi:", phi, "Kappa:", kappa)