# Counting points in a rectilinear grid efficiently

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?

• efficient would be to use geopandas for that, see example here: stackoverflow.com/questions/48097742/geopandas-point-in-polygon 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. 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 Jul 12, 2022 at 12:34
• This has a good method using geopandas: james-brennan.github.io/posts/fast_gridding_geopandas Jul 21, 2022 at 8:11

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