12

I am using shapely and Python.

How should I go about creating a circle using shapely whose radius is in meters?

If i use

sampleCircle = Point(1,1).buffer(1)
print(sampleCircle.area)

This returns an area of 3.13654849055. How do I interpret the unit of this area? I read that this can be interpreted in square degrees but seek a little more explanation.

Also, how do I create a circle with a specified radius (in meters) around a point (defined using lat,long)?

The closest article I found was this http://comments.gmane.org/gmane.comp.python.gis/1582 but it lacks any concrete solution

2 Answers 2

18

By design Shapely is unaware of coordinate reference systems or units. To use it to solve real world problems, you must learn to transform longitude and latitude to an approximately planar local reference system (like a US State Plane) and you must keep track of your own units. I use pyproj to do this, but you can use whatever you want. Once you've transformed your long and lat degrees to x and y meters, Point(x, y).buffer(1.0).area is the area in m^2 of your circle (64-sided polygon to be precise).

6

Based on the link you provided (almost all work is done there), you can get a polygon representing that circle (64-side polygon).

from functools import partial

import pyproj
from shapely import geometry
from shapely.geometry import Point
from shapely.ops import transform

lon, lat = -122.431297, 37.773972 # lon lat for San Francisco
radius = 10000  # in m

local_azimuthal_projection = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(
    lat, lon
)
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"),
)

center = Point(float(lon), float(lat))
point_transformed = transform(wgs84_to_aeqd, center)
buffer = point_transformed.buffer(radius)
# Get the polygon with lat lon coordinates
circle_poly = transform(aeqd_to_wgs84, buffer)

2
  • I'm trying to understand the resolution argument in the .buffer()-method. Resolution defaults to 16, and this creates a 64 sided polygon. Why?
    – Olsgaard
    Commented Jun 15, 2022 at 9:37
  • The +R=6371000 parameter in the local azimuthal projection gave me incoherent results when combined with a geodesic distance function. Omitting the parameter makes the projection default to "GRS-80" and gave me expected results. Note: I did not plot near San Francisco. Commented Mar 20, 2023 at 19:34

Your Answer

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

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