import mathfrom enum import Enumfrom typing import Tuple, Unionclass CoordTransformer: """ 坐标系统转换器:统一处理3D坐标系转换 + GIS投影转换 与C++版本逻辑完全对齐,转换结果100%一致 """ # 外部坐标系配置(位域组合) class EConfig(Enum): NO_CONVERSION = 0 RIGHT_HANDED = (1 << 0) # 外部是右手系 CW_WINDING = (1 << 1) # 外部三角面是CW TEXCOORD_YUP = (1 << 2) # 外部纹理Y向上 UNIT_MM = (1 << 3) # 单位:毫米 UNIT_CM = (2 << 3) # 单位:厘米(内部默认) UNIT_M = (3 << 3) # 单位:米 UNIT_INCH = (4 << 3) # 单位:英寸 UNIT_FEET = (5 << 3) # 单位:英尺 # GIS坐标系WKID class EWKID(Enum): WGS84 = 4326 # WGS84地理坐标系 BJZ54 = 4214 # 北京54地理坐标系 def __init__(self, external_config: EConfig, default_wkid: EWKID = EWKID.WGS84): """ 创建坐标转换器 :param external_config: 外部坐标系的配置(位域组合) :param default_wkid: 默认GIS投影类型 """ self.external_config = external_config.value self.default_wkid = default_wkid # 预计算缓存值,避免重复解析 self.unit_scale = self._get_unit_scale(external_config) self.axis_flip_sign = 1 if self._is_left_handed(external_config) else -1 # ===================== 内部配置解析 ===================== @staticmethod def _get_unit_scale(config: EConfig) -> float: unit_config = EConfig(config.value & (7 << 3)) if unit_config == CoordTransformer.EConfig.UNIT_CM: return 1.0 elif unit_config == CoordTransformer.EConfig.UNIT_MM: return 0.1 elif unit_config == CoordTransformer.EConfig.UNIT_M: return 100.0 elif unit_config == CoordTransformer.EConfig.UNIT_INCH: return 2.54 elif unit_config == CoordTransformer.EConfig.UNIT_FEET: return 30.48 return 1.0 @staticmethod def _is_left_handed(config: EConfig) -> bool: return not (config.value & (1 << 0)) # ===================== 高斯克吕格内部辅助 ===================== class _GKParam: def __init__(self, a: float, f: float): self.a = a self.f = f self.b = a * (1 - f) self.e2 = 2 * f - f * f self.e12 = self.e2 / (1 - self.e2) e4 = self.e2 * self.e2 e6 = e4 * self.e2 e8 = e6 * self.e2 self.A1 = 1 + 3.0/4.0*self.e2 + 45.0/64.0*e4 + 175.0/256.0*e6 + 11025.0/16384.0*e8 self.A2 = 3.0/4.0*self.e2 + 15.0/16.0*e4 + 525.0/256.0*e6 + 2205.0/1024.0*e8 self.A3 = 15.0/64.0*e4 + 105.0/256.0*e6 + 2205.0/4096.0*e8 self.A4 = 35.0/512.0*e6 + 315.0/2048.0*e8 self.A5 = 315.0/16384.0*e8 @staticmethod def _init_gk_param(wkid: EWKID) -> _GKParam: if wkid == CoordTransformer.EWKID.WGS84: return CoordTransformer._GKParam(6378137.0, 1.0 / 298.257223563) else: return CoordTransformer._GKParam(6378245.0, 1.0 / 298.3) # ===================== 3D坐标转换接口 ===================== def PosToInternal(self, src: Union[Tuple[float, float, float], list]) -> Tuple[float, float, float]: """位置:外部 → 内部""" x, y, z = src return ( x * self.unit_scale, y * self.axis_flip_sign * self.unit_scale, z * self.unit_scale ) def PosToExternal(self, src: Union[Tuple[float, float, float], list]) -> Tuple[float, float, float]: """位置:内部 → 外部""" x, y, z = src return ( x / self.unit_scale, y * self.axis_flip_sign / self.unit_scale, z / self.unit_scale ) def RotToInternal(self, src: Union[Tuple[float, float, float, float], list]) -> Tuple[float, float, float, float]: """旋转:外部 → 内部(四元数 x,y,z,w)""" x, y, z, w = src return ( x, y * self.axis_flip_sign, z, w * self.axis_flip_sign ) def RotToExternal(self, src: Union[Tuple[float, float, float, float], list]) -> Tuple[float, float, float, float]: """旋转:内部 → 外部(四元数 x,y,z,w)""" x, y, z, w = src return ( x, y * self.axis_flip_sign, z, w * self.axis_flip_sign ) def TransformToInternal(self, pos: Union[Tuple[float, float, float], list], rot: Union[Tuple[float, float, float, float], list], scale: Union[Tuple[float, float, float], list]) -> Tuple: """变换:外部 → 内部,输入pos/rot/scale,输出转换后的三元组""" internal_pos = self.PosToInternal(pos) internal_rot = self.RotToInternal(rot) internal_scale = ( scale[0] * self.unit_scale, scale[1] * self.unit_scale, scale[2] * self.unit_scale ) return internal_pos, internal_rot, internal_scale def TransformToExternal(self, pos: Union[Tuple[float, float, float], list], rot: Union[Tuple[float, float, float, float], list], scale: Union[Tuple[float, float, float], list]) -> Tuple: """变换:内部 → 外部,输入pos/rot/scale,输出转换后的三元组""" external_pos = self.PosToExternal(pos) external_rot = self.RotToExternal(rot) external_scale = ( scale[0] / self.unit_scale, scale[1] / self.unit_scale, scale[2] / self.unit_scale ) return external_pos, external_rot, external_scale def DirToInternal(self, src: Union[Tuple[float, float, float], list]) -> Tuple[float, float, float]: """方向向量:外部 → 内部""" x, y, z = src return ( x, y * self.axis_flip_sign, z ) def DirToExternal(self, src: Union[Tuple[float, float, float], list]) -> Tuple[float, float, float]: """方向向量:内部 → 外部""" x, y, z = src return ( x, y * self.axis_flip_sign, z ) def UVToExternal(self, src: Union[Tuple[float, float], list]) -> Tuple[float, float]: """纹理UV:内部 → 外部""" u, v = src if self.external_config & self.EConfig.TEXCOORD_YUP.value: return u, 1.0 - v return u, v # ===================== GIS投影转换接口 ===================== def LngLatToXY(self, b_deg: float, l_deg: float, wkid: EWKID = None) -> Tuple[float, float]: """ 经纬度 → 高斯克吕格投影平面坐标 :param b_deg: 纬度(角度) :param l_deg: 经度(角度) :param wkid: 指定坐标系,默认用构造时的配置 :return: (X北方向, Y东方向) 单位:米,Y自带带号+假东 """ if wkid is None: wkid = self.default_wkid param = self._init_gk_param(wkid) b = math.radians(b_deg) l = math.radians(l_deg) # 3度带计算 zone = int(math.floor((l_deg + 1.5) / 3.0)) l0_deg = 3.0 * zone - 1.5 l0 = math.radians(l0_deg) dl = l - l0 sin_b = math.sin(b) cos_b = math.cos(b) tan_b = math.tan(b) N = param.a / math.sqrt(1 - param.e2 * sin_b * sin_b) eta2 = param.e12 * cos_b * cos_b # 子午线弧长 x0 = param.a * (1 - param.e2) * ( param.A1 * b - param.A2 * math.sin(2*b) + param.A3 * math.sin(4*b) - param.A4 * math.sin(6*b) + param.A5 * math.sin(8*b) ) dl2 = dl * dl dl4 = dl2 * dl2 dl6 = dl4 * dl2 t = tan_b t2 = t * t t4 = t2 * t2 # 投影计算 x = x0 + N * t * cos_b*cos_b * ( 0.5 * dl2 + (5 - t2 + 9*eta2 + 4*eta2*eta2) / 24.0 * dl4 + (61 - 58*t2 + t4) / 720.0 * dl6 ) y = N * dl * cos_b * ( 1 + (1 - t2 + eta2) / 6.0 * dl2 + (5 - 18*t2 + t4 + 14*eta2 - 58*eta2*t2) / 120.0 * dl4 ) # 带号+假东 y_out = zone * 1000000.0 + y + 500000.0 return x, y_out def XYToLngLat(self, x: float, y: float, wkid: EWKID = None) -> Tuple[float, float]: """ 高斯克吕格投影平面坐标 → 经纬度 :param x: 北方向坐标(米) :param y: 东方向坐标(米,自带带号+假东) :param wkid: 指定坐标系,默认用构造时的配置 :return: (纬度, 经度) 单位:角度 """ if wkid is None: wkid = self.default_wkid param = self._init_gk_param(wkid) # 解析带号和假东 zone = int(math.floor(y / 1000000.0)) dy = y % 1000000.0 - 500000.0 dx = x l0_deg = 3.0 * zone - 1.5 # 反算纬度 bf = dx / (param.a * (1 - param.e2) * param.A1) b = bf + ( param.A2 * math.sin(2*bf) - param.A3 * math.sin(4*bf) + param.A4 * math.sin(6*bf) - param.A5 * math.sin(8*bf) ) / param.A1 sin_b = math.sin(b) cos_b = math.cos(b) tan_b = math.tan(b) N = param.a / math.sqrt(1 - param.e2 * sin_b * sin_b) eta2 = param.e12 * cos_b * cos_b t = tan_b t2 = t * t t4 = t2 * t2 dy2 = dy * dy dy4 = dy2 * dy2 b = b - t * dy2 / (2*N*N) * ( 1 - (5 + 3*t2 + eta2 - 9*eta2*t2) / 12.0 * dy2/N/N + (61 + 90*t2 + 45*t4) / 360.0 * dy4/N/N/N/N ) dl = dy / (N*cos_b) * ( 1 - (1 + 2*t2 + eta2) / 6.0 * dy2/N/N + (5 + 28*t2 + 24*t4 + 6*eta2 + 8*eta2*t2) / 120.0 * dy4/N/N/N/N ) return math.degrees(b), l0_deg + math.degrees(dl)