# Creating select High-Density Polygons from GPS Points in QGIS

I'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? Dec 22, 2023 at 9:46
• 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 Dec 22, 2023 at 16:06
• See here: gis.stackexchange.com/a/470268/88814 Dec 23, 2023 at 11:04

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

Create a model (or execute the tools manually):

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]):
polygon_top_counts.append(polygon)

if len(polygon_top_counts)>4: #if there are 5 polygons in the list stop searching
break

#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)}",
'OUTPUT':'TEMPORARY_OUTPUT'})


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

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,
QgsProcessingUtils,
QgsProcessingMultiStepFeedback,
QgsFeatureSink,
QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink,
QgsWkbTypes,
QgsField)
from qgis.PyQt.QtCore import QVariant
import processing

class YourModel(QgsProcessingAlgorithm):

def initAlgorithm(self, config=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():
break
# 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():

# 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']

feedback.setCurrentStep((num+1)+20)

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