Como usar el dialecto SQLite en OGR
4 de Abril de 2014
Más uso OGR, más me parece genial. Hace unos días tenia que escribir para un cliente un script Python que extraía datos texto desde una base de datos, creaba capas en formato shapefile, hacía una unión sobre los objetos de estas capas y luego creaba un mapa.
La parte más difícil fue la unión de los objetos. Empecé utilizando el módulo OGR y luego probé con Shapely (que es muy bueno). El resultado era más o menos satisfactorio pero era muy complicado por algo que me parecía bastante simple. Buscando un poco más, encontré la opción del dialecto SQLite de OGR. Basicamente, permite utilizar la sintaxis y las funciones de SQLite (y más importante, de Spatialite) dentro del comando ogr2ogr con cualquier formato de OGR.
Entonces podía usar la función st_union de Spatialite. A continuación la sintaxis:
| ogr2ogr capa_union.shp capa.shp -dialect sqlite -sql "SELECT st_union(st_buffer(capa.geometry, 0.0001)) FROM capa" -explodecollections |
capa.shp es el shape inicial. capa_union.shp es la capa resultado. Antes de usar la función st_union, hago un buffer (con la función st_buffer) casi nulo (de 0.0001 metro), eso permite que la unión funcione aún si hay errores de geometría en capa.shp.
Después de la sintaxis SQL, pusé la opción explodecollections. Eso es porque el resultado de st_union nos da un solo objeto final. Y Spatialite al contrario de PostGIS no tiene una función St_Dump. explodecollections permite convertir los objetos multiparts en singleparts. Es decir un objeto para cada polígono independiente.
La imagen siguiente muestra las 2 capas (capa.shp y abajo capa_union.shp):
Para usar este comando en Python solamente tengo que ponerlo como parametro (como texto, entonces entre ' ') de la función os.system:
| os.system('ogr2ogr capa_union.shp capa.shp -dialect sqlite -sql "SELECT st_union(st_buffer(capa.geometry, 0.0001)) FROM capa" -explodecollections') |
En termino de velocidad, usar ogr2ogr no es la mejor opción (para algunas de mis capas demora unos minutos) pero es tan simple a desarrollar que sería un pecado no usarlo.
En este ejemplo utilizamos una sola capa. Pero como hacemos cuando queremos aplicar una función SQLite/Spatialite sobre varias capas?
Por ejemplo, quiero obtener la diferencia entre 2 capas. Para eso utilizaremos la función symdifference. En la imagen a continuación, se ven estas 2 capas. Cada una contiene un solo objeto multiparts y ha sido obtenida con el mismo comando que arriba pero sin la opción explodecollections.
Las capas se llaman capa_union_mp.shp y capa_2_union_mp.shp (mp por multipart) y se encuentran en una carpeta llamada dif. La sintaxis es la siguiente:
| ogr2ogr dif/capa_dif.shp dif -dialect sqlite -sql "SELECT symdifference(capa_union_mp.geometry, capa_2_union_mp.geometry) FROM capa_union_mp, capa_2_union_mp" -explodecollections |
capa_dif.shp es la capa resultado. Para las capas iniciales, solo tengo que indicar la carpeta en la cual se encuentran (dif). Indico sus nombres en la parte de la sintaxis SQL. Decidí poner la opción explodecollections para tener objetos singleparts en vez de un solo objeto, pero obviamente es opcional y el resultado es visualmente lo mismo:
Como lo ven, usar el dialecto SQLite en OGR es bastante simple y nos permite aprovechar de la fuerza de Spatialite con la lista impresionante de sus funciones disponibles.