1

There are about a million questions on this website about counting points within polygons. There are also a few on counting points within a grid, but most of the answers use a count-points-in-polygons method.

The problem is that this method is stupidly inefficient with a rectilinear large grid if the grid is aligned with the projection coordinates. Point-in-polygon calculations are not cheap, and if the grid is very large, there are a lot of them - we're in O(n_points * n_grid_points) territory.

A much more efficient method is to simply round the X and Y values of each point to the nearest grid line (this is as easy as taking the coordinate index of the last grid coordinate smaller than the point X value, but could be even more efficient on large grids using a bisection search), and then count each unique set of values.

There doesn't seem to be a way to do this in QGIS/GRASS though, unless I am missing something. Is there a straight-forward way to do this efficient points-to-grid counting?

5
  • efficient would be to use geopandas for that, see example here: stackoverflow.com/questions/48097742/geopandas-point-in-polygon
    – eurojam
    Jul 12, 2022 at 6:35
  • With QGIS expression for example, you can compute a new field with $x and $y and some other expressions to create a grid code just from the points, so after you can do statistics on this field. Jul 12, 2022 at 7:34
  • You'd need a supporting binary structure to implement a bisection mapping, while ordinate truncation is merely applying a set of bitwise operations for multiplication by a constant factor (you may want to work in unit integer space). In e.g. Python, loop over your points and increment i.e. dict[ (<rounded_x>, <rounded_y>) ] += 1 - then, for very large grids, maybe build a grid layer from the resulting map rather than joining to an existing layer.
    – geozelot
    Jul 12, 2022 at 7:52
  • @J.Monticolo: that sounds great. Ultimately I would want to get this into a raster. An example of what you mean would make a great answer
    – naught101
    Jul 12, 2022 at 12:34
  • This has a good method using geopandas: james-brennan.github.io/posts/fast_gridding_geopandas
    – naught101
    Jul 21, 2022 at 8:11

1 Answer 1

2

I do this in GRASS GIS as follows. (Use case: to examine how many Lidar pointcloud points are in each grid cell). I save the vector points as a CSV file, with X and Y coordinates. This can be done, of course, in GRASS, or QGIS, whatever is comfortable. Then in GRASS, I set the computational region to the required grid, and the resolution to the required size of the grid cells

So the GRASS "computational region" represents, (and replaces) your vector grid. It is defined using the g.region module, with the extents set to extents of your original vector grid, and the resolution to the size of your original grid cells. For example, if your vector grid was created with 100m spacing, then:

g.region -ap vector=<original vector grid> resolution=100

Now I run r.in.xyz pointing to the CSV file for input. I use the method=n parameter, such that each cell in the output raster obtains the count of points in that cell.That raster can be polygonized, if necessary, etc....

2
  • How do you define your grid?
    – naught101
    Jul 14, 2022 at 3:37
  • See edits to the answer...
    – Micha
    Jul 14, 2022 at 8:56

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.