How to calculate points on a circle on the globe centred on GPS coordinates?

Question

Draw a circle in KML

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?

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

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

Example of data inputs:

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

26.128477, -80.105149, 500, 20

Answer 1

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

Define:

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

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

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

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

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 is the corresponding point "below" C (i.e. where the plane of your circle intersects C) so K=cos(t)C

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

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.

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.

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.

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

Answer 2

I have written Polycircles, small open-source package in Python that does it. It uses geographiclib for the geospatial calculation.