# Creating buffer circle x kilometers from point using Python?

I am wanting to create a circle with an inputted radius around a point of predetermined latitude and longitude. I've been trying to use the formula found in this post which uses an algorithm documented by Ed Williams found here - but the outcome always creates a rectangle (seen below after plotted on google earth).

My code looks like:

``````bufferRange = input("enter radius of buffer zone")
bufferRange = float(bufferRange)
angle = 0
newCoords = []

for i in range(360):

lonDD = abs(lonDD)
d = bufferRange/6371

if 0 < angle < 90:
lat = latDD + lat
lon = lonDD + lon
elif 90 < angle < 180:
lat = latDD - lat
lon = lonDD + lon
elif 180 < angle < 270:
lat = latDD - lat
lon = lonDD - lon
elif 270 < angle < 360:
lat = latDD + lat
lon = lonDD - lon
elif angle == 90:
lat = latDD
lon = lonDD + lon
elif angle == 180:
lat = latDD - lat
lon = lonDD
elif angle == 270:
lat = latDD
lon = lonDD - lon
else:
lat = latDD + lat
lon = lonDD

lon = -lon

if angle < 360:
newCoords.append((lon, lat), )
else:
newCoords.append((lon, lat))

angle = angle + 1
``````

Where latDD and lonDD are the points at the center of the circle (Ottawa for this test).

Update:

After using xunilk's code (see below) the results improved, but the buffer is still not circular. I have a feeling its how I converted my longitude to meters, but I can't seem to find a better way to convert them online.

The following test was done with a buffer radius of 100km. The distance from the centre to the 'top' and 'bottom' of the buffer is exactly 100km. The distance from the centre to the sides of the buffer is 70.5km

``````polygonSides = 360
lat = (latDD*111.320)*1000
lonNum = 111320*math.cos(latDD)
lon = abs(lonDD)*111320

points_list = [ (-1*((lon + np.sin(angle)*bufferRange)/111320), (lat + np.cos(angle)*bufferRange)/111320)
for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ]
``````

Image of my new buffer:

• Where do use the value of `i`? Commented Jul 10, 2018 at 21:30
• I don't. I suppose I could use it by saying 'angle = i' and get the same result. I just want to iterate through the for loop 360 times (once for each degree used to create the buffer circle). Commented Jul 11, 2018 at 4:29
• if you display a 100 km radius circle in geographic coordinates, it will look like an oval. The reverse is true. Think about the purpose of your buffer, then select your type of distance and the projected coordinate system that suits you Commented Jul 19, 2018 at 7:40
• Do you need clear Python solution or you are free to use 3rd party libraries? Commented Jul 20, 2018 at 4:11
• How do you plot the points on google earth? I ask because the way you calculate the coordinates looks right, it will give just a bit different result to azimutal equidistant projection usage (because it uses the WGS84 ellipsoid model by default in PROJ), but your way is far more simple. Commented Jul 23, 2018 at 0:31

Use a spatial projection library to do the hard work. Adapting from a previous answer, use a dynamic azimuthal equidistant projection to do a geodesic buffer.

### Updated solution with PyProj 2.1+

``````from pyproj import CRS, Transformer
from shapely.geometry import Point
from shapely.ops import transform

def geodesic_point_buffer(lat, lon, km):
# Azimuthal equidistant projection
aeqd_proj = CRS.from_proj4(
f"+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0")
tfmr = Transformer.from_proj(aeqd_proj, aeqd_proj.geodetic_crs)
buf = Point(0, 0).buffer(km * 1000)  # distance in metres
return transform(tfmr.transform, buf).exterior.coords[:]
``````

### Original solution with PyProj 1.x

``````from functools import partial
import pyproj
from shapely.ops import transform
from shapely.geometry import Point

proj_wgs84 = pyproj.Proj('+proj=longlat +datum=WGS84')

def geodesic_point_buffer(lat, lon, km):
# Azimuthal equidistant projection
aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
project = partial(
pyproj.transform,
pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
proj_wgs84)
buf = Point(0, 0).buffer(km * 1000)  # distance in metres
return transform(project, buf).exterior.coords[:]
``````

## Example

``````b = geodesic_point_buffer(45.4, -75.7, 100.0)

print(b)
# [(-74.42290765358695, 45.39286001598599),
#  (-74.43102886629593, 45.304749544147974),
#  ...
# (-74.42290765358695, 45.392860015985995),
# (-74.42290765358695, 45.39286001598599)]
``````

