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;
Best Answer
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.
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.