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:
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:
Here's a quick PyQGIS script which should do the trick
from qgis.core import QgsFeature, QgsVectorFileWriter, QgsGeometry
def create_points(feat,writer):
geometry = feat.constGeometry()
if not geometry:
return
length = geometry.length()
# -----------------
# change 'num_points' to match your field name for the number of points field
# -----------------
num_points = feat['num_points']
delta = length / ( num_points + 1.0 )
distance = 0.0
for i in range(num_points):
distance += delta
output_feature = QgsFeature(feat)
output_feature.setGeometry( geometry.interpolate(distance) )
writer.addFeature(output_feature)
layer = iface.activeLayer()
# ---------------
# change 'd:/test_points.shp' to desired output file name
# ---------------
writer = QgsVectorFileWriter('d:/test_points.shp',None, layer.fields(), QGis.WKBPoint, layer.crs())
for f in layer.getFeatures():
create_points(f,writer)
del writer
Just change the num_points
field name and output file name to match your data, select the input layer, and run it in the python console.
Best 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.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.That will create a new point layer with points at every 30 meters along the selected lines, like so:
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:and you can remove the