gps - 如何计算以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();
            }
        }
    }
}