Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Vectorization: Exploding a linestring or polygon into individual vectors in PostGIS

Thursday May 14 2009 at 08:03

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 ST_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_Vectorize(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 ) 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 ST_PointN(p_geom, generate_series(1, ST_NPoints(p_geom)-1)) as sp, ST_PointN(p_geom, generate_series(2, ST_NPoints(p_geom) )) as ep WHERE ST_GeometryType(p_geom) = ‘ST_LineString’ UNION ALL 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 WHERE ST_GeometryType(p_geom) = ‘ST_MultiLineString’ ) as b UNION ALL 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 UNION ALL SELECT ST_InteriorRingN(p_geom,generate_series(1,ST_NumInteriorRings(p_geom))) as geom ) a WHERE ST_GeometryType(p_geom) = ‘ST_Polygon’ UNION ALL 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) as b 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) as c ) a WHERE ST_GeometryType(p_geom) = ‘ST_MultiPolygon’ ) as f;
Begin If ( p_geometry is NULL ) Then return; End If; v_GeometryType := ST_GeometryType(p_geometry); — RAISE NOTICE ‘ST_GeometryType %’, v_GeometryType; IF ( v_GeometryType in (‘ST_Point’,‘ST_MultiPoint’,‘ST_Geometry’) ) Then return; END IF; IF ( v_GeometryType IN (‘ST_GeometryCollection’) ) THEN FOR v_geom IN 1..ST_NumGeometries(p_geometry) LOOP FOR v_rec IN SELECT * FROM ST_Vectorize(ST_GeometryN(p_geometry,v_geom)) LOOP RETURN NEXT v_rec; END LOOP; END LOOP; ELSE OPEN c_points(p_geometry); 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’;

CREATE OR REPLACE FUNCTION ST_VectorizeSQL(p_geometry geometry) RETURNS SETOF VectorType IMMUTABLE
AS
$$ SELECT * FROM ST_Vectorize($1);
$$
LANGUAGE ‘sql’;
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;

Query returned successfully with no result in 16 ms.

CREATE TABLE v_polygons
(
  gid serial NOT NULL primary key
)
WITH (
  OIDS=FALSE
);

NOTICE:  CREATE TABLE will create implicit sequence "v_polygons_gid_seq" for serial column "v_polygons.gid"
NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "v_polygons_pkey" for table "v_polygons"
Query returned successfully with no result in 125 ms.

ALTER TABLE v_polygons OWNER TO postgres;

Query returned successfully with no result in 16 ms.

SELECT addGeometryColumn('public','v_polygons','geom','28355','MULTIPOLYGON','2'); 
Total query runtime: 78 ms.
1 rows retrieved.

addgeometrycolumn
text
public.v_polygons.geom SRID:28355 TYPE:MULTIPOLYGON DIMS:2

ALTER TABLE v_polygons DROP CONSTRAINT enforce_geotype_geom;

Query returned successfully with no result in 16 ms.

ALTER TABLE v_polygons ADD  CONSTRAINT enforce_geotype_geom CHECK ((geometrytype(geom) = ANY (ARRAY['MULTIPOLYGON'::text, 'POLYGON'::text])) OR geom IS NULL);

INSERT INTO v_polygons (geom) 
VALUES 
(ST_GeomFromText('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)))',28355)),
(ST_PolygonFromText('POLYGON((326667.8 5455760.3, 326669.7 5455755.6, 326688.1 5455766.0, 326685.3 5455770.5, 326667.8 5455760.3))',28355));

SELECT v_polygons.gid, 
       (ST_GetVectorSQL(v_polygons.geom)).*
  FROM v_polygons;

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;

        IF az > pi()/2 AND az < pi() OR az > 3 * pi()/2 THEN
                line := ST_Translate(centerline, delx, dely) ;
        ELSE
                line := ST_Translate(centerline, -delx, dely);
        END IF;

        RETURN line;
  END;
  $$
  LANGUAGE 'plpgsql' IMMUTABLE;
  COMMENT ON FUNCTION upgis_lineshift(geometry, double precision) IS 'Takes a 2D line string and shifts it dist units along  the perpendicular defined by the straight line between the start and end point.
Convention: (right is positive and left is negative. Right being defined as to right of observer standing at start point and looking down the end point)';

