如何计算以 GPS 坐标为中心的地球上的圆上的点? [英] How to calculate points on a circle on the globe centred on GPS coordinates?

查看:18
本文介绍了如何计算以 GPS 坐标为中心的地球上的圆上的点?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

在 KML 中画一个圆

Draw a circle in KML

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

How do you take the GPS coordinates of a point on the globe (say in decimal degree format) and generate the coordinates for a polygon approximating a circle centred on that point?

包含 20 多个数据点的多边形看起来像一个圆.数据点越多 - 圆圈越漂亮.

A polygon with 20+ data points looks like a circle. The more data points - the better looking the circle.

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

I am writing a program that will generate KML and dont know how to calculate the coordinates of the polygon vertices.

数据输入示例:

纬度、经度、圆半径(以英尺为单位)、NumberOfDataPoints

Latitude, Longitude, Circle radius (in feet), NumberOfDataPoints

26.128477, -80.105149, 500, 20

26.128477, -80.105149, 500, 20

推荐答案

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

I don't know if this is the simplest solution and it assumes the world is a sphere.

定义:

R 是球体(即地球)的半径.

R is the radius of the sphere (i.e. the earth).

r 是圆的半径(单位相同).

r is the radius of the circle (in the same units).

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

t is the angle subtended by a great-circle arc of length r at the centre of the sphere so t=r/R radians.

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

Now suppose the sphere has radius 1 and is centred at the origin.

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

C is a unit vector representing the centre of the circle.

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

Imagine a circle round the North pole and consider the point where the plane of the circle intersects the line from the centre of the earth to the North pole. Clearly this point will be somewhere below the North pole.

K 是低于"C 的对应点(即圆的平面与 C 相交的位置),因此 K=cos(t)C

K is the corresponding point "below" C (i.e. where the plane of your circle intersects C) so K=cos(t)C

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

s is the radius of the circle measured in 3D space (i.e. not on the sphere) so s=sin(t)

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

Now we want points on the circle in 3D space with centre K, radius s and lying in the plane passing through and perpendicular to K.

这个答案(忽略旋转的东西)解释了如何找到平面的基向量(即正交向量到正常的 K 或 C).使用叉积求第二个.

This answer (ignore the rotation stuff) explains how to find a basis vector for the plane (i.e. a vector orthogonal to the normal K or C). Use the cross product to find a second.

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

Call these basis vectors U and 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))
}

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

Convert each point to spherical coordinates and you are done.

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

Being bored, I coded this up in C#. The results are plausible: they are in a circle and lie on the sphere. Most of the code implements a struct representing a vector. The actual calculation is very simple.

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();
            }
        }
    }
}

这篇关于如何计算以 GPS 坐标为中心的地球上的圆上的点?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