Friday, 19 October 2018

postgresql - Get the lines between all points of a polygon in PostGis : avoid nested loop?


I am trying to get all the lines between the points of a polygon, by writing PostGreSQL functions. I have written the two following functions that do the job fine :




  • the first function to get all the points from the polygon



CREATE OR REPLACE FUNCTION PointsFromPolygon(polygon geometry)
RETURNS SETOF geometry AS
$$
DECLARE
point geometry;
BEGIN

FOR point IN SELECT DISTINCT points.geom FROM ( SELECT (ST_DumpPoints(polygon)).* ) AS points LOOP
RETURN NEXT point;
END LOOP;
END;
$$
LANGUAGE plpgsql ;


  • the second function to create all segments between these points.




CREATE OR REPLACE FUNCTION AllSegmentsFromPoints(polygon geometry)
RETURNS SETOF geometry AS
$$
DECLARE
point1 geometry;
point2 geometry;
i integer;
BEGIN
i:=1;

FOR point1 IN SELECT * FROM PointsFromPolygon(polygon) LOOP
FOR point2 IN SELECT * FROM PointsFromPolygon(polygon) OFFSET i LOOP
RETURN NEXT ST_MakeLine(point1,point2);
END LOOP;
i:= i+1;
END LOOP;
END;
$$
LANGUAGE plpgsql;


Is there a better way to perform this ? Perhaps by doing a cartesian product between the points ?



Answer



I think I found my answer :


 WITH poly_geom1 AS (SELECT way FROM planet_osm_polygon WHERE id=1),
poly_geom2 AS (SELECT way FROM planet_osm_polygon WHERE id=1),
points1 AS (SELECT (ST_DumpPoints(poly_geom1.way)).* FROM poly_geom1),
points2 AS (SELECT (ST_DumpPoints(poly_geom2.way)).* FROM poly_geom2)
SELECT DISTINCT ST_MakeLine(points1.geom, points2.geom)
FROM points1,points2
WHERE points1.path <> points2.path;


Execution time for a 40 points polygon goes down from 135ms to 115ms.


Update


Thank to simplexio, the following query does exactly the same and is shorter :


WITH poly_geom AS (SELECT way FROM planet_osm_polygon WHERE id=1),
points1 AS (SELECT (ST_DumpPoints(poly_geom.way)).* FROM poly_geom),
points2 AS (SELECT (ST_DumpPoints(poly_geom.way)).* FROM poly_geom)
SELECT DISTINCT ST_MakeLine(points1.geom, points2.geom)
FROM points1,points2
WHERE points1.path <> points2.path;

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...