Creating buffer circle x kilometers from point using Python?

Question

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
    Radian = angle * 0.017453 
    lat = math.asin(math.sin(latDD)*math.cos(d)+math.cos(latDD)*math.sin(d)*math.cos(Radian))
    lon=math.fmod(lonDD - math.asin(math.sin(Radian) * math.sin(d)/math.cos(lat))+math.pi,2*math.pi)-math.pi

    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:

Answer 1

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

Answer 2

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

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

it can be corroborated that buffer was properly produced:

Answer 3

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

def func_radius_around_point(point_list,m_list,crs_from,crs_to):

    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)]
 buffers_df = func_radius_around_point(warehouse_locs,
 [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])
 map_.add_child(poly)
 poly = folium.GeoJson(buffers_df.buffered_150000[0])
 map_.add_child(poly)

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.

The above code should give you the following output

Answer 4

the non-round shape is probably a function of the projection you use and the latitude. A circle in wgs coordinates at the equator wi