Build Density Map for large dataset

Question

I've got a postgre/postgis database with about 120 million positions (geography type). I would like to create a density map (e.g. a heatmap) and serve this via geoserver.

I did create a one-band raster. I loop through all points and increase the value of the raster element containing this point. (see code below) This leads to my database exploding (querery ran about 15 hours, database increased by 60 gigabye).

So I need a faster and more efficient way of doing that.

CREATE TABLE raster_table (
  rid         SERIAL PRIMARY KEY,
  name        VARCHAR(30),
  description VARCHAR(255),
  rast        RASTER
);

/* insert new raster for whole germany */
DELETE FROM raster_table
WHERE name = 'test';
INSERT INTO raster_table (rid, name, description, rast)
VALUES (1, 'count', 'count raster',
        ST_AddBand(ST_MakeEmptyRaster(2500, 4500, 5.5, 47, 0.004, 0.002, 0, 0, 4326), '8BUI' :: TEXT, 0));

CREATE OR REPLACE
FUNCTION update_grid()
  RETURNS VOID AS $$
DECLARE
  pl  RECORD;
  val INT;
BEGIN
  FOR pl IN SELECT positions.id, positions.position
            FROM positions LOOP

    SELECT ST_Value(rast, 1, pl.position :: GEOMETRY)
    INTO val
    FROM raster_table;
    val = val + 1;

    UPDATE raster_table
    SET rast = ST_SetValue(rast, 1, pl.position :: GEOMETRY, val)
    WHERE raster_table.rid = 1;
  END LOOP;
END;
$$ LANGUAGE plpgsql;

SELECT update_grid();

When this is done, I would create the colormap like so

INSERT INTO raster_table (rid, name, description, rast)
  WITH ref AS (SELECT raster_table.rast
               FROM raster_table
               WHERE rid = 1)
  SELECT 2, 'rgba', 'colorful map', ST_ColorMap(ref.rast, 'bluered') FROM ref;

Answer 1

I did a bit of testing and came up with the following solution based on this blog post from Vipin Raj. I think this is the best answer until now, as everything is done on the postgis server.

CREATE OR REPLACE FUNCTION density_raster(_rid INTEGER)
  RETURNS VOID AS
$BODY$
DECLARE
  n          INTEGER = 128;
  query_part VARCHAR = '';
  query      VARCHAR = '';
  metadata   RECORD;
  i          INTEGER = 0;
  j          INTEGER = 0;
  point_x    DOUBLE PRECISION;
  point_y    DOUBLE PRECISION;
  value      DOUBLE PRECISION;
BEGIN
  /* get the metadata of the raster*/
  SELECT (F1.md).*
  INTO metadata
  FROM (SELECT ST_MetaData(rasters.rast) AS md
        FROM rasters
        WHERE rasters.rid = _rid) AS F1;
  RAISE NOTICE 'Found Raster from (%, %) to (%, %) with % elements', metadata.upperleftx, metadata.upperlefty,
  metadata.upperleftx + metadata.width * metadata.scalex,
  metadata.upperlefty + metadata.height * metadata.scaley,
  metadata.width * metadata.height;


  RAISE NOTICE 'Raster from has a width of % and a height of %', metadata.width, metadata.height;

  /* loop through all elements of the raster */
  WHILE (i < metadata.width) LOOP
    WHILE (j < metadata.height) LOOP
      point_x = metadata.upperleftx + i * metadata.scalex + 0.5 * metadata.scalex;
      point_y = metadata.upperlefty + j * metadata.scaley + 0.5 * metadata.scaley;

      SELECT count(fixes.*)
      INTO value
      FROM fixes
      WHERE ST_INTERSECTS(fixes.position, ST_MakeEnvelope(metadata.upperleftx + i * metadata.scalex,
                                                          metadata.upperlefty + j * metadata.scaley,
                                                          metadata.upperleftx + (i + 1) * metadata.scalex,
                                                          metadata.upperlefty + (j + 1) * metadata.scaley,
                                                          4326));

      query_part =
      query_part || ',' || 'ROW(st_setsrid(''POINT(' || point_x || ' ' || point_y || ')''::geometry, 4326), ' || value
      || ')';

      j = j + 1;
      n = n + 1;
      IF n > N
      THEN
        query_part = substring(query_part FROM 2 FOR (length(query_part) - 1));
        query =
        'UPDATE rasters SET rast = ST_SetValues(rast, 1, ARRAY[' || query_part || ']::geomval[]) WHERE rid = ' || _rid
        ||
        '';
        EXECUTE query;
        query_part = '';
        query = '';
        n = 0;
      END IF;
    END LOOP;
    j = 0;
    i = i + 1;


  END LOOP;
END;
$BODY$
LANGUAGE plpgsql;

The functions builds a query and updates several pixels of the raster at the same time. My benchmark showed an update rate (or the inverse of it) of 0,69ms per pixel. That means, that my raster of germany would take about 3 hours, which is fair enough for me.