Construir el mapa de densidad para la gran base de datos

Question

Tengo un postgre/postgis base de datos con alrededor de 120 millones de puestos de trabajo (tipo geography). Me gustaría crear un mapa de densidad (por ejemplo, un mapa de calor) y servir esto a través de geoserver.

Me hizo crear una banda de ráster. Yo bucle a través de todos los puntos y aumentar el valor de la trama elemento que contiene este punto. (véase el código de abajo) Esto lleva a mi base de datos de la explosión (querery corrió alrededor de 15 horas, de la base de datos aumentaron un 60 gigabye).

Así que necesito más rápido y más eficiente manera de hacer eso.

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();

Cuando esto se hace, me gustaría crear el mapa de colores como

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

Hice un poco de prueba y llegó a la siguiente solución basada en esta entrada del blog de Vipin Raj. Creo que esta es la mejor respuesta hasta ahora, ya que todo se hace en el postgis servidor.

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;

Las funciones que se basa una consulta y actualizaciones de varios píxeles de la trama en el mismo tiempo. Mi punto de referencia mostró una tasa de actualización (o a la inversa) de 0,69 ms por píxel. Eso significa, que mi ráster de alemania podría tomar alrededor de 3 horas, que es justo lo suficiente para mí.