On Qgis 2.18, I'm trying to write a python code that allows me to create a 1km, 2km, 3km, 4km, 5km buffer around a selected point feature. To avoid having to reproject things, I'm using the reverse haversine formula. My workflow is as follows:
- Get lat long of selected point feature.
- Use reverse haversine formula to get coordinates of new points 1km away from origin at specified angle intervals.
- Join new points to form buffer ring polygon feature.
Repeat for 2km, 3km, 4km, 5km.
from math import *
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
layer = iface.activeLayer()
for feat in layer.selectedFeatures():
geom = feat.geometry()
start_x = geom.asPoint().x()
start_y = geom.asPoint().y()
start_lon = start_x * pi / 180
start_lat = start_y * pi / 180
sides = 64
radius = 6378137.0 # meters
points = []
distance = [1000, 2000, 3000, 4000, 5000]
for i in range(len(distance)):
dist = distance[i]
degrees = 0
while degrees < 360:
degrees = degrees + 360/sides
bearing = degrees * pi / 180
end_lat = asin((sin(start_lat) * cos(dist / radius))+ (cos(start_lat) * sin(dist / radius) * cos(bearing)))
end_lon = start_lon + atan2(sin(bearing) * sin(dist / radius) * cos(start_lat), cos(dist / radius) - (sin(start_lat) * sin(end_lat)))
points.append(QgsPoint(end_lon * 180 / pi, end_lat * 180 / pi))
geometry = QgsGeometry.fromPolygon([points])
Two places where I'm stuck:
- When I add this script on the Qgis processing toolbox (Add script from file) and try to run it, I get an error which says "iface is not defined".
- How can I save each iteration of the distance buffer as a shp file?
Answer
To return the xy
coordinates of a selected feature(s), you could use the following:
layer = iface.activeLayer()
for feat in layer.selectedFeatures():
geom = feat.geometry()
print geom.asPoint().x(), geom.asPoint().y()
No comments:
Post a Comment