0

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.

9
  • 3
    Welcome to Geographic Information Systems! Please do not ever use EPSG:3857 (Web Mercator) to carry out any sort of distance related analysis. The longitude (X) axis is distorted to infinity as you move north or south of the Equator.
    – Ian Turton
    Sep 7 at 14:11
  • To compliment Ian's comment, your dataset needs to be in metres or feet. So reprojected your data.
    – Hornbydd
    Sep 7 at 14:23
  • @IanTurton the calculation isnt being done in EPSG:3857 though - its in lat-long using Haversine distance. Which should return lat-long points of a circle of 1km around the centre. Which should look like a circle in this projection - otherwise the fountain in the centre would look oval too, wouldn't it? I'd see what you get for a point on the equator, and make sure the formula is correctly coded - test some points of well-known separation in lat and in long.
    – Spacedman
    Sep 7 at 14:51
  • but Google Maps is using EPSG:3857
    – Ian Turton
    Sep 7 at 14:57
  • 3
    A circle almost never looks like a circle in most projections. Web Mercator is an awful projection for the UK because the high latitude distorts significantly. The images used in Maps are processed so they don't look awful, so your basis of comparison is faulty.
    – Vince
    Sep 7 at 15:09
0

To generate the points of a circle in decimal degrees from a center point and radius in meters, you can use shapely and pyproj for the reprojection, and export the circle to KML using simplekml.

The azimuthal equidistant projection (aka aeqd) is an azimuthal map projection that provides points on the map that are at proportionally correct distances from the center point. Create circle in aeqd coordinate space then reproject to WGS-84.

from shapely.geometry import Point
from pyproj import Transformer
from shapely.ops import transform
import simplekml

# Buckingham Palace, London
lat = 51.5014
lng = -0.1419
radius = 1000   # radius in meters

local_azimuthal_projection = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(
    lat, lng
)
wgs84_to_aeqd = Transformer.from_proj('+proj=longlat +datum=WGS84 +no_defs',
                                      local_azimuthal_projection)
aeqd_to_wgs84 = Transformer.from_proj(local_azimuthal_projection,
                                      '+proj=longlat +datum=WGS84 +no_defs')
# Get polygon with lat lon coordinates
point_transformed = Point(wgs84_to_aeqd.transform(lng, lat))

buffer = point_transformed.buffer(radius)
circle = transform(aeqd_to_wgs84.transform, buffer)

# Next export as KML

kml = simplekml.Kml()
ls = kml.newlinestring(tessellate=1,
                       coords=list(circle.exterior.coords),
                       altitudemode=simplekml.AltitudeMode.clamptoground)
ls.style.linestyle.color = 'ffff00aa'  # purple
ls.style.linestyle.width = 3
kml.save("circle.kml")

Output:

View of circle in Google Earth

1
  • nice solution - thanks!
    – c0nt0s0
    Sep 16 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.