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

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.

• Given the answer I've provided below, can I get this question re-instated? I think my answer will be very useful to others, but I'd also like to get it reviewed. – Simon Willison Jan 17 '18 at 1:04
• Placing a question on hold does not prevent any of those actions but since it has been improved I have re-opened it so that more answers may be added. – PolyGeo Jan 17 '18 at 2:49
• Twenty points is far too few for a 100km radius polygon. Your function prototype is a very different shape, with a 1km radius. – Vince Jan 17 '18 at 4:46

## 2 Answers

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.)

• That's brilliant, thank you very much. I had to modify the code a tiny bit to add some missing imports but it's doing exactly what I needed now. – Simon Willison Jan 17 '18 at 19:39
• Excellent. By the way, shapely's docs are now at shapely.readthedocs.io/en/stable. – sgillies Jan 17 '18 at 21:53
• Thanks - this is very helpful. – Connor Dibble Jan 31 '20 at 4:12
• Does somebody know how to then clip gml shapefile with roads (linestrings in geodataframe) with such circle? – Vilq Feb 25 '20 at 15:41
• @bugmenot123 Because it seems that this code generates data invalid to current GeoJSON standard? It is quite old answer, so maybe it used to be correct. – reducing activity Mar 9 at 14:09

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.