• Thanks for the response Mike. I got a slightly different result than you with the same code (for some reason). When b printed for me the result was ((-74.42290765358695 45.39286001598599, -74.43102886629593 45.30474954414797 ... -74.42290765358695 45.392860015986, -74.42290765358695 45.39286001598599)). Additionally how did you plot the object into Google Earth? The simpleKML library wont accept b as an input with the typeerror: 'Polygon' object is not iterable. Commented Jul 19, 2018 at 15:16
• @DanielMutton I forgot to dump the coords at the end of the function; fixed now. It looks like your result should be the same now for my example lat/lon. I cobbled together a kml file in a non-elegant way, but simplekml should now understand the list of coordinate pairs. Commented Jul 19, 2018 at 19:58
• Sorry for the delayed response but this worked perfectly! Thanks so much Mike Commented Jul 24, 2018 at 17:25
• IDK, you probably saved me like 40 hours. Thanks for that! Commented Jun 9, 2021 at 23:29
• @Chang yes, it's possible, although challenging to get the maths correct. Libraries like GeographicLib can be used to use geodesics for distance, where it only depends on an ellipsoid (e.g. WGS84), and not a projection. Commented May 29, 2022 at 9:35

By using linspace method, from numpy python module, you can use following more concise code:

``````import numpy as np

bufferLength = 100  # 0.1 km
polygonSides = 360

x = 915884
y = 5042490

angles = np.linspace(0, 2 * np.pi, polygonSides, endpoint=False)
points_list = [(x + np.sin(a) * bufferLength,
y + np.cos(a) * bufferLength)
for a in angles]

print(points_list)
``````

where x, y represents an arbitrary point in Ottawa (26917 EPSG code; NAD83/UTM zone 17N)

By using following PyQGIS code (with only 50 points):

``````import numpy as np

bufferLength = 100
polygonSides = 50

layer = qgis.utils.iface.activeLayer()

points = [feat.geometry().asPoint() for feat in layer.getFeatures()]

epsg = layer.crs().postgisSrid()

angles = np.linspace(0, 2 * np.pi, polygonSides, endpoint=False)
buffer_points = [(points[0][0] + np.sin(a) * bufferLength,
points[0][1] + np.cos(a) * bufferLength)
for a in angles]

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
'buffer_points',
'memory')

prov = mem_layer.dataProvider()

feats = [QgsFeature() for i in range(len(buffer_points))]

for i, feat in enumerate(feats):
feat.setAttributes([i])
feat.setGeometry(QgsGeometry.fromPoint(
QgsPoint(buffer_points[i][0], buffer_points[i][1])
))

``````

it can be corroborated that buffer was properly produced:

• I'm getting a better result than I was before but my zone still isn't perfectly round, it's oval with the desired radius (100km for the test) at the Y axis, and a 70km radius along the X axis - similar in shape to a 0. This squishing becomes more severe the further north I test. Commented Jul 11, 2018 at 14:20
• Sounds like a coordinate system issue. What CRC are you using, have you tried others? Commented Jul 12, 2018 at 0:01
• Thanks for the response Keagan. I'm plotting these points in Google Earth using the simpleKML library in python. As far as I was aware Google Earth didn't use a coordinate system, but I will look into it more. Commented Jul 12, 2018 at 13:30

I might have a simpler albeit hacky solution that utilizes geopandas' transformation capabilities.

``````import pandas as pd
import geopandas as gpd
import folium

from shapely import geometry

points_df = pd.DataFrame(point_list).reset_index()
points_df = points_df.rename(columns={0:'geometry'})
points_gpd = gpd.GeoDataFrame(points_df,geometry='geometry',crs=crs_from)

for m in m_list:
buffered= points_gpd.to_crs(epsg=crs_to).buffer(m).to_crs(epsg=crs_from)
col_name = 'buffered_{m}'.format(m=m)
points_gpd[col_name] = buffered

return points_gpd

warehouse_loc = (74.175929,32.134755) #x,y
warehouse_locs = [geometry.Point(warehouse_loc)]
[100000,150000,200000],crs_from=4326,crs_to=3857)

map_ = folium.Map(location=(warehouse_loc[-1],warehouse_loc[0]))
poly = folium.GeoJson(buffers_df.buffered_100000[0])
The function above takes a list of `shapely.geometry.Point` objects in the `point_list` parameter, and a list of radii `m_list` and the respective CRS inputs to return a dataframe where every radius buffer is a new column.