Validación¶

Note

Fecha	Autores
1 Noviembre 2012	Micho García (micho.garcia@geomati.co)
15 Octubre 2013	Jorge Arévalo(jorge.arevalo@geomati.co)

Excepto donde quede reflejado de otra manera, la presente documentación se halla bajo licencia : Creative Commons (Creative Commons - Attribution - Share Alike: http://creativecommons.org/licenses/by-sa/3.0/deed.es)

Validar geometrías¶

Una operación común cuando se trabaja con datos vectoriales es validar que dichos datos cumplen ciertas condiciones que los hacen óptimos para realizar análisis espacial sobre los mismos. O de otra forma, que cumplen ciertas condiciones topológicas.

Los puntos y las líneas son objetos muy sencillos. Intuitivamente, podemos afirmar que no hay manera de que sean topológicamente inválidos. Pero un polígono es un objeto más complejo, y debería cumplir ciertas condiciones. Y debe cumplirlas porque muchos algoritmos espaciales son capaces de ejecutarse rápidamente gracias a que asumen una consistencias de los datos de entrada. Si tuviéramos que forzar a que esos algoritmos revisaran las entradas, serían mucho más lentos.

Veamos un ejemplo de porqué esto es importante. Supongamos que tenemos este polígono sencillo:

# POLYGON((0 0, 0 1, 2 1, 2 2, 1 2, 1 0, 0 0));

Gráficamente:

Podemos ver el límite exterior de esta figura como un símbolo de infinito cuadrado. O sea, que tiene un lazo en el medio (una intersección consigo mismo). Si quisiéramos calcular el área de esta figura, podemos ver intuitivamente que tiene 2 unidades de área (si hablamos de metros, serían 2 metros cuadrados).

Veamos qué piensa PostGIS del área de esta figura:

# SELECT ST_Area(ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))'));

El resultado será:

# st_area
 ---------
0

¿Qué es lo que ha sucedido aquí?

El algoritmo de cálculo de áreas de PostGIS (muy rápido) asume que los anillos no van a intersectar consigo mismos. Un anillo que cumpla las condiciones adecuadas para el análisis espacial, debe tener el área que encierra siempre en el mismo lado. Sin embargo, en la imagen mostrada, el anillo tiene, en una parte, el área encerrada en el lado izquierdo. Y en la otra, el área está encerrada en el lado derecho. Esto causa que las áreas calculadas para cada parte del polígono tengan valores opuestos (1 y -1) y se anulen entre si.

Este ejemplo es muy sencillo, porque podemos ver rápidamente que el polígono es inválido, al contener una intersección consigo mismo (algo que ESRI permite en un SHP, pero PostGIS no, porque implementa SFSQL: http://www.opengeospatial.org/standards/sfs). Pero, ¿qué sucede si tenemos millones de polígonos? Necesitamos una manera de detectar si son válidos o inválidos. Afortunadamente, PostGIS tiene una función para esto: ST_IsValid, que devuelve TRUE o FALSE:

# SELECT ST_IsValid(ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))'))

Devuelve:

# st_isvalid
 ------------
         f

Incluso tenemos una función que nos dice la razón por la que una geometría es inválida:

# SELECT ST_IsValidReason(ST_GeometryFromText('POLYGON((0 0, 0 1, 1 1, 2 1, 2 2, 1 2, 1 1, 1 0, 0 0))'));

Que devuelve:

# st_isvalidreason
------------------------
Self-intersection[1 1]

Práctica¶

Vamos a comprobar la validez de las geometrías del shapefile world_borders:

# SELECT gid, name, ST_IsValidReason(geom) FROM tm_world_borders WHERE ST_IsValid(geom)=false;

Obtenemos el resultado:

#  gid |  name  |                  st_isvalidreason
  -----+--------+-----------------------------------------------------
    24 | Canada | Ring Self-intersection[-53.756367 48.5032620000001]
    33 | Chile  | Ring Self-intersection[-70.917236 -54.708618]
   155 | Norway | Ring Self-intersection[5.33694400000002 61.592773]
   175 | Russia | Ring Self-intersection[143.661926 49.31221]

Observamos que hay 4 polígonos con intersecciones consigo mismos. Esto es un ejemplo del aspecto que tienen estas auto-intersecciones:

Para resolver estos errores topológicos, tenemos a nuestra disposición la función ST_MakeValid. Esta función es nueva en PostGIS 2.0. Hasta entonces, estos problemas se resolvían con técnicas como hacer un buffer de tamaño 0 alrededor de la geometría inválida, y dejar que la función ST_Buffer la arreglara. Esto es así porque ST_Buffer en realidad construye una nueva geometría réplica de la antigua y construyendo un buffer alrededor de ella. Si este buffer es de tamaño 0, el resultado es solo la réplica de la anterior geometría. Pero al ser construida siguiendo las reglas topológicas de OGC, solucionaba muchos problemas como éste.

La función ST_MakeValid es más apropiada para arreglar geometrías. Únicamente requiere GEOS 3.3.0 o superior para funcionar (GEOS 3.3.4) si estamos usando PostGIS 2.1). Para saber qué versión de GEOS tenemos instalada basta con ejecutar:

# SELECT postgis_full_version()

Si se tiene una versión de GEOS inferior a la 3.3.0, se pueden seguir los consejos de Paul Ramsey: http://blog.opengeo.org/2010/09/08/tips-for-the-postgis-power-user/

Para comprobar el funcionamiento de ST_MakeValid vamos a crear una tabla nueva donde almacenemos únicamente uno de los polígonos conflictivos, marcado como erroneo. A continuación, crearemos un nuevo registro en dicha tabla con el polígono corregido.

Para hacerlo, ejecutemos esta query, que es algo compleja. Como sabemos que el problema es una auto-intersección que forma un anillo, vamos a desmontar la geometría en su lista de anillos y quedarnos solo con aquel que intersecta con el punto donde se detectó el error:

# SELECT * INTO invalid_geometries
FROM (
SELECT 'broken'::varchar(10) as status,
ST_GeometryN(geom, generate_series(1, ST_NRings(geom)))::geometry(Polygon,4326) as the_geom
FROM tm_world_borders
WHERE name = 'Chile') AS foo
WHERE ST_Intersects(the_geom, ST_SetSRID(ST_Point(-70.917236,-54.708618), 4326));

Con eso hemos creado la tabla invalid_geometries y añadido el anillo que contiene el error. Ahora añadamos un nuevo registro con el resultado de llamar a ST_MakeValid sobre el polígono erróneo:

# INSERT INTO invalid_geometries
VALUES ('repaired', (SELECT ST_MakeValid(the_geom) FROM invalid_geometries));

La función ST_MakeValid, realmente solo ha añadido un anillo más a la geometría inválida, para hacerla válida. Lo podemos comprobar con:

# SELECT status, ST_NRings(the_geom) FROM invalid_geometries;

Que devuelve:

# status  | st_nrings
----------+-----------
 broken   |         1
 repaired |         2

Ahora que ya hemos comprobado cómo funciona ST_MakeValid, podemos arreglar todas las geometrías inválidas:

# UPDATE tm_world_borders
SET the_geom = ST_MakeValid(the_geom)
WHERE ST_IsValid(the_geom) = false;

Una manera de evitar tener tablas con geometrías inválidas es definir una constraint que lo impida:

# ALTER TABLE tm_world_borders
ADD CONSTRAINT geometry_valid_check
CHECK (ST_IsValid(geom));