Re: Simple Point Density Surface

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Re: 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.
_______________________________________________
postgis-users mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/postgis-users