Query returned successfully with no result in 78 ms.

Then we need a table to hold the vectors we create and shift.

DROP   TABLE polygon_vectors;

Query returned successfully with no result in 31 ms.

CREATE TABLE polygon_vectors
(
  gid serial NOT NULL primary key
)
WITH (
  OIDS=FALSE
);

NOTICE:  CREATE TABLE will create implicit sequence "polygon_vectors_gid_seq" for serial column "polygon_vectors.gid"
NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "polygon_vectors_pkey" for table "polygon_vectors"
Query returned successfully with no result in 125 ms.

ALTER TABLE v_polygons OWNER TO postgres;

Query returned successfully with no result in 16 ms.

SELECT addGeometryColumn('public','polygon_vectors','geom','28355','LINESTRING','2'); 

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;

Query returned successfully: 16 rows affected, 15 ms execution time.

SELECT COUNT(*)
  FROM polygon_vectors;

Total query runtime: 16 ms.
1 rows retrieved.

count
bigint
16


Finally, let’s visualise the polygons and the shifted vectors.

Polygons and Vectors

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

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment [9]

Hello Simon Greener!

I tried to include your function to explode a multilinestring table in Postgis, but I’ve got the message ‘ERROR: type “vectortype” do not exists’. Do you know what is this error, or wath I’m forgetting to do? Thank you!

— Felipe · 28 June 2011, 10:20 · #

Felipe,

At the beginning of the article two types have to be created:

  • coordType
  • vectorType

Create them and then compile again.

regards
Simon

Simon Greener · 8 July 2011, 06:55 · #

thank you Simon for a great article. May I ask you how you did the visualization (at the bottom of the blog)?
Many thanks
George

— george wash · 19 January 2012, 21:33 · #

George,

That’s a very good question! ‘Twas a while ago. I suspect that what I did was link Manifold GIS to the created tables (which is why I created the tables). Equally you could use QGIS without any difficulty. Finally, the new PostGISViewer plugin for pgAdmin III would IMHO be best though I have not been able to get it to work on Windows 7 64 bit.

See: http://www.postgis.us/page_desktop_tools

regards
Simon

Simon Greener · 20 January 2012, 00:22 · #

Thank you Simon for the quick reply. i installed the postgisviewer on my Windows 7/64 and I get a flash on the screen and no more. But what I really wanted to ask you,as a novice with Postgis, how to call your function with a table name eg:

SELECT * from ST_GetVector(select the_geom from rivers where name=‘wrl24’);

But I only get an error:

ERROR: syntax error at or near “select”
LINE 1: select * from ST_GetVector(select the_geom from rivers where… ^

  • Error **********

ERROR: syntax error at or near “select”
SQL state: 42601
Character: 28

— george wash · 20 January 2012, 06:36 · #

George,

To execute the function against a whole table you would do the following:

SELECT (ST_GetVectorSQL(rivers.the_geom)).* FROM rivers WHERE name=‘wrl24’;

regards
Simon

Simon Greener · 20 January 2012, 06:54 · #

Hi thanks for this very useful.

In the function:
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

I believe pointn( should be ST_PointN, and npoints( should be ST_Npoints.
Running it against a line dataset I also removed the other three types, as I got an error about ST_ExteriorRing not existing in a line geometry.

Thanks.

— Heikki · 21 September 2015, 02:18 · #

Heikki,

The PostGIS API must have changed since I wrote this function. I have updated the function and tested it against:

POSTGIS=“2.1.7 r13414” GEOS=“3.4.2-CAPI-1.8.2 r3924” PROJ=“Rel. 4.8.0, 6 March 2012” GDAL=“GDAL 1.11.1, released 2014/09/24” LIBXML=“2.7.8” LIBJSON=“UNKNOWNTOPOLOGY RASTER

And it now performs correctly.

Thanks.

Regards
Simon

— Simon Greener · 21 September 2015, 21:06 · #

Hi Simon,
I used a modified version of your function in this article: http://community.jaspersoft.com/wiki/how-generate-map-data-postgis-or-shapefiles-use-reports
I hope that’s OK with you?

Thanks!

— Kamal · 20 November 2015, 02:25 · #