Last active
March 29, 2016 11:17
-
-
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!)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
##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 |
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
Update: For last and first points of each polygon, distance and bearing to a reference point is calculated.