[GIS] How to create equidistant points in QGIS along a polygon

pythonqgis

I am looking for a solution to distribute with a python script equidistant points in QGIS along a polygon. I tried to use the following code example:

http://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/

Using it in my example (running over each polygon from a polygon shape file):

# Reading a polygon layer
pp_layer = QgsVectorLayer("poly.shp" , "polylyr", "ogr")

# create layer to store result"
pntl = QgsVectorLayer("Point", "name", "memory") 
pr = pntl.dataProvider()  

for feature in pp_layer.selectedFeatures():
   geom = feature.geometry()

   currentdistance = self.distance
   feats = []    
   # geom.asPolyline()
   while currentdistance < length: 
      point = geom.interpolate(currentdistance)
      print point
      fet = QgsFeature()
      fet.setGeometry(point)
      feats.append(fet)
      currentdistance = currentdistance + self.distance
   pr.addFeatures(feats)

However, when running the code, if fails. I got an error message in the python console

Traceback (most recent call last):
File "<input>", line 1, in <module>
File "/home2/fun.py", line 92, in run_fun
fet.setGeometry(point)
TypeError: QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType'

an I got an errormessage from qgis

    ERROR 1: IllegalArgumentException: LinearIterator only supports lineal geometry components

Does anyone has an idea, if i have an error in my code or is it not possible to solve the question with this approach? I am using QGIS 2.0.1.

Best Answer

You cannot use interpolate() with polygons:

If you have a polygon:

geom.exportToGeoJSON()
u'{ "type": "Polygon", "coordinates": [ [ [167546.39129160347511061, 156935.60758662305306643], [167464.55340823522419669, 154998.77768024150282145], [169183.148958968144143, 154698.70544122465071268], [167901.02211953248479404, 153498.41648515724227764], [167109.92258030621451326, 152734.59624038706533611], [165364.04773511723033153, 154725.98473568074405193], [167546.39129160347511061, 156935.60758662305306643] ] ] }'  
geom.asPolygon()
[[(167546,156936), (167465,154999), (169183,154699), (167901,153498), (167110,152735), (165364,154726), (167546,156936)]]
# interpolate
geom.interpolate(100).asPoint()
Traceback (most recent call last):
File "<input>", line 1, in <module>
AttributeError: 'NoneType' object has no attribute 'asPoint'

The LinearRing of the Polygon is:

LinearRing = QgsGeometry.fromPolyline(geom.asPolygon()[0])
LinearRing.exportToGeoJSON()
u'{ "type": "LineString", "coordinates": [ [167546.39129160347511061, 156935.60758662305306643], [167464.55340823522419669, 154998.77768024150282145], [169183.148958968144143, 154698.70544122465071268], [167901.02211953248479404, 153498.41648515724227764], [167109.92258030621451326, 152734.59624038706533611], [165364.04773511723033153, 154725.98473568074405193], [167546.39129160347511061, 156935.60758662305306643] ] }'
# interpolate
LinearRing.interpolate(100).asPoint()
(167542,156836)

Polygons:
enter image description here

Part of a class to create memory layers from the Python console:

class Crea_layer(object):
   def __init__(self,nom,type):
        self.type=type
        self.name = nom
        self.layer =  QgsVectorLayer(self.type, self.name , "memory")
        self.pr =self.layer.dataProvider() 
   def create_point(self,geometry):
        self.seg = QgsFeature()
        self.seg.setGeometry(QgsGeometry.fromPoint(geometry))
        self.pr.addFeatures([self.seg])
        self.layer.updateExtents()
   @property
   def show_layer(self):
        QgsMapLayerRegistry.instance().addMapLayers([self.layer])

Creation of the equidistant points (step = 100m)

polys = qgis.utils.iface.activeLayer()
points = Crea_layer("equidistant", "Point")
for feat in polys.getFeatures():
    geom = feat.geometry()
    line = QgsGeometry.fromPolyline(geom.asPolygon()[0])
    for distance in xrange(0,line.length(),100):
         point = line.interpolate(distance)
         point = point.asPoint()
         point.create_point(QgsPoint(point.x(),point.y()))

 points.show_layer

Result:

enter image description here

Related Solutions

QGIS – Creating Equidistant Points in QGIS

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

[GIS] Creating set of equally spaced points along line using QChainage and exporting these as latitude and longitude coordinates in degrees

The lines to buffer should be contiguos.

I usually take the borders from Natural Earth world boundaries, which comes as polygons. The border lines from gadm.org should be closed as well.

I described a workflow to combine lines here: How to visually sort the order of Points with QGIS, GDAL, ...?

Vector -> Geoprocessing -> Dissolve using dissolve all lets you create a single buffer around all segments.

Best Answer

Related Solutions

QGIS – Creating Equidistant Points in QGIS

[GIS] Creating set of equally spaced points along line using QChainage and exporting these as latitude and longitude coordinates in degrees

Related Question