Tuesday, 17 March 2015

linestring - QGIS Extract Nodes with M-Values for Linear Referencing



I have a MultiLineStringZM layer in a sqlite database, and I'm trying to visualize the measures or m-values at the vertexes. I've tried looking up information on how to do this in QGIS, and about all I've been able to gather is that this is not possible directly from the linestring layer and that the points need to be extracted to a separate layer.


I've used Vector --> Geometry Tools --> Extract nodes to create a multipoint layer representing the vertexes of my multlinestring layer, but the process loses the m-values of the vertexes. I need the m-values preserved by either saving the m-value as an attribute of the point, or something else?


Internally we have a command line tool that converts linestrings to a point shapefile with the m-values stored as an attribute on each point, and I've used that to verify that there are m-values assigned to the vertexes, and I could use that if I have to, but if possible it would be nice if this could be done directly inside of QGIS.


EDIT - Repeating what I've said above, but re-emphasizing the fact that we do have a command line tool that can achieve the results that I'm looking for that uses the GDAL libraries, so a solution showing just a partial answer in PyQGIS is not the answer I'm looking for. I'm looking for a built in tool, a plugin ready made for QGIS, or a complete script that can extract (not create/generate) and visualize m-values from a MultiLineStringZM or LineStringZM geometry.



Answer



From what I can find there doesn't appear to be an existing solution for this exact situation, but I still wanted to be able to do this in QGIS, so I took the plunge into python scripting.


A guide for writing processing algorithms can be found here https://docs.qgis.org/2.18/en/docs/user_manual/processing/scripts.html


To use this code open up the Processing toolbox, then expand Scripts, then expand Tools. Select "Create new script" and copy and paste the code below into the script window (use caution when copying and pasting python code since whitespace is syntactically significant. If you are having problems put the code into a text editor that shows whitespace and make sure that it copied correctly). Save it wherever you want and there is an execute script button at the top of the window. After you save it you can "Add script from file" and permanently have the script under "User scripts".


When the processing window comes up select the layer that contains the vector geometry and select run. The script behaves the same way as "Extract Nodes" except that it adds a column called MValues and or ZValues depending on what is available in the input geometry.


##input_layer=vector

##output_layer=output vector

from qgis.core import QgsWKBTypes, QgsField, QgsVectorFileWriter, QgsFeature, QgsGeometry
from PyQt4.QtCore import QVariant

def addVertices( geometry, writer, inFeature ):
coordinateSequence = geometry.coordinateSequence()
for rings in coordinateSequence:
for points in rings:
for point in points:

feature = QgsFeature( fields )
feature.setGeometry( QgsGeometry( point ) )
type = point.wkbType()
attributes = inFeature.attributes()
if QgsWKBTypes.hasM( type ):
attributes.append( point.m() )
if QgsWKBTypes.hasZ( type ):
attributes.append(point.z())
feature.setAttributes( attributes )
writer.addFeature( feature )

return

inlayer = processing.getObject( input_layer )
provider = inlayer.dataProvider()
fields = provider.fields()
geomType = QgsWKBTypes.Type(inlayer.wkbType())
outputGeomType = QgsWKBTypes.Point

if QgsWKBTypes.hasM( geomType ):
outputGeomType = QgsWKBTypes.addM( outputGeomType )

fields.append( QgsField( "MValue", QVariant.Double ) )

if QgsWKBTypes.hasZ( geomType ):
outputGeomType = QgsWKBTypes.addZ( outputGeomType )
fields.append( QgsField( "ZValue", QVariant.Double ) )

layer_options = 'SHPT=' + QgsWKBTypes.displayString(outputGeomType)
writer = QgsVectorFileWriter( output_layer, 'UTF-8', fields, outputGeomType , inlayer.crs(), layerOptions=[layer_options] )

features = inlayer.getFeatures()

featureCount = inlayer.featureCount()
featureIndex = 0

for f in features:
percent = ( featureIndex/float( featureCount ) ) * 100
progress.setPercentage( percent )
g = f.geometry().geometry()
addVertices( g, writer, f )
featureIndex +=1


del writer

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