What's the best way to produce a density map from weighted points in QGIS?

Question

I'd really like to repeat the functionality of the ArcGIS Kernel Density tool - how could I do that?

I would like to limit answers to easy-to-install free-and-open-source software – i.e. QGIS installs easily with GRASS on all platforms so QGIS+GRASS would be fine, but SAGA does not (unfortunately, since it has what may be the ideal tool).

I am trying to produce density maps of wildlife across protected areas in QGIS. The density map in the following example was produced in ArcGIS from point observations of wildlife, with raw counts of group size (a field in the vector file) used to weight the density of each grid cell in the Kernel Density (SpatialAnalyst) tool, with a chosen search radius and grid cell size:

In a previous question on density, it was suggested to use the GRASS v.kernel tool to mimic the ArcGIS Kernel Density tool, but v.kernel doesn't do the same job. After looking at the manual and (successfully) producing a density map, it seems that the v.kernel function only works with point density, and there is no opportunity to give a variable for each point (such as raw counts) to weight each point.

UPDATE

There seem to be various *.surf.* tools in GRASS that may help create a density surface - and these do accept a weighting column or z-value, or are performed on rasters. @underdark suggested v.surf.rst - and the 'zcolumn' would be my weighting (count) variable - but I cannot figure out how to ask the tool to create a certain grid size or use a certain radius.

Suggestions on how to use v.surf.rst or any other method?

---

Sample data

x,y,count
431250,8707500,0
418750,8707500,5
413750,8707500,3
411250,8707500,1
408750,8685000,0
411250,8685000,0
416250,8685000,0
416250,8682500,6
411250,8682500,3
418750,8680000,0
433750,8677500,3
421250,8677500,0
423750,8675000,1
431250,8672500,0
428750,8672500,2
426250,8672500,2
423750,8670000,0

Answer 1

According to its manual page, the GRASS command r.resamp.filter will do for rasters representing point data exactly what ArcGIS will do for point layers: use the filter=box option for a "simple" raster and use the filter=gauss option for the other ArcGIS kernel. Use the -n flag to avoid propagating nulls.

---

Note that kernel density estimates (aka "heat maps") are not interpolations of the data. The value of a KDE at a location x estimates the amount of a value "Z" per unit area near x. (The radius or "bandwidth" quantifies what "near" means.) The values of Z do not need to be defined at every possible location on the map. For example, Z might represent the presence of something such as a person, in which case the KDE gives population density. Nor do the values of Z need to vary continuously across the map. For interpolation, it is assumed that Z is defined at all locations and that the data are observations of the values of Z at specified points. The interpolator attempts to predict the unobserved values of Z at all other points. This would make sense when Z is, say, a temperature or pressure, but is usually nonsensical when Z records the presence of something or when the data are a complete census. (In the latter case, contemplate what a road density map for a region might mean and how one could possibly make sense of "interpolating" roads across non-road areas.)

Answer 2

SAGA's module 'Kernel Density Estimation' is what you are looking for.

Install SAGA Module interface in QGIS (in menu: Plugins --> Fetch Python plugins..) and use the module. Good luck!

Answer 3

A really simplistic method with GRASS GIS that is closer to Point Density in ArcMap than Kernel Density:

1. Import the points to a raster map with r.in.xyz using method=sum at a specified raster resolution (set with g.region).

2. Use r.neighbors to smooth the map with method=average (which is default) and use the option size to set the search radius.

(I don't have access to GRASS at the moment so I didn't actually try this!)

Answer 4

As you requested more guidance on v.surf.rst, here are my inputs

First, about the grid size - you can use Plugin -> GRASS -> Edit Current GRASS region and set the output resolution. Your output from v.surf.rst will have that resolution.

For the radius, the 'tension' seems to be the parameter. I am no expert on this algorithm but reading from the manual, this seems to be the relevant bit

"... High tension "increases the distances between the points" and reduces the range of impact of each point, low tension "decreases the distance" and the points influence each other over longer range).."

So you could use the tension parameter roughly as you would use the radius parameter.

From your sample data, the result of v.surf.rst looks like below and seems reasonable given it is using the counts as weights for interpolation

Answer 5

While I haven't tested it, in the QGIS Contributors repository there is a plugin called 'Home range estimation with R'. This includes Kernel (density?) calculations. I reckon, if it works, that this would be the best method. R will do the real statistical method of calculating kernel density.

If you have R installed you should be able to just install the plugin and give it a go.

Answer 6

If you accept to do a bit of java programming outside of qgis, you could simply use this density map builder library.

Using the constructor HeatMapBuilder(int w, int h, int[][] pts, int[] weights), it is possible to give a weight for each point as you need. The output image can be retrieved with the getImage() method and saved on your disk with a ImageIO.write("mymap.png").

Here is an example of output:

It is possible to change the smoothing strength and the colour palette.