Getting multipolygon vertexes using PostGIS

Getting multipolygon vertexes using PostGIS

EN | PT

Today I needed to create a view in PostGIS that returned the vertexes of a multi-polygon layer. Besides, I needed that they were numerically ordered starting in 1, and with the respective XY coordinates.

Screenshot from 2015-11-05 23:58:19

It seemed to be a trivial task – all I would need was to use the ST_DumpPoints() function to get all vertexes – if it wasn’t for the fact that PostGIS polygons have a duplicate vertex (the last vertex must be equal to the first one) that I have no interess in showing.

After some try-and-fail, I came up with the following query:

CREATE OR REPLACE VIEW public.my_polygons_vertexes AS
WITH t AS -- Transfor polygons in sets of points
    (SELECT id_polygon,
            st_dumppoints(geom) AS dump
     FROM public.my_polygons),
f AS -- Get the geometry and the indexes from the sets of points 
    (SELECT t.id_polygon,
           (t.dump).path[1] AS part,
           (t.dump).path[3] AS vertex,
           (t.dump).geom AS geom
     FROM t)
-- Get all points filtering the last point for each geometry part
SELECT row_number() OVER () AS gid, -- Creating a unique id
       f.id_polygon,
       f.part,
       f.vertex,
       ST_X(f.geom) as x, -- Get point's X coordinate
       ST_Y(f.geom) as y, -- Get point's Y coordinate
       f.geom::geometry('POINT',4326) as geom -- make sure of the resulting geometry type
FROM f 
WHERE (f.id_polygon, f.part, f.vertex) NOT IN
      (SELECT f.id_polygon,
              f.part,
              max(f.vertex) AS max
       FROM f
       GROUP BY f.id_polygon,
                f.part);

The interesting part occurs in the WHERE clause, basically, from the list of all vertexes, only the ones not included in the list of vertexes with the maximum index by polygon part are showed, that is, the last vertex of each polygon part.

Here’s the result:

Screenshot from 2015-11-05 23:58:40

The advantage of this approach (using PostGIS) instead of using “Polygons to Lines” and “Lines to points” processing tools is that we just need to change the polygons layer, and save it, to see our vertexes get updated automatically. It’s because of this kind of stuff that I love PostGIS.

Advertisements

11 thoughts on “Getting multipolygon vertexes using PostGIS

  1. Hi!
    Thanks for valuable info. I’m not good at sql stuff – I tried doing that with my data but no luck – is it possible to know, how exactly look your public.my_polygons layer? – I mean attributes. In my sql query I tried just to change all “id_polygon” to “gid” as that’s the field in my multipolygon layer, that is unique, but that doesn’t work.
    Anyway, thanks for article.

    Like

  2. Hi!
    I managed to do it 🙂 – here are my steps for someone like me:))) :
    first I create empty multipoly layer:

    CREATE TABLE my_polygons
    (
    gid serial NOT NULL,
    geom geometry(MultiPolygon,2180),
    CONSTRAINT my_polygons_pkey PRIMARY KEY (gid)
    )
    WITH (
    OIDS=FALSE
    );
    ALTER TABLE my_polygons
    OWNER TO postgres;

    — Index: my_polygons_geom_idx

    — DROP INDEX my_polygons_geom_idx;

    CREATE INDEX my_polygons_geom_idx
    ON my_polygons
    USING gist
    (geom);

    — then I create view as you, I change just “gid” to “rn” (sth like “row number”) – you have “gid” in 13 line:
    CREATE OR REPLACE VIEW my_polygons_vertexes AS
    WITH t AS (
    SELECT my_polygons.gid,
    st_dumppoints(my_polygons.geom) AS dump
    FROM my_polygons
    ), f AS (
    SELECT t.gid,
    (t.dump).path[1] AS part,
    (t.dump).path[3] AS vertex,
    (t.dump).geom AS geom
    FROM t
    )
    SELECT row_number() OVER () AS rn,
    f.gid,
    f.part,
    f.vertex,
    st_x(f.geom) AS x,
    st_y(f.geom) AS y,
    f.geom::geometry(Point,2180) AS geom
    FROM f
    WHERE NOT ((f.gid, f.part, f.vertex) IN ( SELECT f_1.gid,
    f_1.part,
    max(f_1.vertex) AS max
    FROM f f_1
    GROUP BY f_1.gid, f_1.part));

    Like

    1. For lines it’s easier because you don’t have repeated vertex at the end. You should be able to use the same query removing the WHERE clause part, and replacing the «(t.dump).path[3] AS vertex,» line by «(t.dump).path[3] AS vertex,». Good luck!

      Like

  3. Thank you so much for this post… this is nowhere elsewhere!!!
    I am a pretty newbie so I face a problem….
    The view is being created but when I doubleclick on it appears empty. The fields are correct

    When I choose in QGIS DB manager to “Load as a new layer” , I get this message..:
    syntax error at or near “CREATE”
    LINE 1: SELECT * FROM (CREATE OR REPLACE VIEW public.my_polygons_ver…

    Can someone advise please?

    Thank you very much and a happy new year!

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s