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-boun...@postgis.refractions.net [mailto:postgis-users-
> boun...@postgis.refractions.net] On Behalf Of JamesH
> Sent: Wednesday, April 25, 2012 1:15 PM
> To: postgis-users@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@postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
