Creating buffer circle x kilometers from point using Python?
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.
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)]
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: