Friday, 24 March 2017

Creating equidistant points in QGIS?


I am trying to create points (new layer) at specific distance along the road (existing layer) in QGIS. Creating Regular points every meter at county level using ArcGIS Desktop? gives solution for ArcGIS. How to achieve this in QGIS? Adding points to point vector layer using QGIS? explains how to create points but does not do anything about the distance.





(I applied the proposed solutions with different measures of lengths for I did not know the conversion) @Nathans's solution worked to some extent, I got...


enter image description here. Here, the projection of these equidistant points is different from the original line.


With @underdark's suggestion, I got


this image where the points does not seem to be equidistant. I guess there is some projection issue with both of these which I am not understanding.



Answer



Note: There is now a QGIS plugin QChainage. It does all this and more. The code below is out of date with QGIS 2.0 and above.


Here is some Python code that you can stick in a file and use inside QGIS:


QGIS does have a method in it API to do liner referencing however I couldn't get it to work correctly, but I will contact the author of the code and see if I was doing something wrong.


For now you will need the shapely Python library, which you should install anyway because it's handy to have around. It also has great documentation at http://toblerity.github.com/shapely/manual.html


This is the section I am using in the following example http://toblerity.github.com/shapely/manual.html#interoperation.



Most of the following code is QGIS boilerplate code just creating the features, layers, converting from wkb and wkt and back. The core bit is the point = line.interpolate(currentdistance) which returns a point at a distance along a line. We just wrap this in a loop until we run out of line.


import qgis
from qgis.core import *
from PyQt4.QtCore import QVariant
from shapely.wkb import loads
from shapely.wkt import dumps

vl = None
pr = None


def createPointsAt(distance, geom):
if distance > geom.length():
print "No Way Man!"
return

length = geom.length()
currentdistance = distance
feats = []

while currentdistance < length:

line = loads(geom.asWkb())
point = line.interpolate(currentdistance)
fet = QgsFeature()
fet.setAttributeMap( { 0 : currentdistance } )
qgsgeom = QgsGeometry.fromWkt(dumps(point))
fet.setGeometry(qgsgeom)
feats.append(fet)
currentdistance = currentdistance + distance

pr.addFeatures(feats)

vl.updateExtents()

def pointsAlongLine(distance):
global vl
vl = QgsVectorLayer("Point", "distance nodes", "memory")
global pr
pr = vl.dataProvider()
pr.addAttributes( [ QgsField("distance", QVariant.Int) ] )
layer = qgis.utils.iface.mapCanvas().currentLayer()
for feature in layer.selectedFeatures():

geom = feature.geometry()
createPointsAt(distance, geom)

QgsMapLayerRegistry.instance().addMapLayer(vl)

Copy and paste the above code into file, I called my locate.py, in ~./qgis/python directory (because it is in the Python path) and just do the this in the Python console inside QGIS.


 import locate
locate.pointsAlongLine(30)

That will create a new point layer with points at every 30 meters along the selected lines, like so:



enter image description here


Note: Code is pretty rough and might need some clean up.


EDIT: The lastest QGIS dev build can now do this natively.


Change the while loop in createPointsAt to:


 while currentdistance < length: 
point = geom.interpolate(distance)
fet = QgsFeature()
fet.setAttributeMap( { 0 : currentdistance } )
fet.setGeometry(point)
feats.append(fet)

currentdistance = currentdistance + distance

and you can remove the


from shapely.wkb import loads
from shapely.wkt import dumps

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