Calculating the number of points per raster pixel

Question

I'm working with QGIS 3.14.16, where I have a large dataset of points (in a shapefile) and a raster of my study area.

I want to calculate a new raster with the number of points in each cell. I've done a bit of research, and saw some posts that suggested using the GRASS function r.in.xyz.

However, I cannot find that function in my toolbox. Does anyone have any suggestions about how I could perform this calculation?

And, as an extra, does someone know why I don't have the r.in.xyz tool?

Answer 1

you could try the function raster pixels to polygons. With the output of that, run count points in polygons where you end up with a poly dataset that has a field indicating the count of point that fall in each pixels.
Finally, rasterize the polygon dataset.

Answer 2

The reason that r.in.xyz is not in the QGIS Processing toolbox rests in the way the QGIS works with GRASS modules. Each module runs in a temporary GRASS Location/mapset, and these are setup based on the input date to the module. However the input to r.in.xyz - a table of x-y-z locations - does not specify the raster resolution, so there's no easy way for QGIS Processing to setup the region in the temporary Location.

May I suggest to run this natively in GRASS? Here's how to do it:

1. First you need to export your shapefile to a table with, at least, the X-Y coordinates. You can do this easily in QGIS: Right click on the point layer, Export->Save Feature as..., then in the top Format dropdown, choose Comma separated layer. Next in the Filename text box, click the ellipsis (at right end) to choose output location and filename. Below in the Geometry section, choose Point, and further down, under Layer Options select in the GEOMETRY dropdown the AS_XYZ option. Now click OK to export.
2. Next step is to start GRASS GIS and use your raster to determine the Location coordinate reference system, and region. When you first start GRASS, the Location/Mapset wizard opens. In the left Location panel, click New to start a new Location, then Next. Enter a directory for the GRASS database, and a name for your new Location, and again Next. Here select the second option Read projection from georeferenced file and browse to your raster and select it for the georeferenced file. After clicking Next you will be asked if you want to import that file. Select Yes. And finally, back to the wizard, and click Start GRASS .
3. Now you're in a running GRASS session, and you can prepare for the r.in.xyz command. First is to insure that the current computational region matches your raster. This is always the first step working with GRASS. In the Layer Manager, from the Menu choose Settings->Computational Region->Set region. Then in the g.region window, choose the second option Set region to match raster map and choose your raster map from the drop down. and click Run
4. Now, finally, for r.in.xyz, Select from the Menu File->Import raster->ASCII x,y,z point. In the r.in.xyz module window, click browse to get to the XYZ file (that you prepared in step 1 above), and enter a name for the output raster. Go next to the Statistic tab, and choose n (instead of the default mean) so you'll get a count of points in each output raster cell. Go on to the Input tab and switch to comma for the field separator, instead of pipe. Make sure that the column numbers Column for x, Column for y, Column for data match the actual columns in the XYZ file. 1,2,3 should be right. That's it, you should be good to go. Click run...
5. If you want to export the new raster, with point counts, to a Geotiff, then from the Menu File->Export raster->Common export formats will handle that.

Answer 3

I know you are asking of a Qgis or a grass answer, however this can be done easily in R by rasterizing your point shape file. You use the terra package to do that (would work with the raster package as well, but terra is more recent and better)

library(terra)
the_raster_from_your_study_area <- rast("path/to/your/raster/file")

your_point_shapefile <- vect("path/to/your/shapefile")

result <- rasterize(your_point_shapefile, the_raster_from_your_study_area, fun="length")

It is fast and easy.

I'm pretty sure you would get the same result by using the rasterize tool from gdal directly in QGIS, but I haven't tried.