2

I am trying to find a method to, as efficiently as possible, fill irregular shape polygons with points (defined spacing) or alternatively, polygons of defined size ( which i can find the point centroid).

I am using QGis, Python

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Thomas
  • 127
  • 2
  • 8
  • Have you considered QGIS Geoalgorithms - Vector creation tools - Random points inside polygons? – user30184 Nov 02 '16 at 15:27
  • 1
    What about creating a grid and creating points at its nodes? – Victor Nov 02 '16 at 15:48
  • Arcgis solution http://gis.stackexchange.com/questions/185889/arcmap-fill-polygon-with-points-highest-possible-number-bin-packing-problem/186025#186025 – FelixIP Nov 02 '16 at 19:50

1 Answers1

3

Next PyQGIS code uses a defined spacing of 50 m for generating points inside each feature of polygon layers.

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures()]

points = []

spacing_y = 50
spacing_x = 50

for feat in feats:
    extent = feat.geometry().boundingBox()
    xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

    rows = int(ymax - ymin)/spacing_y
    cols = int(xmax - xmin)/spacing_x

    x = xmin
    y = ymax

    geom_feat = feat.geometry()

    for i in range(rows+1):
        for j in range(cols+1):
            pt = QgsPoint(x,y)
            tmp_pt = QgsGeometry.fromPoint(pt)
            if tmp_pt.within(geom_feat):
                points.append(tmp_pt.asPoint())
            x += spacing_x
        x = xmin
        y -= spacing_y

epsg = layer.crs().postgisSrid()

#points
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running the code at the Python Console of QGIS, with a two features polygon layer, I got:

enter image description here

xunilk
  • 29,891
  • 4
  • 41
  • 80