SpatialDB Advisor
I find having a vectorising function in my database kit-bag almost indispensable in my work. I have written and revised my Oracle PL/SQL functions that do this many times over the years and so, when I first started using PostGIS, I decided to learn how to program in PL/pgSQL by implementing something familiar: a GetVector() function.
Now, my PostgreSQL/PostGIS is not bad, but I am always willing to learn from others who are more expert than me: so don’t take the following function (which , for example, doesn’t handle Curves as they are still in their infancy in PostGIS) as being perfect or correct in every manner.
Firstly, we need to create two data types that we can use.
A Coordinate Type
CREATE TYPE coordtype AS (
x double precision,
y double precision,
z double precision,
m double precision
);
Query returned successfully with no result in 78 ms.
A Vector Type based on the Coordinate Type
CREATE TYPE vectortype AS (
startcoord coordtype,
endcoord coordtype
);
Query returned successfully with no result in 16 ms.
What I tried to do in the following GetVector() function is handle all geometry types in the one function. Initially I used one cursor per geometry type (and a single refCursor for processing) but in this implementation I have chosen to use a single cursor with UNION ALLs between the SQL that processes the vertices of each geometry type. I am hoping that the PostgreSQL query optimiser does its equivalent of Oracle’s partition pruning. Perhaps someone may test this aspect and report if it is not performing as fast as one would hope.
CREATE OR REPLACE FUNCTION ST_GetVector(p_geometry geometry)
RETURNS SETOF VectorType IMMUTABLE
AS $$
DECLARE
v_GeometryType varchar(1000);
v_rec RECORD;
v_vector VectorType;
v_start CoordType;
v_end CoordType;
c_points CURSOR ( p_geom geometry,
p_Geometrytype text )
IS
SELECT ST_X(sp) as sx,ST_Y(sp) as sy,ST_Z(sp) as sz,ST_M(sp) as sm,
ST_X(ep) as ex,ST_Y(ep) as ey,ST_Z(ep) as ez,ST_M(ep) as em
FROM ( SELECT pointn(p_geom, generate_series(1, npoints(p_geom)-1)) as sp,
pointn(p_geom, generate_series(2, npoints(p_geom) )) as ep
WHERE ( p_GeometryType = 'ST_LineString' )
) AS linestring
UNION ALL
SELECT ST_X(sp),ST_Y(sp),ST_Z(sp),ST_M(sp),ST_X(ep),ST_Y(ep),ST_Z(ep),ST_M(ep)
FROM ( SELECT ST_PointN(b.geom, generate_series(1, ST_NPoints(b.geom)-1)) as sp,
ST_PointN(b.geom, generate_series(2, ST_NPoints(b.geom) )) as ep
FROM (SELECT ST_GeometryN(p_geom,generate_series(1,ST_NumGeometries(p_geom))) as geom) as b
WHERE ( p_GeometryType = 'ST_MultiLineString' )
) AS multiline
UNION ALL
SELECT ST_X(sp),ST_Y(sp),ST_Z(sp),ST_M(sp),ST_X(ep),ST_Y(ep),ST_Z(ep),ST_M(ep)
FROM ( SELECT ST_PointN(a.geom, generate_series(1, ST_NPoints(a.geom)-1)) as sp,
ST_PointN(a.geom, generate_series(2, ST_NPoints(a.geom) )) as ep
FROM ( SELECT ST_ExteriorRing(p_geom) as geom
WHERE ( p_GeometryType = 'ST_Polygon' )
UNION ALL
SELECT ST_InteriorRingN(p_geom,generate_series(1,ST_NumInteriorRings(p_geom))) as geom
WHERE ( p_GeometryType = 'ST_Polygon' )
) a
) as polygon
UNION ALL
SELECT ST_X(sp),ST_Y(sp),ST_Z(sp),ST_M(sp),ST_X(ep),ST_Y(ep),ST_Z(ep),ST_M(ep)
FROM ( SELECT ST_pointn(a.geom, generate_series(1, ST_npoints(a.geom)-1)) as sp,
ST_pointn(a.geom, generate_series(2, ST_npoints(a.geom) )) as ep
FROM ( SELECT ST_exteriorring(b.geom) as geom
FROM (select ST_geometryn(p_geom,generate_series(1,ST_numgeometries(p_geom))) as geom ) b
WHERE ( p_GeometryType = 'ST_MultiPolygon' )
UNION ALL
SELECT ST_InteriorRingn(c.geom,generate_series(1,ST_NumInteriorRings(c.geom))) as geom
FROM (select ST_GeometryN(p_geom,generate_series(1,ST_NumGeometries(p_geom))) as geom ) c
WHERE ( p_GeometryType = 'ST_MultiPolygon' )
) a
) as multipoly;
Begin
If ( p_geometry is NULL ) Then
return;
End If;
v_GeometryType := ST_GeometryType(p_geometry);
If ( v_GeometryType in ('ST_Point','ST_MultiPoint') ) Then
return;
End If;
IF ( v_GeometryType = 'ST_Geometry' ) THEN
-- Could be anything... use native PostGIS function
v_GeometryType = GeometryType(p_geometry);
IF ( v_geometryType = 'GEOMETRYCOLLECTION' ) THEN
FOR v_geom IN 1..ST_NumGeometries(p_geometry) LOOP
FOR v_rec IN SELECT * FROM ST_GetVector(ST_GeometryN(p_geometry,v_geom)) LOOP
RETURN NEXT v_rec;
END LOOP;
END LOOP;
ELSE
-- Probably CURVED something...
RETURN;
END IF;
END IF;
IF ( v_geometryType NOT IN ('ST_Geometry','GEOMETRYCOLLECTION') ) THEN
OPEN c_points(p_geometry,v_Geometrytype);
LOOP
FETCH c_points INTO
v_start.x, v_start.y, v_start.z, v_start.m,
v_end.x, v_end.y, v_end.z, v_end.m;
v_vector.startcoord := v_start;
v_vector.endcoord := v_end;
EXIT WHEN NOT FOUND;
RETURN NEXT v_vector;
END LOOP;
CLOSE c_points;
END IF;
end;
$$ LANGUAGE 'plpgsql';
Query returned successfully with no result in 156 ms.
Now, so that we can use this function in nested table Select statements we construct a “wrapper”.
CREATE OR REPLACE FUNCTION ST_GetVectorSQL(p_geometry geometry)
RETURNS SETOF VectorType IMMUTABLE
AS
$$
SELECT * FROM ST_GetVector($1);
$$
LANGUAGE 'sql';
Query returned successfully with no result in 31 ms.
Now let’s conduct some simple tests.
select * from ST_GetVector('GEOMETRYCOLLECTION(POINT(2 3 4),LINESTRING(2 3 4,3 4 5))'::geometry);
| startcoord coordtype |
endcoord coordtype |
|---|---|
| (2,3,4,) | (3,4,5,) |
select * from ST_GetVector('LINESTRING(0 0, 1 1, 2 2, 3 3)'::geometry) as geom;
| startcoord coordtype |
endcoord coordtype |
|---|---|
| (0,0,,) | (1,1,,) |
| (1,1,,) | (2,2,,) |
| (2,2,,) | (3,3,,) |
select * from ST_GetVector('MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))'::geometry) As GV;
| startcoord coordtype |
endcoord coordtype |
|---|---|
| (0,0,,) | (1,1,,) |
| (1,1,,) | (1,2,,) |
| (2,3,,) | (3,2,,) |
| (3,2,,) | (5,4,,) |
select * from ST_GetVector('POLYGON((326454.7 5455793.7,326621.3 5455813.7,326455.4 5455796.6,326454.7 5455793.7))'::geometry);
| startcoord coordtype |
endcoord coordtype |
|---|---|
| (326454.7,5455793.7,,) | (326621.3,5455813.7,,) |
| (326621.3,5455813.7,,) | (326455.4,5455796.6,,) |
| (326455.4,5455796.6,,) | (326454.7,5455793.7,,) |
select * from ST_GetVector('MULTIPOLYGON(((326454.7 5455793.7,326621.3 5455813.7,326455.4 5455796.6,326454.7 5455793.7)),((326771.6 5455831.6,326924.1 5455849.9,326901.9 5455874.2,326900.7 5455875.8,326888.9 5455867.3,326866 5455853.1,326862 5455851.2,326847.4 5455845.8,326827.7 5455841.2,326771.6 5455831.6)))'::geometry);
| startcoord coordtype |
endcoord coordtype |
|---|---|
| (326454.7,5455793.7,,) | (326621.3,5455813.7,,) |
| (326621.3,5455813.7,,) | (326455.4,5455796.6,,) |
| (326455.4,5455796.6,,) | (326454.7,5455793.7,,) |
| (326771.6,5455831.6,,) | (326924.1,5455849.9,,) |
| (326924.1,5455849.9,,) | (326901.9,5455874.2,,) |
| (326901.9,5455874.2,,) | (326900.7,5455875.8,,) |
| (326900.7,5455875.8,,) | (326888.9,5455867.3,,) |
| (326888.9,5455867.3,,) | (326866,5455853.1,,) |
| (326866,5455853.1,,) | (326862,5455851.2,,) |
| (326862,5455851.2,,) | (326847.4,5455845.8,,) |
| (326847.4,5455845.8,,) | (326827.7,5455841.2,,) |
| (326827.7,5455841.2,,) | (326771.6,5455831.6,,) |
Table Tests
Now we can test our function against some real tables and SQL.
DROP TABLE v_polygons;
| addgeometrycolumn text |
|---|
| public.v_polygons.geom SRID:28355 TYPE:MULTIPOLYGON DIMS:2 |
ALTER TABLE v_polygons DROP CONSTRAINT enforce_geotype_geom;
| oid integer |
startcoord coordtype |
endcoord coordtype |
|---|---|---|
| 1 | (326454.7,5455793.7,,) | (326621.3,5455813.7,,) |
| 1 | (326621.3,5455813.7,,) | (326455.4,5455796.6,,) |
| 1 | (326455.4,5455796.6,,) | (326454.7,5455793.7,,) |
| 1 | (326771.6,5455831.6,,) | (326924.1,5455849.9,,) |
| 1 | (326924.1,5455849.9,,) | (326901.9,5455874.2,,) |
| 1 | (326901.9,5455874.2,,) | (326900.7,5455875.8,,) |
| 1 | (326900.7,5455875.8,,) | (326888.9,5455867.3,,) |
| 1 | (326888.9,5455867.3,,) | (326866,5455853.1,,) |
| 1 | (326866,5455853.1,,) | (326862,5455851.2,,) |
| 1 | (326862,5455851.2,,) | (326847.4,5455845.8,,) |
| 1 | (326847.4,5455845.8,,) | (326827.7,5455841.2,,) |
| 1 | (326827.7,5455841.2,,) | (326771.6,5455831.6,,) |
| 2 | (326667.8,5455760.3,,) | (326669.7,5455755.6,,) |
| 2 | (326669.7,5455755.6,,) | (326688.1,5455766,,) |
| 2 | (326688.1,5455766,,) | (326685.3,5455770.5,,) |
| 2 | (326685.3,5455770.5,,) | (326667.8,5455760.3,,) |
Let’s visualise this by moving the resultant vectors parallel 10 meters from the original polygon boundaries.
We will use Regina Obe’s upgis_lineshift function to move individual polygon vectors to the right by 10m.
CREATE OR REPLACE FUNCTION upgis_lineshift(centerline geometry, dist double precision)
RETURNS geometry AS
$$
DECLARE
delx float;
dely float;
x0 float;
y0 float;
x1 float;
y1 float;
az float;
dir integer;
line geometry;
BEGIN
az := ST_Azimuth (ST_StartPoint(centerline), ST_EndPoint(centerline));
dir := CASE WHEN az < pi() THEN -1 ELSE 1 END;
delx := ABS(COS(az)) * dist * dir;
dely := ABS(SIN(az)) * dist * dir;
Then we need a table to hold the vectors we create and shift.
DROP TABLE polygon_vectors;
| addgeometrycolumn text |
|---|
| public.polygon_vectors.geom SRID:28355 TYPE:LINESTRING DIMS:2 |
Finally, let’s create, shift and write the vectors.
INSERT INTO polygon_vectors (geom)
SELECT upgis_lineshift(ST_SetSRID(ST_MakeLine(ST_MakePoint((startcoord).x,(startcoord).y)::geometry,
ST_MakePoint((endcoord).x, (endcoord).y)::geometry),
28355),
10)
FROM (SELECT (ST_GetVector(v_polygons.geom)).*
FROM v_polygons ) as v;
| count bigint |
|---|
| 16 |
Finally, let’s visualise the polygons and the shifted vectors.

I hope this function, and article, is of use to someone.


















<<How to apply spatial constraints to PostGIS tables >>FOSSLIC - PostGIS presentations