Skip to content

Instantly share code, notes, and snippets.

@johanez
Last active March 29, 2016 11:17
Show Gist options
  • Save johanez/ad11bb47bebe4d600c1f to your computer and use it in GitHub Desktop.
Save johanez/ad11bb47bebe4d600c1f to your computer and use it in GitHub Desktop.
Extract polygon nodes from vector layer and creates a point layer with four attribuites: the id of the origin polygon, the node numebr (per polygon), bearing to the follwing point, and distance to follow up node (if in projected CRS!)
##input_polygons=vector
##ref_x=number 5
##ref_y=number 4
##out_pnt_file=output vector
##out_ln_file=output vector
# Extract polygon nodes from polygons in a vector layer
# creates a point layer for the nodes, and a lien layer with
# with lies to the reference point from the closest node and the following node.
# Both out puts have four attributes:
# the id of the origin polygon,
# the node number (per polygon, starting at closest node to reference),
# bearing to the follwing point or reference (degeree), and
# distance to follow up node or reference (CRS units),
# if input has projected CRS.
# Author: Johannes Eberenz 2016
# License: GNU GPL
from qgis.core import *
from PyQt4.QtCore import *
from processing.tools.vector import VectorWriter
from math import sqrt
import sys
# import custom modules
sys.path.append('/home/jo/landmapp/dev/qgis/scripts')
import cert_features
# get layer
layer = processing.getObjectFromUri(input_polygons)
# prepare output
out_pnt_fields = [QgsField("poly_id", QVariant.Int),
QgsField("node_id", QVariant.Int),
QgsField("bear", QVariant.Double),
QgsField("dist", QVariant.Double),
QgsField("has_ref", QVariant.Int),
QgsField("bear_ref", QVariant.Double),
QgsField("dist_ref", QVariant.Double)]
out_pnt_writer = VectorWriter(out_pnt_file, None, out_pnt_fields,
QGis.WKBPoint, layer.crs())
out_ln_writer = VectorWriter(out_ln_file, None, out_pnt_fields,
QGis.WKBLineString, layer.crs())
# ref pint
ref_pnt = QgsPoint(ref_x, ref_y)
# check projection
progress.setInfo(layer.crs().authid())
crs_geo = layer.crs().geographicFlag()
if crs_geo:
progress.setInfo("Layer is in LatLong, distances will not be calculated!")
# iterate over polygons
iter = layer.getFeatures()
for feature in iter:
# retrieve every feature with its geometry and attributes
geom = feature.geometry()
if geom.type() == QGis.Polygon:
poly = geom.asPolygon()
# nodes list of first ring of poly, assuming there is only one!
nodes = poly[0][:-1]
print(nodes)
# calc dists to ref to select entry node
ref_dists = [sqrt(ref_pnt.sqrDist(nodes_pnt)) for nodes_pnt in nodes]
print ref_dists
#for i in range(len(poly[0]) - 1):
# ref_dist = sqrt(nodes[i].sqrDist(ref_pnt))
# reorder
index_entry = ref_dists.index(min(ref_dists))
nodes_o = nodes[index_entry:]
nodes_o.extend(nodes[:index_entry+1])
print(index_entry)
print(nodes_o)
for i, node in enumerate(nodes_o[:-1]):
print(i, node)
# get dist and bearing
dist = sqrt(node.sqrDist(nodes_o[i+1]))
azimuth = node.azimuth(nodes_o[i + 1])
# calc reference bearing and distance for last and first pnt
ref_dist = sqrt(node.sqrDist(ref_pnt))
ref_bear = node.azimuth(ref_pnt)
if azimuth < 0:
azimuth = 360 + azimuth
nodes_o_text = ("node==%i: x=%2d, y=%2d, azim=%6.1f, dist=%6.1f"
% (i, node.x(), node.y(), azimuth, dist))
print(nodes_o_text)
# write line features for entry and exit
if (i == 0 or i == (len(nodes_o) - 2)):
has_ref = 1
print ('line!: ', i)
line = QgsFeature()
ln = QgsGeometry.fromPolyline([ref_pnt, node])
line.setGeometry(ln)
line.setAttributes([feature.id(), i, azimuth, dist,
has_ref, ref_bear, ref_dist])
print([ref_pnt, node])
out_ln_writer.addFeature(line)
else:
has_ref = 0
# write node features
point = QgsFeature()
pnt = QgsGeometry.fromPoint(node)
point.setGeometry(pnt)
point.setAttributes([feature.id(), i, azimuth, dist,
has_ref, ref_bear, ref_dist])
out_pnt_writer.addFeature(point)
else:
progress.setInfo("Skipped non polygon feature!")
# write file
del out_pnt_writer
del out_ln_writer
@johanez
Copy link
Author

johanez commented Mar 29, 2016

Update:
adds lines to ref. point from closest nodes

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment