[postgis-users] Simple Point Density Surface

[Liglio Cavalcante]

Pierre,

I used your code to make a point density surface. But It´s not working very
well. The points used are from a table tb_lightninig and I fixed a period of
time in the function ST_PointsInPixel just for test. All the pixels in the
created table tb_densitysurface have value equal to zero, and for sure there
are points in this surface. For test used the RAISE NOTICE, and I also
received 1.001 notices during the execution (My raster has 34.225 pixels,
why only 1.001 querys ? It is not one query per pixel ?). example:

NOTICE:  query = SELECT Count(*) 
        FROM public.tb_lightning 
        WHERE ST_Intersects(ST_GeomFromText('POLYGON((-32.5499940183796
-36.1707879633483,-32.3247687931544 -36.1707879633483,-32.3247687931544
-36.3960131885735,-32.5499940183796 -36.3960131885735,-32.5499940183796
-36.1707879633483))', 4326), point_lightning) 
        AND din_lightning BETWEEN '2016-09-01 00:00:00' AND '2016-11-01
00:00:00' 
CONTEXT:  SQL function "st_mapalgebrafct" during inlining

This query returns value 1, so I don't understand why my raster has only
values equal to zero.

My Code:

CREATE OR REPLACE FUNCTION ST_PointsInPixel(pixel FLOAT, pos INTEGER[],
VARIADIC args TEXT[]) RETURNS INTEGER AS $$ 
DECLARE 
pixelgeom text; 
result int4; 
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 ' || 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]) || ') ' || '
        AND din_lightning BETWEEN ''2016-09-01 00:00:00'' AND ''2016-11-01
00:00:00'' ';

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 tb_densitysurface AS 
SELECT ST_MapAlgebraFct( 
                        rast, 
                       
'ST_PointsInPixel(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', 'tb_lightning', 'point_lightning',
'Count(*)' 
                        ) rast 
FROM (SELECT ST_AddBand(ST_MakeEmptyRaster(185, 185, -73.991435459821,
5.27065347809308, 0.2252252252252252, -0.2252252252252252, 0, 0, 4326),
'32BSI'::text, 0, -1) rast) foo;




--
View this message in context: http://postgis.17.x6.nabble.com/Simple-Point-Density-Surface-tp4917343p5010775.html
Sent from the PostGIS - User mailing list archive at Nabble.com.