Generating polygon representing rough 100km circle around latitude/longitude point using Python?

Question

Given a latitude/longitude point anywhere on earth, I would like to assemble a polygon that roughly resembles a circle with a 100km radius around that point.

It doesn't have to be pin-point accurate, just good to within a few hundred meters. I'd like to control the number of points returned (20 is probably fine).

I'm already using the Shapely library for Python, but it looks like that can't quite get the job done on its own.

My ideal function would look something like this:

polygon = circle_around_point(
    latitude=29.9500,
    longitude=-90.0667,
    radius_in_meters=1000,
    num_points=20
)

polygon would then be a list of 20 lat/lon points around the edge of that circle.

Answer 1

I think I've found a working solution, using the Python geog library: https://pypi.python.org/pypi/geog

import numpy as np
import json
import geog
import shapely.geometry
p = shapely.geometry.Point([-90.0667, 29.9500])

n_points = 20
d = 10 * 1000  # meters
angles = np.linspace(0, 360, n_points)
polygon = geog.propagate(p, angles, d)
print(json.dumps(shapely.geometry.mapping(shapely.geometry.Polygon(polygon))))

Then paste the resulting GeoJSON into http://geojson.io/ to preview the result. This seems to do what I'm looking for.

Answer 2

An alternative solution would be to create a local metric projection and create the buffer in that projection.

https://shapely.readthedocs.io/en/stable/manual.html#other-transformations tells you how to transform Shapely geometries

https://en.wikipedia.org/wiki/Map_projection#Azimuthal_(projections_onto_a_plane) has some information about the properties of an azimuthal equidistant projection.

Using a proj string we can construct such a projection centered on the point. Below is an example for Python 3.6+.

import pyproj
import json
from shapely.geometry import Point, mapping
from functools import partial
from shapely.ops import transform

point = Point(12, 34)

local_azimuthal_projection = f"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}"

wgs84_to_aeqd = partial(
    pyproj.transform,
    pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),
    pyproj.Proj(local_azimuthal_projection),
)

aeqd_to_wgs84 = partial(
    pyproj.transform,
    pyproj.Proj(local_azimuthal_projection),
    pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),
)

point_transformed = transform(wgs84_to_aeqd, point)

buffer = point_transformed.buffer(100_000)

buffer_wgs84 = transform(aeqd_to_wgs84, buffer)
print(json.dumps(mapping(buffer_wgs84)))

(EPSG codes did not work on my local system so I used the full proj string for EPSG:4326 instead.)