Calcule el ángulo de los anillos exteriores PostGIS (polígonos y multipolígonos)

enter image description here

DATA DE MUESTRA:

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

Tengo una tabla de multipolígonos y polígonos y estoy tratando de calcular los angularjs interiores de los anillos exteriores en la tabla (es decir, no hay anillos interiores …) usando ST_Azimuth . ¿Hay alguna manera de modificar la consulta adjunta para usar ST_Azimuth en sp y ep de las cadenas lineales ?

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

El uso de la function ST_Azimuth en el ejemplo aengus está dando un error porque no puede usar dos series_de_generación como parameters para una sola function. Si usa valores ep / sp en otra subconsulta, entonces no tendrá ningún problema.

Aquí está mi solución:

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

Agregué cierta complejidad ya que quería identificar cada polígono para evitar mezclar líneas para el cálculo del ángulo.

Como sqlfiddle no tiene soporte para PostGIS, he cargado este ejemplo con un código complementario a github aquí.

Puede agregar el cálculo de acimut en su subconsulta de la siguiente manera. Tenga en count que ST_Azimuth calcula el ángulo en el sentido de las agujas del reloj de abajo hacia arriba, por lo que se necesitará más trabajo para garantizar que az es realmente el ángulo interior.

  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; 

Necesitaba los angularjs internos para solo los polígonos. La solución aceptada por vamaq no me dio consistentemente los angularjs internos, pero a veces los angularjs externos, para un polígono de forma arbitraria. Esto se debe a que calcula el ángulo entre dos líneas sin considerar qué lado de esas dos líneas es interno, por lo que si el ángulo interno es obtuso, la solución de vamaq da el ángulo externo.

En mi solución, uso la regla del coseno para determinar el ángulo entre dos líneas sin reference al Azimut y luego uso ST_Contains para verificar si el ángulo del polígono es agudo u obtuso, de modo que pueda garantizar que el ángulo sea interno.

 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; 

Uso:

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