enter image description hereI'm diving into a project in QGIS and have a dataset of GPS points uploaded in Guatemala. My goal is to create five 10x10km polygons strategically dropped in areas with the highest pin density.

Here's the plan:

Box Generation for Highest Density (Polygon 1): Drop the first 10x10 polygon where it captures the most pins. Once I have the correct placement, delete all the points within this polygon.

Box Generation for Second Highest Density (Polygon 2): Drop the second 10x10 polygon where it would have the most pins, excluding the ones deleted in polygon 1. Delete all the points within this polygon.

Repeat for Boxes 3, 4, and 5:

Continue the process until I have five polygons, each placed in areas with decreasing pin density.

I've tried to use heatmaps to help with this task, but it's still proving difficult.

How do I tackle this in QGIS?

The desired outcome even would be a centroid point for 10km radius of highest density points.

  • @BERA I have added an example. Thank you!
    – ALO
    Dec 20, 2023 at 20:04
  • You write "The desired outcome even would be a centroid point for 10km radius of highest density points." which is different from the previous stated goal of the 10 km box. Which one do you need? And: in your scenary are these boxes allowed to overlap?
    – Vincé
    Dec 22, 2023 at 9:46
  • 1
    How accurate do you need it to be? I would draw 10km rectangles around every point, then count the points within each rectangle, keep the one with the most points inside, take those points out of the original set and repeat
    – Andre Geo
    Dec 22, 2023 at 16:06
  • See here: gis.stackexchange.com/a/470268/88814
    – Babel
    Dec 23, 2023 at 11:04

2 Answers 2


This will probably not find the optimal solution but maybe good enough.

Create a model (or execute the tools manually): enter image description here

The model creates random points in the point layer extent, then rectangles around each point with width and height 10000, and join the points (summary) to the rectangles to count how many points there are in each rectangle.

I create 5000 random points. You can probably create more, the execution time for model and script is < 30 s.

Make sure to tick Discard any records which could not be joined

Then modify and execute the code below.

The script lists all random rectangles and sort the list by point count.

Then iterates over each rectangle from highest point count to lower, and check if each rectangle overlaps any of the already found rectangles with high point count. If not it will save it.

layer = QgsProject.instance().mapLayersByName("join_output")[0] #The name of the Join output. Adjust

#The field in the joined layer with number of points in the polygon.
count_field = "id_count" 

polygon_features = [f for f in layer.getFeatures()] #List all random polygons
polygon_features.sort(key=lambda x: x[count_field], reverse=True) #Sort from max to min point count

polygon_top_counts = [] #A list to hold the five polygons with the highest point count
polygon_top_counts.append(polygon_features.pop(0)) #Append the polygon with the hightest point count (index 0)

for polygon in polygon_features: #For each randomly placed polygon, sorted from max to min point count
    print(len(polygon_top_counts)) #Print the number of features in the list
    if not any([polygon.geometry().intersects(x.geometry()) for x in polygon_top_counts]):
    if len(polygon_top_counts)>4: #if there are 5 polygons in the list stop searching

#Create a layer of the five polygons
ids_with_highest_point_count = [f.id() for f in polygon_top_counts]
    {'INPUT':layer,'EXPRESSION': f"$id IN{tuple(ids_with_highest_point_count)}",

enter image description here


I would start with a model like this, that works for the first 10km block:

base processing model

However, to repeat the process five times, you would need to adapt the code of the model by adding a for loop (Export model as Script Algorithm):

from qgis.core import (QgsProcessing,
from qgis.PyQt.QtCore import QVariant
import processing

class YourModel(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFeatureSource('points_layer', 'Point Layer', types=[QgsProcessing.TypeVectorPoint], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Polygons', 'polygons', type=QgsProcessing.TypeVectorPolygon, createByDefault=True, defaultValue=None))
        #self.addParameter(QgsProcessingParameterFeatureSink('Nextpoints', 'nextPoints', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(100, model_feedback)
        results = {}
        outputs = {}

        source = self.parameterAsSource(parameters, 'points_layer', context)
        points_layer = parameters['points_layer']
        output_fields = source.fields()
        output_fields.append(QgsField('NUMPOINTS', QVariant.Int))
        (output_polygons, dest_id) = self.parameterAsSink(parameters, 'polygons', context, output_fields, QgsWkbTypes.Polygon, source.sourceCrs())
        for num in range(5): # get the 5 densest polygons
            if feedback.isCanceled():
            # Rectangles, ovals, diamonds
            alg_params = {
                'HEIGHT': 10000,
                'INPUT': points_layer,
                'ROTATION': 0,
                'SEGMENTS': 36,  # not used
                'SHAPE': 0,  # Rectangle
                'WIDTH': 10000,
                'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
            rectangles = processing.run('native:rectanglesovalsdiamonds', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

            # Count points in polygon
            alg_params = {
                'CLASSFIELD': '',
                'FIELD': 'NUMPOINTS',
                'POINTS': points_layer,
                'POLYGONS': rectangles['OUTPUT'],
                'WEIGHT': '',
                'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
            rectangles_with_pointcount = processing.run('native:countpointsinpolygon', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

            # Extract by expression
            alg_params = {
                'EXPRESSION': '"NUMPOINTS" =  maximum( "NUMPOINTS")',
                'INPUT': rectangles_with_pointcount['OUTPUT'],
                'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
            fullest_polygon = processing.run('native:extractbyexpression', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
            fullest_polygon_layer = QgsProcessingUtils.mapLayerFromString(fullest_polygon['OUTPUT'], context=context)
            # if more than one polygon has max count:
            for feat in fullest_polygon_layer.getFeatures():
                output_polygons.addFeature(feat, QgsFeatureSink.FastInsert)

            # Extract by location
            alg_params = {
                'INPUT': points_layer,
                'INTERSECT': fullest_polygon['OUTPUT'],
                'PREDICATE': [2],  # disjoint
                'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
            outputs['ExtractByLocation'] = processing.run('native:extractbylocation', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
            # set new point layer for next iteration:
            points_layer = outputs['ExtractByLocation']['OUTPUT']
        results['Polygons'] = output_polygons
        return results

    def name(self):
        return 'dense_blocks'

    def displayName(self):
        return 'Draw High-Density Blocks'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return YourModel()

Here the points inside the highest-density 10km block are excluded completely from the next iteration, meaning that the next highest-density block is calculated from the remaining points. This also means that the output blocks might overlap, but points will never be counted twice.

Execution time for 1300 points is less than 3 seconds. Since each block has a point in the middle, no blocks are created where no points are...

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.