Calculate the angle of exterior rings PostGIS (polygons & multipolygons)

Question

SAMPLE DATA:

CREATE TABLE poly_and_multipoly (
  "id" SERIAL NOT NULL PRIMARY KEY,
  "name" char(1) NOT NULL,
  "the_geom" geometry NOT NULL
);
-- add data, A is a polygon, B is a multipolygon
INSERT INTO poly_and_multipoly (name, the_geom) VALUES (
    'A', 'POLYGON((7.7 3.8,7.7 5.8,9.0 5.8,7.7 3.8))'::geometry
    ), (
    'B',
    'MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)), ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))'::geometry
);

I have a table of multipolygons and polygons and I am trying to calculate the interior angles of the exterior rings in the table (i.e. no interior rings...) using ST_Azimuth. Is there any way to modify the attached query to use ST_Azimuth on the sp and ep of the linestrings?

    SELECT id, name, ST_AsText( ST_MakeLine(sp,ep) )
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT id, name,
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
       -- extract the individual linestrings
      (SELECT id, name, (ST_Dump(ST_Boundary(the_geom))).geom
       FROM poly_and_multipoly
       ) AS linestrings
    ) AS segments;

1;"A";"LINESTRING(7.7 3.8,7.7 5.8)"
1;"A";"LINESTRING(7.7 5.8,9 5.8)"
1;"A";"LINESTRING(9 5.8,7.7 3.8)"
2;"B";"LINESTRING(0 0,4 0)"

Answer 1

The use of ST_Azimuth function at aengus example is giving an error because you can't use two generate_series as parameters for a single function. If you use ep / sp values in another subquery then you won't have any trouble.

Here is my solution:

-- 3.- Create segments from points and calculate azimuth for each line.
--     two calls of generate_series for a single function wont work (azimuth).
select id,
       name,
       polygon_num,
       point_order as vertex,
       --
       case when point_order = 1
         then last_value(ST_Astext(ST_Makeline(sp,ep))) over (partition by id, polygon_num)
         else lag(ST_Astext(ST_Makeline(sp,ep)),1) over (partition by id, polygon_num order by point_order)
       end ||' - '||ST_Astext(ST_Makeline(sp,ep)) as lines,
       --
       abs(abs(
       case when point_order = 1
         then last_value(degrees(ST_Azimuth(sp,ep))) over (partition by id, polygon_num)
         else lag(degrees(ST_Azimuth(sp,ep)),1) over (partition by id, polygon_num order by point_order)
       end - degrees(ST_Azimuth(sp,ep))) -180 ) as ang
from (-- 2.- extract the endpoints for every 2-point line segment for each linestring
      --     Group polygons from multipolygon
      select id,
             name,
             coalesce(path[1],0) as polygon_num,
             generate_series(1, ST_Npoints(geom)-1) as point_order,
             ST_Pointn(geom, generate_series(1, ST_Npoints(geom)-1)) as sp,
             ST_Pointn(geom, generate_series(2, ST_Npoints(geom)  )) as ep
      from ( -- 1.- Extract the individual linestrings and the Polygon number for later identification
             select id,
                    name,
                    (ST_Dump(ST_Boundary(the_geom))).geom as geom,
                    (ST_Dump(ST_Boundary(the_geom))).path as path -- To identify the polygon
              from poly_and_multipoly ) as pointlist ) as segments;

I've added some complexity as I wanted to identify each polygon to avoid mixing lines for angle calculation.

As sqlfiddle don't have support for PostGIS I've uploaded this example with some complementary code to github here

Answer 2

I needed the internal angles for only polygons. The accepted solution by vamaq did not consistently give me the internal angles, but sometimes external angles, for an arbitrarily shaped polygon. This is because it calculates the angle between two lines without any consideration of which side of those two lines is internal, so if the internal angle is obtuse vamaq's solution gives the external angle instead.

In my solution I use the cosine rule to determine the angle between two lines without reference to the Azimuth and then use ST_Contains to check whether the angle of the polygon is acute or obtuse so that I can ensure the angle is internal.

CREATE OR REPLACE FUNCTION is_internal(polygon geometry, p2 geometry, p3 geometry)
RETURNS boolean as
$$
    BEGIN
        return st_contains(polygon, st_makeline(p2, p3));
    END
$$ language plpgsql;


CREATE OR REPLACE FUNCTION angle(p1 geometry, p2 geometry, p3 geometry)
RETURNS float AS
$$
    DECLARE
        p12 float;
        p23 float;
        p13 float;
    BEGIN
        select st_distance(p1, p2) into p12;
        select st_distance(p1, p3) into p13;
        select st_distance(p2, p3) into p23;
        return acos((p12^2 + p13^2 - p23^2) / (2*p12*p13));
    END
$$ language plpgsql;

CREATE OR REPLACE FUNCTION internal_angle(polygon geometry, p1 geometry, p2 geometry, p3 geometry)
RETURNS float as
$$
    DECLARE
        ang float;
        is_intern boolean;
    BEGIN
        select angle(p1, p2, p3) into ang;
        select is_internal(polygon, p2, p3) into is_intern;
        IF not is_intern THEN
            select 6.28318530718 - ang into ang;
        END IF;
        return ang;
    END
$$ language plpgsql;


CREATE OR REPLACE FUNCTION corner_triplets(geom geometry)
RETURNS table(corner_number integer, p1 geometry, p2 geometry, p3 geometry) AS
$$
    DECLARE
        max_corner_number integer;
    BEGIN
        create temp table corners on commit drop as select path[2] as corner_number, t1.geom as point from (select (st_dumppoints($1)).*) as t1 where path[1] = 1;
        select max(corners.corner_number) into max_corner_number from corners;
        insert into corners (corner_number, point) select 0, point from corners where corners.corner_number = max_corner_number - 1;
        create temp table triplets on commit drop as select t1.corner_number, t1.point as p1, t2.point as p2,  t3.point as p3 from corners as t1, corners as t2, corners as t3 where t1.corner_number = t2.corner_number + 1 and t1.corner_number = t3.corner_number - 1;
        return QUERY TABLE triplets;
    END;
$$
LANGUAGE plpgsql;

CREATE OR REPLACE FUNCTION internal_angles(geom geometry)
RETURNS table(corner geometry, angle float)
AS $$
    BEGIN
        create temp table internal_angs on commit drop as select p1, internal_angle(geom, p1, p2, p3) from (select (c).* from (select corner_triplets(geom) as c) as t1) as t2;
        return QUERY TABLE internal_angs;
    END;
$$
LANGUAGE plpgsql;

Usage:

select (c).* into temp from (select internal_angles(geom) as c from my_table) as t;

Answer 3

You can add the azimuth calculation into your subquery as follows. Note that ST_Azimuth computes the angle clockwise from down to up, so more work will be required to ensure that az is actually the interior angle.

     SELECT id, name, ST_AsText( ST_MakeLine(sp,ep) ), az
    FROM
      -- extract the endpoints for every 2-point line segment for each      linestring
 (SELECT id, name,
  ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
  ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep,
    ST_Azimuth(ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)),
    ST_PointN(geom, generate_series(2, ST_NPoints(geom)  ))) as az
FROM
   -- extract the individual linestrings
  (SELECT id, name, (ST_Dump(ST_Boundary(the_geom))).geom
   FROM poly_and_multipoly
   ) AS linestrings
) AS segments;