[postgis-users] Simple Point Density Surface

[Pierre Racine]

Use ST_MapAlgebraFct:

http://postgis.refractions.net/documentation/manual-svn/RT_ST_MapAlgebraFct.html

Create a function that take all the georeference parameters of the raster (width, height, ulx, uly, scalex, scaley, skewx, skewy, srid), the schema, the table and the name of the geometry column of the point table and a radius parameter, much like the following example, but returning the number of intersecting geometries instead:

-- Custom function, to be executed, returning the value of the last geometry intersecting with the pixel shape
CREATE OR REPLACE FUNCTION ST_FirstGeomValue4ma(pixel FLOAT, pos INTEGER[], VARIADIC args TEXT[])
    RETURNS FLOAT
    AS $$ 
    DECLARE
        pixelgeom text;
        result float4;
        query text;
    BEGIN
        -- Reconstruct the pixel shape
        pixelgeom = ST_AsText(ST_PixelAsPolygon(ST_MakeEmptyRaster(args[1]::integer, args[2]::integer, args[3]::float, args[4]::float, args[5]::float, args[6]::float, args[7]::float, args[8]::float, args[9]::integer), pos[1]::integer, pos[2]::integer));
        -- Intersects
        query = 'SELECT ' || quote_ident(args[13]) || '
                 FROM ' || quote_ident(args[10]) || '.' || quote_ident(args[11]) || ' 
                 WHERE ST_Intersects(ST_GeomFromText(' || quote_literal(pixelgeom) || ', '|| args[9] || '), ' || quote_ident(args[12]) || ') LIMIT 1';
--RAISE NOTICE 'query = %', query;
        EXECUTE query INTO result;
        RETURN result;
    END; $$
    LANGUAGE 'plpgsql' IMMUTABLE;


-- Actual creation of the raster. The year has to be provided.
CREATE TABLE snow_2010 AS
SELECT ST_MapAlgebraFct(
			rast,
			'ST_FirstGeomValue4ma(float,integer[],text[])'::regprocedure, 
			ST_Width(rast)::text, 
			ST_Height(rast)::text,
			ST_UpperLeftX(rast)::text, 
			ST_UpperLeftY(rast)::text, 
			ST_ScaleX(rast)::text, 
			ST_ScaleY(rast)::text,
			ST_SkewX(rast)::text, 
			ST_SkewY(rast)::text,
			ST_SRID(rast)::text,
			'public', 'snow', 'geom', '2010'
                       ) rast
FROM (SELECT ST_AddBand(ST_MakeEmptyRaster(89, 89, -3474442 - (26 * 198500) - (198500/2), 8833946 + (198500/2), 198500, -198500, 0, 0, 99999), '16BSI'::text, -9, -9) rast) foo;

So in your case only the last (year) parameter has to change and the query inside the custom function.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
> bounces at postgis.refractions.net] On Behalf Of JamesH
> Sent: Wednesday, April 25, 2012 1:15 PM
> To: postgis-users at postgis.refractions.net
> Subject: [postgis-users] Simple Point Density Surface
> 
> Hi all,
> 
> I am looking to take a point dataset and generate a point density raster in
> PostGIS, similar to the ArcGIS Point Density tool in Spatial Analyst.
> 
> What I believe I need to do is to define a neightbourhood around each output
> raster cell centre and calculate a density value for each pixel based on the
> number of points that fall within that neighbourhood.
> 
> Ideally to recreate the same results, I would like to define a neighbourhood
> with a radius of 8 cells - most likely a square neighbourhood.
> 
> I have been considering this for a while but I really cannot understand how
> to achieve this.
> 
> Any advise or help would be greatly appreciated.
> 
> Kind Regards,
> James Holmes
> 
> -----
> GIS Undergraduate
> --
> View this message in context: 
and create a
http://postgis.17.n6.nabble.com/Simple-Point-
> Density-Surface-tp4917343p4917343.html
> Sent from the PostGIS - User mailing list archive at Nabble.com.
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users