Python – How to Generate a Circle of Coordinates Around a Point Correctly

haversinepython

I'm using the following code to attempt to generate a circle of coordinates a fixed distance (in this case 1km) around a point. I'm basing the formula on Haversine, but the output I'm seeing in Google Earth is an Oval, and not a Circle.

import simplekml
import math

kml = simplekml.Kml()

def GetCirclePointCoordinate(point, bearing):
    
    lat = point[0]
    lon = point[1]
    d = point[2]

    R = 6371
    brng = bearing * (math.pi / 180)
    lat1 = lat * (math.pi / 180)
    lon1 = lon * (math.pi / 180)

    lat2 = (math.asin(math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng)))*(180/math.pi)
    lon2 = (lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),math.cos(d/R)-math.sin(lat1)*math.sin(lat2)))*(180/math.pi)

    return lon2,lat2

point = (51.5014,-0.1419,1) # Buckingham Palace, London

circle_coords = []
for deg in range(360):
    circle_coords.append(GetCirclePointCoordinate(point,deg))

ls = kml.newlinestring(name="radius", description="radius")
ls.coords=circle_coords
kml.save("circle.kml")

The output I'm seeing is depicted below:

Oval around Buckingham Palace

I was looking at a solution in this post but wondered if there was a way to do this using a more native solution without the various library imports.

Best Answer

Another direct geometrical approach, without having to play with projections is to rely on libraries intended to do geodesics computation.

The geographiclib is the go to solution for this. It has bindings in many languages.

Python binding

The geodesic Direct and Inverse problem descriptions here: https://geographiclib.sourceforge.io/html/python/geodesics.html#solution-of-geodesic-problems

Quoting:

Traditionally two geodesic problems are considered:

the direct problem — given φ1, λ1, α1, s12, determine φ2, λ2, and α2; this is solved by Geodesic.Direct.

the inverse problem — given φ1, λ1, φ2, λ2, determine s12, α1, and α2; this is solved by Geodesic.Inverse.

(φ_i, λ_i, are spherical angles, s_ij segments length, α_i are azimuth)

A good example is this python library to plot circles: https://polycircles.readthedocs.io/en/latest/

Implementation details:

  • determine nb of circle vertices N for the circle
  • convert that into N bearings (from 0 to 360°)
  • for loop: compute N place/bearing/distance point from center (known as the geodesic Direct problem (geodesic.Geodesic.WGS84.Direct)
  • voilà: you've got your N points that are truely equidistant from center
  • now, put that in a line or polygon in your favorite plotting tool (eg: LineString in google earth)

see the code source from: https://github.com/adamatan/polycircles/blob/main/polycircles/polycircles.py

quoting:

class Polycircle(Shape):
    """A polygonial approximation of a circle in WGS84 coordinates.
    >>> import polycircles
    >>> number_of_vertices = 20
    >>> polycircle = polycircles.Polycircle(latitude=31.611878, longitude=34.505351, radius=100, number_of_vertices=number_of_vertices)
    >>> len(polycircle.to_lat_lon()) == number_of_vertices + 1
    True
    """

    def __init__(self, latitude, longitude, radius, number_of_vertices=DEFAULT_NUMBER_OF_VERTICES):
        """
        Arguments:
        latitude -- WGS84 latitude, between -90.0 and 90.0.
        longitude -- WGS84 longitude, between -180.0 and 180.0.
        radius -- Circle radius in meters.
        number_of_vertices -- Number of vertices for the approximation
        polygon.
        """
        # Value assertions
        assert number_of_vertices >= 3, "The minimal number of vertices in a polygon is 3."
        assert radius > 0, "Radius can only have positive values."
        assert -180 <= longitude <= 180, "Longitude must be between -180 and 180 degrees."
        assert -90 <= latitude <= 90, "Latitude must be between -90 and 90 degrees."

        self.latitude = latitude
        self.longitude = longitude
        self.radius = radius
        self.number_of_vertices = number_of_vertices

        vertices = []
        for i in range(number_of_vertices):
            degree = 360.0/number_of_vertices*i
            vertex = geodesic.Geodesic.WGS84.Direct(latitude, longitude, degree, radius)
            lat = vertex['lat2']
            lon = vertex['lon2']
            vertices.append((lat, lon)) # because it is intended for KML, longitudes come first
        vertices.append(vertices[0])
        self.vertices = tuple(vertices) # tuplification of coords (e.g. for simplekml)

usage here:

https://polycircles.readthedocs.io/en/latest/kmls.html

quoting the usage example:

from polycircles import polycircles

polycircle = polycircles.Polycircle(latitude=40.768085,
                                    longitude=-73.981885,
                                    radius=200,
                                    number_of_vertices=36)
kml = simplekml.Kml()
pol = kml.newpolygon(name="Columbus Circle, Manhattan",
                                         outerboundaryis=polycircle.to_kml())
pol.style.polystyle.color = \
        simplekml.Color.changealphaint(200, simplekml.Color.green)
kml.save("test_kml_polygon_3_manhattan.kml")

NB: polycircles has a small limitation in google earth near poles and date change meridian, this is discussed here https://github.com/adamatan/polycircles/issues/2 with a partial fix proposal (not included in polycircle yet ...)

Related Question