14

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.

3
  • 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. Jan 17, 2018 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, 2018 at 2:49
  • 1
    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, 2018 at 4:46

2 Answers 2

10

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.

1
  • This threw an error for some reason and I had to change the code slightly such that in line 5: p = [-90.0667, 29.9500] Jan 14 at 8:05
8

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

8
  • 1
    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. Jan 17, 2018 at 19:39
  • 1
    Excellent. By the way, shapely's docs are now at shapely.readthedocs.io/en/stable.
    – sgillies
    Jan 17, 2018 at 21:53
  • Thanks - this is very helpful. Jan 31, 2020 at 4:12
  • Does somebody know how to then clip gml shapefile with roads (linestrings in geodataframe) with such circle?
    – Vilq
    Feb 25, 2020 at 15:41
  • 1
    @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. Mar 9, 2021 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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