26 Sep 2011 Análisis espacial con PostGIS

En este post vamos a mostrar cómo utilizar varias funciones espaciales de PostGIS para obtener la línea de costa de una provincia mediterránea.  Los datos de partida constarán de  una única tabla espacial llamada ‘municipios’ . Esta tabla, con todos los municipios de la provincia, contiene los  siguientes atributos: código municipal, nombre del municipio  y geometría de cada municipio.

Fig 1. Datos de partida, municipios de Girona.

Para obtener la línea de costa vamos a utilizar, entre otras, dos funciones espaciales de postgis:  ST_Line_Substring y ST_Line_Locate_Point.

ST_Line_Substring devuelve un subconjunto de una linestring previamente conocida. Para ello requiere, además de la geometría inicial (linestring),  de una posición de inicio y de una posición final. Tanto la posición inicial como la posición final representan un valor real comprendido entre 0 y 1. De este modo podemos indicar que el subconjunto deseado va desde el 33% (0.33)  hasta el 66% (0.66) de la geometría inicial.

ST_Line_Locate_Point por su parte devuelve la posición que ocupa un punto concreto dentro de una geometría. El valor devuelto es un valor numérico y está comprendido entre cero (inicio) y uno (final)

Dada la geometría LinestringA, y dado el puntoA=(2,2) obtenemos el siguiente resultado.

Puntos de inicio y final de la línea de costa.

La línea de costa para el municipio de Girona alcanza, de norte a sur, desde el municipio de Portbou (codigo_municipio=17138)  hasta el municipio de Blanes (codigo_municipio=17023).

En la  imagen anterior podemos ver dónde empieza (punto A) y dónde termina (punto B) la línea de costa. Por lo tanto obteniendo esos puntos podremos aplicar la función ST_Line_Locate_Point para determinar en qué posición, respecto a la longitud total de la provincia, empieza y termina la línea de costa. Después solo faltará combinar el resultado de la función ST_Line_Locate_Point con ST_Line_Substring.

El comando SQL final tendrá el siguiente aspecto.

select st_line_substring(
geometria_perimetro_provincia ,
st_line_locate_point ( geometria_perimetro_provincia , puntoA ) ,
st_line_locate_point ( geometria_perimetro_provincia , puntoB )
)
)
from municipios

Obteniendo los puntos A y B.

Punto A.

La intersección del municipio con el perímetro de su bounding box nos devolverá una sola geometría de tipo MULTIPOINT que contendrá los 4 puntos que aparecen en  la siguiente imagen.  Con la función ST_Dump podemos descomponer esa geometría en 4 geometrías de tipo POINT.

El comando SQL quedaría como:

select (ST_Dump(
    ST_Intersection(st_boundary(the_geom), st_boundary(envelope(the_geom)) )
    )).geom as puntoA
from town_geom where code in (17138)

Finalmente podemos utilizar la consulta anterior y ordenarlos descendientemente por la componente ‘X’ y limitando el resultado a un solo punto.

select * from
(
select (ST_Dump(
        ST_Intersection(st_boundary(the_geom), st_boundary(envelope(the_geom)) )
        )).geom as puntoA
    from town_geom where code in (17138)
) as foo
order by x(puntoA) desc
limit 1

Punto B.

Podemos hacer ahora lo mismo para el punto B. La diferencia está tanto en el municipio (Blanes, 17023) como en la ordenación del punto. En este caso debemos ordenar por la componente  ‘Y’ de manera ascendente.

select * from
(
select (ST_Dump(
        ST_Intersection(st_boundary(the_geom), st_boundary(envelope(the_geom)) )
        )).geom as puntoB
    from town_geom where code in (17023)
) as foo
order by y(puntoB) asc
limit 1

Volviendo al comando SQL general…

En primer lugar vamos a recordar el esquema general del comando SQL que hemos visto al principio:

select st_line_substring(
    geometria_perimetro_provincia,
    st_line_locate_point(geometria_perimetro_provincia, puntoA),
    st_line_locate_point(geometria_perimetro_provincia, puntoB)
    )
)
from municipios

Por lo tanto, llegados a este punto, solo nos falta obtener el límite (linestring) de la geometría de la provincia. Para ello debemos unir las geometrías de los municipios con la función ST_Union y aplicar la función ST_Boundary:
st_boundary(st_union(the_geom))

Juntando todas las partes tenemos.

select asbinary( st_line_substring(
     st_boundary(st_union(the_geom))
    ,
    st_line_locate_point(
        st_boundary(st_union(the_geom))
        ,
        (
            select * from
            (
                select (ST_Dump(
                    ST_Intersection(st_boundary(the_geom), st_boundary(envelope(the_geom)) )
                    )).geom as puntoA
                from municipios where codigo_municipio in (17138)
            ) as foo
            order by x(puntoA) desc
            limit 1       
        )
    ),
    st_line_locate_point(
        st_boundary(st_union(the_geom))
        ,
        (
            select * from
            (
                select (ST_Dump(
                    ST_Intersection(st_boundary(the_geom), st_boundary(envelope(the_geom)) )
                    )).geom as puntoB
                from municipios where codigo_municipio in (17023)
            ) as foo
            order by y(puntoB) asc
            limit 1       
        )
    )
     )
)
from municipios

Obteniendo finalmente el siguiente resultado: