Crear un polígonos a partir de puntos

Hace un par de días que empecé con esta tarea, que al principio parecía muy fácil. Es más seguro que ustedes piensan en este momento que es fácil. Pero es realmente algo bastante difícil de lograr, de hecho no encontré nada que sirva más que un par de select anidados con un makepolygon al final. Creo que este es el triunfo del mes, o al menos hasta ahora. Por tal motivo, y porque no quiero volver a ver esto por un tiempo, por favor no me pregunten sobre esto, si lo entienden me alegro mucho y espero que les ayude, caso contrario lean de nuevo.

El escenario es el siguiente, imaginense un par de puntos desornedados y potencialmente repetidos con los que nosotros tenemos que “adivinar” el ordén, y crear un polígono. Para explicarlo mejor:

Puntos desordenados

Puntos desordenados

Ahora bien, el problema viene porque, si bien para un humano es fácil darse cuenta que el polígono que quiero formar es el formado por los segmentos A-D-E-F-G-H-C-B-A para postgis no es tan sencillo darse cuenta, ya que al no estar ordenados los polígonos podría ser cualquier combinación de ellos. Hasta incluso armar varios polígonos como por ejemplo A-D-B por un lado, y C-H-G-F-E-C por el otro.

Si siguen pensando que esto es sencillo de resolver, piénsenlo de nuevo.

El escenario es el siguiente:

tengo una tabla zonas_usuales de la siguiente forma:

CREATE TABLE zonas_usuales
(
  nombre character varying(40),
  x double precision,
  y double precision
)
WITH (
  OIDS=TRUE
);

Bien el principal problema es que los puntos no están ordenados, o sea, nada impide usar el A-D-E-F-G-c-h-B-A (cambié dos letras) que no da como resultado ningún polígono.

Entonces, a crear funciones:

CREATE OR REPLACE FUNCTION fneiff(geometry,geometry,geometry) RETURNS float as $$
DECLARE
angle decimal;
BEGIN
select abs((ST_Azimuth($1,$2)) - (ST_Azimuth($1,$3))) into angle;
if (angle > pi()) then
select (2*pi())-angle into angle;
end if;
RETURN (select angle / ((distance($1,$2) * distance($1,$3))));
END;
$$ language 'plpgsql';

Esta primera función es la que hace la mágina de seleccionar los dos puntos, que están más cerca y en un ángulo correcto a un punto dado.

CREATE OR REPLACE FUNCTION esunico(geometry, text) RETURNS boolean as $$
BEGIN
RETURN exists(select line from
	(select distinct MakeLine(point1,point2) as line from
	(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point1 from zonas_usuales where nombre = $2) as p1,
	(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point2 from zonas_usuales where nombre = $2) as p2
	where point1 <> point2 and MakeLine(point1,point2) = $1 limit 1) as tabla
	where astext(line) = astext($1));
END;
$$ language 'plpgsql';

Esta segunda función me dirá si un segmento es único o no, porque ambos segmentos A-B y B-A son segmentos válidos, por lo cual uso esta función para quitar uno de ellos.

CREATE OR REPLACE FUNCTION nocruza(geometry, geometry, text) RETURNS boolean as $$
BEGIN
return exists(select point1 from (
select distinct fneiff ($1,t2.the_geom, t3.the_geom) as func, t2.the_geom as point1, t3.the_geom as point2 from
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as the_geom from zonas_usuales where nombre = $3) as t2,
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as the_geom from zonas_usuales where nombre = $3) as t3
where fneiff($1,t2.the_geom, t3.the_geom) is not null and
(essolucion(t2.the_geom,$1,$3) or essolucion(t3.the_geom,$1,$3)) and astext(t3.the_geom) <> astext(t2.the_geom) order by func desc limit 2) as asdf
where astext(point1) = astext($2));
END;
$$ language 'plpgsql';
CREATE OR REPLACE FUNCTION essolucion(geometry,geometry, text) RETURNS boolean as $$
BEGIN
return exists(select point1 from (
select distinct fneiff ($1,t2.the_geom, t3.the_geom) as func, t2.the_geom as point1, t3.the_geom as point2 from
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as the_geom from zonas_usuales where nombre = $3) as t2,
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as the_geom from zonas_usuales where nombre = $3) as t3
where fneiff($1,t2.the_geom, t3.the_geom) is not null and astext(t3.the_geom) <> astext(t2.the_geom) order by func desc limit 2) as asdf
where astext(point1) = astext($2));
END;
$$ language 'plpgsql';

no cruza, llama a mi superfunción fneiff y se asegura de tomar puntos que no sean bucles por ejemplo.

CREATE OR REPLACE FUNCTION dossegmentosmascortos(geometry,geometry, text) RETURNS boolean as $$
BEGIN
RETURN exists(select line from
(select distinct MakeLine(point1,point2) as line,length(MakeLine(point1,point2)) from
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point1 from zonas_usuales where nombre = $3) as p1,
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point2 from zonas_usuales where nombre = $3) as p2
where point1 <> point2 and point1 = $1 and esunico(MakeLine(point1,point2),$3) and nocruza(point1,point2,$3) order by length(MakeLine(point1,point2)) limit 2) as tabla where astext(line) = astext(MakeLine($1,$2)));
END;
$$ language 'plpgsql';

Esta función me busca:
+ Los puntos que son únicos
+ Solo dos puntos
+ y que no se crucen entre sí (otra vez con la función mágica)

CREATE OR REPLACE FUNCTION nombre2polygon(text) RETURNS geometry as $$
BEGIN
RETURN (select MakePolygon(ST_LineMerge(ST_Collect(line))) as arrlin from
(select distinct MakeLine(point1,point2) as line from
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point1 from zonas_usuales where nombre = $1) as p1,
(select distinct transform(geomfromtext('POINT ('|| x || ' ' || y || ')',22175),4326) as point2 from zonas_usuales where nombre = $1) as p2
where point1 <> point2 and dossegmentosmascortos(point1,point2,$1)) as tabla
);
END;
$$ language 'plpgsql';

Solo queda el broche de oro, la función que une todos los segmentos en uno solo y llama a makepolygon, ahora que estoy terminando el post no puedo creer que me llevó 15m escribirlo y unos cuatro días codearlo. Espero les ahorre un par de horas de sueño. Siempre tengan en cuenta que en mi empresa ofrecemos el servicio de hosting gis, en esta comparativa pueden ver una lista completa de precios y servicios.

Tags: , , , , , , , , ,

