Generating polygon representing rough 100km circle around latitude/longitude point using Python?
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.)
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.