Tuesday, 7 June 2016

Pyqgis Fixed Buffer Creation in Meters with Haversine


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:




  1. Get lat long of selected point feature.

  2. Use reverse haversine formula to get coordinates of new points 1km away from origin at specified angle intervals.

  3. Join new points to form buffer ring polygon feature.


  4. 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:




  1. 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".

  2. 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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...