如何计算以 GPS 坐标为中心的地球上的圆上的点？

Question

在 KML 中画一个圆

如何获取地球上某个点的 GPS 坐标（例如十进制度格式）并生成一个多边形的坐标，该多边形近似于以该点为中心的圆？

具有 20 多个数据点的多边形看起来像一个圆。数据点越多 - 圆圈越漂亮。

我正在编写一个会生成 KML 的程序，但不知道如何计算多边形顶点的坐标。

数据输入示例：

纬度、经度、圆半径（英尺）、NumberOfDataPoints

26.128477, -80.105149, 500, 20

Answer 1

我不知道这是否是最简单的解决方案，它假设世界是一个球体。

定义：

R 是球体（即地球）的半径。

r 是圆的半径（单位相同）。

t 是一个长度为 r 的大圆弧在球体中心所对的角度，因此 t=r/R 弧度。

现在假设球体的半径为 1 并且以原点为中心。

C是表示圆心的单位向量。

想象一个围绕北极的圆，并考虑圆的平面与从地球中心到北极的线相交的点。很明显，这个点将位于北极之下的某个地方。

K 是“低于”C 的对应点（即圆的平面与 C 相交的位置），因此 K=cos(t)C

s 是在 3D 空间（即不在球体上）测量的圆的半径，因此 s=sin(t)

现在我们想要 3D 空间中圆上的点，圆心为 K，半径为 s，并且位于通过并垂直于 K 的平面内。

This answer（忽略旋转的东西）解释了如何找到平面的基向量（即与法线 K 或 C 正交的向量）。使用叉积求第二个。

将这些基向量称为 U 和 V。

// Pseudo-code to calculate 20 points on the circle
for (a = 0; a != 360; a += 18)
{
    // A point on the circle and the unit sphere
    P = K + s * (U * sin(a) + V * cos(a))
}

将每个点转换为球坐标就完成了。

因为无聊，我用 C# 编写了这个代码。结果是合理的：它们在一个圆圈中并位于球体上。大多数代码都实现了一个struct 代表一个向量。实际计算很简单。

using System;

namespace gpsCircle
{
    struct Gps
    {
        // In degrees
        public readonly double Latitude;
        public readonly double Longtitude;

        public Gps(double latitude, double longtitude)
        {
            Latitude = latitude;
            Longtitude = longtitude;
        }

        public override string ToString()
        {
            return string.Format("({0},{1})", Latitude, Longtitude);
        }

        public Vector ToUnitVector()
        {
            double lat = Latitude / 180 * Math.PI;
            double lng = Longtitude / 180 * Math.PI;

            // Z is North
            // X points at the Greenwich meridian
            return new Vector(Math.Cos(lng) * Math.Cos(lat), Math.Sin(lng) * Math.Cos(lat), Math.Sin(lat));
        }
    }

    struct Vector
    {
        public readonly double X;
        public readonly double Y;
        public readonly double Z;

        public Vector(double x, double y, double z)
        {
            X = x;
            Y = y;
            Z = z;
        }

        public double MagnitudeSquared()
        {
            return X * X + Y * Y + Z * Z;
        }

        public double Magnitude()
        {
            return Math.Sqrt(MagnitudeSquared());
        }

        public Vector ToUnit()
        {
            double m = Magnitude();

            return new Vector(X / m, Y / m, Z / m);
        }

        public Gps ToGps()
        {
            Vector unit = ToUnit();
            // Rounding errors
            double z = unit.Z;
            if (z > 1)
                z = 1;

            double lat = Math.Asin(z);

            double lng = Math.Atan2(unit.Y, unit.X);

            return new Gps(lat * 180 / Math.PI, lng * 180 / Math.PI);
        }

        public static Vector operator*(double m, Vector v)
        {
            return new Vector(m * v.X, m * v.Y, m * v.Z);
        }

        public static Vector operator-(Vector a, Vector b)
        {
            return new Vector(a.X - b.X, a.Y - b.Y, a.Z - b.Z);
        }

        public static Vector operator+(Vector a, Vector b)
        {
            return new Vector(a.X + b.X, a.Y + b.Y, a.Z + b.Z);
        }

        public override string ToString()
        {
            return string.Format("({0},{1},{2})", X, Y, Z);
        }

        public double Dot(Vector that)
        {
            return X * that.X + Y * that.Y + Z * that.Z;
        }

        public Vector Cross(Vector that)
        {
            return new Vector(Y * that.Z - Z * that.Y, Z * that.X - X * that.Z, X * that.Y - Y * that.X);
        }

        // Pick a random orthogonal vector
        public Vector Orthogonal()
        {
            double minNormal = Math.Abs(X);
            int minIndex = 0;
            if (Math.Abs(Y) < minNormal)
            {
                minNormal = Math.Abs(Y);
                minIndex = 1;
            }
            if (Math.Abs(Z) < minNormal)
            {
                minNormal = Math.Abs(Z);
                minIndex = 2;
            }

            Vector B;
            switch (minIndex)
            {
                case 0:
                    B = new Vector(1, 0, 0);
                    break;
                case 1:
                    B = new Vector(0, 1, 0);
                    break;
                default:
                    B = new Vector(0, 0, 1);
                    break;
            }

            return (B - minNormal * this).ToUnit();
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            // Phnom Penh
            Gps centre = new Gps(11.55, 104.916667);

            // In metres
            double worldRadius = 6371000;
            // In metres
            double circleRadius = 1000;

            // Points representing circle of radius circleRadius round centre.
            Gps[] points  = new Gps[20];

            CirclePoints(points, centre, worldRadius, circleRadius);
        }

        static void CirclePoints(Gps[] points, Gps centre, double R, double r)
        {
            int count = points.Length;

            Vector C = centre.ToUnitVector();
            double t = r / R;
            Vector K = Math.Cos(t) * C;
            double s = Math.Sin(t);

            Vector U = K.Orthogonal();
            Vector V = K.Cross(U);
            // Improve orthogonality
            U = K.Cross(V);

            for (int point = 0; point != count; ++point)
            {
                double a = 2 * Math.PI * point / count;
                Vector P = K + s * (Math.Sin(a) * U + Math.Cos(a) * V);
                points[point] = P.ToGps();
            }
        }
    }
}

Answer 2

我已经编写了Polycircles，这是一个用 Python 实现的小型开源包。它使用geographiclib 进行地理空间计算。