地球（球体）表面绕点旋转坐标的实现问题求助

我现在卡在一个经纬度坐标转换的问题上，折腾好几天了实在摸不着头脑，来求助各位大佬！

问题场景 #

我要实现的是：把一个点（红色点）基于红色参考向量的相对位置，映射到蓝色参考向量的坐标系中，最终得到目标蓝色点。操作的范围很小，只有几平方公里，坐标在瑞典斯德哥尔摩附近（大概N59.327 E18.055）。

目前平移和缩放的部分已经搞定了——绿色点是只做了平移缩放的结果，完全符合预期，但旋转环节彻底翻车：我想把绿色点绕蓝色向量的起点旋转到蓝色点的位置，实际得到的却是错误的红色点，轨迹还变成了椭圆；但如果把所有点移到赤道附近，旋转出来的就是正常的圆形轨迹，这明显是球面坐标的问题，但这么小的范围为啥旋转会受影响呢？

我本来觉得小范围用2D计算应该没问题，平移缩放都好好的，结果旋转就翻车了，现在怀疑是用了纯2D的旋转公式，没考虑地球球面的曲率，导致在高纬度地区计算偏差变大，但实在不知道怎么改。

现有代码实现 #

我用C#写了这些类和方法，平移缩放逻辑没问题，就是旋转部分出问题：

internal class DoublePoint 
{
    internal double Latitude;
    internal double Longitude;

    internal DoublePoint() 
    {
        Latitude = 0.0;
        Longitude = 0.0;
    }

    internal DoublePoint(double longitude, double latitude) 
    {
        Latitude = latitude;
        Longitude = longitude;
    }

    public override string ToString() 
    {
        return Longitude.ToString().Replace(",", ".") + ", " + Latitude.ToString().Replace(",", ".");
    }
}

internal class DoubleVector 
{
    internal DoublePoint A;
    internal DoublePoint B;

    internal DoubleVector(DoublePoint a, DoublePoint b) 
    {
        A = a;
        B = b;
    }
}

public static DoublePoint RelativeMove(DoublePoint what, DoubleVector vecFrom, DoubleVector vecTo, double angle = 0.0, bool skipRotate = false) 
{
    double scale = Scale(vecFrom, vecTo);
    DoublePoint p = new DoublePoint();
    // 平移缩放逻辑，这部分得到的绿色点是正确的
    p.Latitude = vecTo.A.Latitude + ((what.Latitude - vecFrom.A.Latitude) * scale);
    p.Longitude = vecTo.A.Longitude + ((what.Longitude - vecFrom.A.Longitude) * scale);

    // 旋转部分出问题，返回的是错误的红色点
    return RotateAroundPoint2(p, vecTo.A, RelativeAngle(vecFrom, vecTo));
}

public static DoublePoint RotateAroundPoint2(DoublePoint pointToRotate, DoublePoint centerDoublePoint, double angleInDegrees) 
{
    double angleInRadians = angleInDegrees * (Math.PI / 180.0);
    double cosTheta = Math.Cos(angleInRadians);
    double sinTheta = Math.Sin(angleInRadians);

    DoublePoint temp = new DoublePoint(
        (cosTheta * (pointToRotate.Longitude - centerDoublePoint.Longitude) - sinTheta * (pointToRotate.Latitude - centerDoublePoint.Latitude) + centerDoublePoint.Longitude),
        (sinTheta * (pointToRotate.Longitude - centerDoublePoint.Longitude) + cosTheta * (pointToRotate.Latitude - centerDoublePoint.Latitude) + centerDoublePoint.Latitude)
    );
    return temp;
}

public static double RelativeAngle(DoubleVector src, DoubleVector dst) 
{
    return Angle(dst) - Angle(src);
}

public static double Angle(DoubleVector vec) 
{
    double longDiff = vec.B.Longitude - vec.A.Longitude;
    double latDiff = vec.B.Latitude - vec.A.Latitude;
    return Math.Atan2(latDiff, longDiff) * (180.0 / Math.PI);
}

public static double Scale(DoubleVector src, DoubleVector dst) 
{
    return Distance(dst) / Distance(src);
}

public static double Distance(DoubleVector vec) 
{
    var sCoord = new GeoCoordinate(vec.A.Latitude, vec.A.Longitude);
    var eCoord = new GeoCoordinate(vec.B.Latitude, vec.B.Longitude);
    return sCoord.GetDistanceTo(eCoord) * 1.2; // 这个1.2是我调试时临时加的系数，暂时可以忽略
}

我试过的方向 #

• 把所有测试点移到赤道附近，旋转结果立刻正常（圆形轨迹），这说明高纬度的经纬度不是均匀的平面坐标，直接套用2D旋转公式肯定有问题
• 考虑过换成DbGeography这类专门处理地理空间的类，但不知道具体怎么对接现有逻辑
• 我猜是不是应该先把经纬度转成局部平面坐标（比如UTM投影），做完旋转平移缩放后再转回经纬度？但我对投影完全不懂，不知道怎么实现

有没有大佬能指点一下，我应该怎么修改旋转部分的代码，或者换哪种思路来实现这个需求？感激不尽！

内容来源于stack exchange