4 Responses to “Crear un polígonos a partir de puntos”

  1. Hector Cantor says:

    Que tal, estoy tratando de aprender a usar Sistemas Geograficos, estoy tratando de aprender a programar y se me hace interesante lo que presentas en tu blog, para crear los poligonos. No entiendo mucho lo que haces en el postgresql, pero mencionas que te falta una funcion que se llama makepolygon. Seria posible que me explicaras como integrar el procedimiento para crear poligonos??

    Gracias.

    • Matias Neiff says:

      Hola Hector, perdona si no se entendió, está todo ahí, por el broche de oro me refería a la última función nombre2polygon.

      Un abrazo.
      MN

  2. JORGE says:

    Hola
    He cambiado la funcion a lo siguiente y no funciona:
    cual seria el procedimiento que esta mal
    te agradezco la atencion
    –tabla de ejemplo
    /*
    CREATE TABLE zonas_usuales
    (
    nombre character varying(40),
    x double precision,
    y double precision
    )
    WITH (
    OIDS=TRUE
    );
    */

    –Problema :
    ————————————————————————–
    –Bien el principal problema es que los puntos no están ordenados, o sea,
    –nada impide usar el A-D-E-F-G-c-h-B-A (cambié dos letras) que no da como
    –resultado ningún polígono.
    —————————————————————————-

    ————————————————————–
    –Esta primera función selecciona los dos puntos,
    –que están más cerca y en un ángulo correcto a un punto dado.
    —————————————————————-
    CREATE OR REPLACE FUNCTION fneiff(geometry,geometry,geometry) RETURNS float as $$
    DECLARE
    angle decimal;
    BEGIN
    select abs((ST_Azimuth($1,$2)) – (ST_Azimuth($1,$3))) into angle;
    if (angle > pi()) then
    select (2*pi())-angle into angle;
    end if;
    RETURN (select angle / ((distance($1,$2) * distance($1,$3))));
    END;
    $$ language ‘plpgsql’;

    ——————————————————————————
    –Esta segunda función verifica si un segmento es único o no, porque ambos
    –segmentos A-B y B-A son segmentos válidos,
    ——————————————————————————–
    CREATE OR REPLACE FUNCTION esunico(geometry, text) RETURNS boolean as $$
    BEGIN
    RETURN exists(select line from
    (select distinct ST_ST_MakeLine(point1,point2) as line from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point1 from zonas_usuales where nombre = $2) as p1,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point2 from zonas_usuales where nombre = $2) as p2
    where ST_Equals(point1,point2) ‘t’ and ST_MakeLine(point1,point2) = $1 limit 1) as tabla
    where ST_AsText(line) = ST_AsText($1));
    END;
    $$ language ‘plpgsql’;

    ————————————————————
    –no cruza, llama a mi superfunción fneiff y se asegura
    –de tomar puntos que no sean bucles por ejemplo.
    ————————————————————-

    CREATE OR REPLACE FUNCTION nocruza(geometry, geometry, text) RETURNS boolean as $$
    BEGIN
    return exists(select point1 from (
    select distinct fneiff ($1,t2.the_geom, t3.the_geom) as func, t2.the_geom as point1, t3.the_geom as point2 from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as the_geom from zonas_usuales where nombre = $3) as t2,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as the_geom from zonas_usuales where nombre = $3) as t3
    where fneiff($1,t2.the_geom, t3.the_geom) is not null and
    (essolucion(t2.the_geom,$1,$3) or essolucion(t3.the_geom,$1,$3)) and ST_AsText(t3.the_geom) ST_AsText(t2.the_geom) order by func desc limit 2) as asdf
    where ST_AsText(point1) = ST_AsText($2));
    END;
    $$ language ‘plpgsql’;
    CREATE OR REPLACE FUNCTION essolucion(geometry,geometry, text) RETURNS boolean as $$
    BEGIN
    return exists(select point1 from (
    select distinct fneiff ($1,t2.the_geom, t3.the_geom) as func, t2.the_geom as point1, t3.the_geom as point2 from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as the_geom from zonas_usuales where nombre = $3) as t2,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as the_geom from zonas_usuales where nombre = $3) as t3
    where fneiff($1,t2.the_geom, t3.the_geom) is not null and ST_AsText(t3.the_geom) ST_AsText(t2.the_geom) order by func desc limit 2) as asdf
    where ST_AsText(point1) = ST_AsText($2));
    END;
    $$ language ‘plpgsql’;

    —————————————————————————-
    –Esta función busca:
    –+ Los puntos que son únicos
    –+ Solo dos puntos
    –+ y que no se crucen entre sí
    —————————————————————————-

    CREATE OR REPLACE FUNCTION dossegmentosmascortos(geometry,geometry, text) RETURNS boolean as $$
    BEGIN
    RETURN exists(select line from
    (select distinct ST_MakeLine(point1,point2) as line,ST_Length(ST_MakeLine(point1,point2)) from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point1 from zonas_usuales where nombre = $3) as p1,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point2 from zonas_usuales where nombre = $3) as p2
    where ST_Equals(point1,point2) ‘t’ and point1 = $1 and esunico(ST_MakeLine(point1,point2),$3) and nocruza(point1,point2,$3) order by ST_Length(ST_MakeLine(point1,point2)) limit 2) as tabla where ST_AsText(line) = ST_AsText(ST_MakeLine($1,$2)));
    END;
    $$ language ‘plpgsql’;

    ————————————————————————
    –función que une todos los segmentos
    –en uno solo y llama a makepolygon
    ————————————————————————-
    /*
    CREATE OR REPLACE FUNCTION nombre2polygon(text) RETURNS geometry as $$
    BEGIN
    RETURN (select MakePolygon(ST_LineMerge(ST_Collect(line))) as arrlin from
    (select distinct ST_MakeLine(point1,point2) as line from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point1 from zonas_usuales where nombre = $1) as p1,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point2 from zonas_usuales where nombre = $1) as p2
    where point1 point2 and dossegmentosmascortos(point1,point2,$1)) as tabla
    );
    END;
    $$ language ‘plpgsql’;
    */

    CREATE OR REPLACE FUNCTION nombre2polygon(text) RETURNS geometry as $$
    BEGIN
    RETURN (select ST_MakePolygon(ST_LineMerge(ST_Collect(line))) as arrlin from
    (select distinct ST_MakeLine(point1,point2) as line from
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point1 from zonas_usuales where nombre = $1) as p1,
    (select distinct ST_Transform(ST_GeomFromText(‘POINT (‘|| x || ‘ ‘ || y || ‘)’,22175),4326) as point2 from zonas_usuales where nombre = $1) as p2
    where ST_Equals(point1,point2) ‘t’ and dossegmentosmascortos(point1,point2,$1)) as tabla
    );
    END;
    $$ language ‘plpgsql’;