5

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
  • 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 '18 at 4:46
4

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

  • 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. – Simon Willison Jan 17 '18 at 19:39
  • 1
    Excellent. By the way, shapely's docs are now at shapely.readthedocs.io/en/stable. – sgillies Jan 17 '18 at 21:53
3

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.

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.