4

beginner here I am trying to analyse a large set of points (2 millions) over a set of polygons (3000), and particularly find clusters of points. Doing it in QGIS was taking too long/crashing the compputer so I uploaded my layers on a Spatialite DB and use the interface provided by QGIS.

I have managed to do what I want using a few scripts in a row, but the cluster analysis is based on measuring the area of the dissolved cluster of buffered points, which doesn't take into account density. I have tried to add "weight" to the points and do a count but I can't get the query to work because the query results in a multi-part polygon.

Any idea to simplify or improve this process to a lower number of steps?

Query to calculate the number of points within the polygons:

SELECT b.name AS NameID,
count(a.geom) AS density
FROM 
Point_layer AS a, 
Polygon_layer AS b
WHERE ST_Intersects(a.geom, b.geom) = 1
group by NameID

Query to get the ID of points that intersects polygons (used to reduce the size of the point database to query:

 SELECT 
b.Name AS Name_id,
a.PointID AS Point_ID
FROM 
Point_layer AS a, 
Polygon_layer AS b
WHERE ST_Intersects(a.geom, b.geom) = 1
group by Point_ID

Then I sort out the points I need to a new layer, add a buffer and dissolve:

Select a.id, (st_unaryunion(st_collect(ST_Buffer(a.geom, 7)))) as geometry
FROM point_layer_filtered as a

Clusters

Then I split the resulting polygon into multipart into singlepart, and calculate area.

What I would like is to use something else than ST_unaryunion that allow to retain the number of points per clusters..

Is there a more efficient strategy that doesn't require buffers? Something that just group points together? I guess I could rerun the first script on the dissolved buffers but that would take forever

EDIT, Solution:

Response that might help others in the future.. The easiest and fastest way I found was to switch to PostGIS which has been orders of magnitudes faster and way easier to handle. It didn't end up being all that hard. If you're on a Windows machine, just install PostgreSQL with the PostGIS extension. You'll operate it via PGAdmin. To import geo files to your PostGIS databases you will need the PostGIS importer/exporter which is basically just a window that helps you connect to your Postgres database and upload file to it. There's good youtube videos walking you through the install.Depending what you do it might involves going back in forth between QGIS and PGAdmin. Adding PostGIS layer in QGIS is easy, click on the Elephant icon below the vector and raster "add layer" icon. Then I basically ran two queries back to back, (running in seconds instead of hours!). The first one buffers each points and dissolve to create clusters.

Create table buffered_cluster_20m (gid serial, geom geometry, area numeric, numID text);
Insert into buffered_cluster_20m(geom, area, numID)
Select (st_dump(st_union(st_buffer(a.geom, 7)))).geom as geometry,
st_area((st_dump(st_union(st_buffer(a.geom, 7)))).geom) as area, 
b.numID as numID
From point_samples as a, "donut_20m_cleaned" as b
Where ST_Intersects(a.geom, b.geom)
Group by numID

Problem is, area of cluster is not enough for intended purpose to measure the importance of a cluster (i.e.: 5 points could cover more area than 100 if they are farther apart than each others). So I run a second query that counts the points within the freshly generated areas to get a density instead of just a count:

Create Density_Cluster_20m (bufferID numeric, density numeric);
Insert into Density_Cluster_20m(density, bufferID)
SELECT 
(count(a.geom))/(b.area) as density,
count(a.geom) AS counts,
b.gid as bufferID,
b.numid as numID,
b.geom as geom 
FROM 
"sample_points" AS a, 
"buffered_cluster_20m" AS b
WHERE ST_Intersects(a.geom, b.geom)
group by bufferID, numID, b.area, b.geom
order by density desc

And that does it. I end up with a new layer that has areas with a density column that I can import into QGIS and browse through area of interests.

  • At 2m points, your spatialite query would probably go faster if you used an index. Indexing in spatialite is not automatic - see the manual. – BradHards Jun 14 '17 at 3:47
  • Thanks for your feedback, I didnt know that, will try next time. – Leo Jun 14 '17 at 5:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.