计算PostGIS（多边形和多重多边形）外部环的角度。

Question

样本数据：

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
);

我有一个多边形和单独的多边形表格，想要计算表格中外环的内角（即没有内环...），使用ST_Azimuth。是否有任何方法修改附加查询，以在线串的sp和ep上使用ST_Azimuth？

    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)"

- user14696

Answer 1

我只需要多边形的内角。vamaq 提供的解决方案并没有始终给出任意形状多边形的内角, 有时会给出外角。这是因为它在计算两条线之间的夹角时没有考虑这两条线的哪一侧是内部的, 因此如果内角是钝角，vamaq 的解决方案会给出外角。

我的解决方案使用余弦定理来确定两条线之间的角度而不涉及方位角，然后使用 ST_Contains 检查多边形的角度是锐角还是钝角，以确保角度是内角。

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;

使用方法：

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

- public static void

Answer 2

在aengus示例中使用ST_Azimuth函数会出错，因为你不能将两个generate_series作为单个函数的参数使用。如果在另一个子查询中使用ep / sp值，则不会有任何问题。

以下是我的解决方案：

-- 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;

我增加了一些复杂性，因为我想识别每个多边形以避免在计算角度时混淆线条。

由于sqlfiddle不支持PostGIS，因此我已经上传了这个示例，并补充了一些代码到github这里。

- vamaq

---

哇！这非常有帮助（而且代码逻辑非常易懂）。谢谢。 - user14696

很好的答案。由于某种原因，第一个角度会随机失败（如果我只为该行运行它，则可以工作，但如果我为整个表运行它，则会失败）。我通过将选择从“from poly_and_multipoly”移动到第一级并使用左连接侧来修复它。不确定为什么这会有所帮助。 - Vesanto

Answer 3

您可以将方位角计算添加到子查询中，如下所示。请注意，ST_Azimuth从下到上顺时针计算角度，因此需要更多的工作来确保az实际上是内角。

     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;

- aengus

---

谢谢。嗯，我用那个出错了...看来我得调试一下...“错误：函数和运算符最多只能接受一个集合参数”。 - user14696