Dividir un shapefile complejo en una cuadrícula

Question

Tengo un archivo shapefile decentemente detallado con características poligonales/multipoligonales (el archivo tiene unos 500mb). En realidad es un shapefile de todo el mundo, con las características que representan las líneas costeras. Necesito dividir estos datos mediante una cuadrícula. Para que quede claro, no quiero "ordenar" los datos, sino cortar los polígonos en mosaicos. Me doy cuenta de que esta pregunta se ha hecho antes, pero las soluciones que he encontrado no me han funcionado.

Lo he intentado:

• Utilizando QGIS e intersectando el contenido de mi shapefile con una cuadrícula vectorial, los resultados son terribles. La mayor parte de la masa terrestre principal desaparece por arte de magia, aunque parece que a veces aparecen trozos de tierra más pequeños. Debo señalar que este método funciona muy bien con datos mucho más simples (es decir, menos puntos)

• Utilizando las herramientas de Intersección de OGR. Lo he intentado tanto a través de ogr2ogr como rodando mi propia herramienta en C++. Ambos tienen el mismo problema que QGIS. Tampoco presentan este problema para archivos simples, pero fallan los más complejos. Como referencia, estoy usando un shapefile de Australia y Nueva Zelanda, de menos de 20mb de tamaño, y tanto QGIS como OGR fallan al 'cuadricularlo'.

Alguien sugirió utilizar PostGIS en un momento dado, ya que tiene una función de intersección - pero ST_Intersect de PostGIS utiliza el mismo back end de GEOS que OGR. De hecho, ambos llaman a la misma función por lo que puedo decir, así que no creo que PostGIS produzca resultados diferentes.

Estaba buscando sugerencias sobre qué más podría probar. Necesito una aplicación robusta o un conjunto de herramientas que pueda dividir archivos shapefiles muy detallados en mosaicos.

EDIT: Añadiendo algo más de información

En respuesta a Simbamangu:

• El shapefile es básicamente datos de la línea de costa de OpenStreetMap. Es una versión combinada del archivo 'processed_p' (por lo que no está dividido en mosaicos) que obtuve enviando un correo electrónico a su lista de desarrolladores. Tenga en cuenta que su división de los azulejos (en 100 kilometros x 100 kilometros trozos con la superposición) no es necesariamente lo que quiero - No quiero que la superposición, y quiero la libertad de elegir el tamaño de la cuadrícula, o yo sólo uso el processed_p por defecto.

• Por defecto los datos de la línea de costa tienen errores de geometría reportados por QGIS. Yo corrijo estos errores con una pequeña herramienta que armé usando un código que encontré diseñado para abordar específicamente este problema (reparación de errores de geometría en los datos de la línea de costa: https://github.com/tudelft-gist/prepair ). Repasando los archivos con esta herramienta se solucionan prácticamente todos los errores que detecta QGIS. Sólo intento hacer la intersección después de limpiar los archivos.

• Exactamente lo que hice usando QGIS: Abrir los datos para asegurarse de que se ven bien en QGIS. Intente dividirlos en mosaicos creando una capa de mosaicos usando Vector Grid con un espacio especificado, y luego intersectando las dos capas -- no va. Trate de usar un conjunto de datos más pequeño - seleccione las características en Oceanía (Aus, NZ) para probar un conjunto de datos más pequeño - este archivo de forma es < 20mb de tamaño. Nuevamente intente dividirlo, no funciona.

• Lo que hice con OGR: ogr2ogr directamente usando las opciones '-spat' y '-clipsrc' con spat_extent. También escribí una pequeña herramienta en C++ que funciona en WKT, así que convierto el shapefile a WKT usando ogr2ogr, y luego introduzco el archivo de texto en mi aplicación. Se ejecuta a través del archivo y llama al método Intersection() documentado aquí: http://www.gdal.org/ogr/classOGRGeometry.html . Creo que acaba haciendo exactamente lo mismo que usando ogr2ogr directamente.

En respuesta a Brent:

1. Así es. Todo está en WGS84 Lat/Lon
2. Yo habría pensado que lo contrario es cierto - que para un conjunto dado de azulejos de la cuadrícula, tomaría mucho más tiempo para intersecar un multipolígono gigante en lugar de un montón de características fragmentadas que podrían ser más localizados espacialmente a cada azulejo, pero esta es una sugerencia interesante - Voy a probar y reportar.
3. No se guarda ningún campo de atributos durante el proceso, sólo me interesa la geometría.
4. No estoy seguro, pero creo que estás diciendo que debo seleccionar los polígonos que se superponen a una determinada malla de la cuadrícula y luego realizar la intersección. Esto es demasiado engorroso manualmente con QGIS. Mi herramienta ya lo hace hasta cierto punto con una comprobación del cuadro delimitador. Hay un poco de velocidad, pero el resultado final sigue siendo pobre y no es notablemente diferente.
5. Esto no es una opción. En este momento estoy tratando de dividir los datos de manera que su 1 deg lat x 1 deg lon, y estoy buscando una metodología general / robusto que funciona con todos los casos. He probado a aumentar el tamaño de la cuadrícula (es decir, 10x10) para ver si obtenía mejores resultados y no veo ninguna correlación entre el tamaño de la cuadrícula y la calidad del resultado.

Edición #2:

He intentado jugar más con esto y en general parece que los resultados son poco fiables tanto con GEOS como con QGIS (que utiliza fTools, no sé si a su vez utiliza GEOS de nuevo). Me equivoqué al afirmar que el tamaño de la cuadrícula no tiene nada que ver con los resultados - cuanto más grande es la cuadrícula, mejores son los resultados (es bueno saberlo pero sigue sin ser una solución). Aquí hay una captura de pantalla de una cuadrícula muy espaciada que en su mayor parte funcionó, pero falló parcialmente en un azulejo:

La geometría está limpia -- QGIS muestra 0 errores con la herramienta "Check Validity". No quiero abordar este problema paso a paso; no es práctico verificar si ciertas características fallaron la intersección en un conjunto de datos tan grande cuando no es visualmente aparente (y no lo será con mosaicos más pequeños).

Answer 1

Acabé creando mis propias herramientas para hacerlo.

He utilizado la biblioteca Clipper ( http://www.angusj.com/delphi/clipper.php ) junto con OGR para dividir mi conjunto de datos. Algo que hay que tener en cuenta es que realizar las intersecciones de forma ingenua con esta librería lleva mucho tiempo, así que en su lugar utilicé un enfoque de árbol cuádruple... es decir, dividir en cuatro celdas de la cuadrícula, dividir cada una de ellas en cuatro más, etc, hasta obtener la resolución deseada. La librería funciona muy bien, he adjuntado una captura de pantalla que muestra los resultados en el hemisferio oriental:

El resultado anterior tardó unas 4,5 horas en un procesador de 1,33 GHz.

Aquí están las herramientas por si alguien se encuentra con un problema similar en el futuro. Por favor, ten en cuenta que son pruebas de concepto pirateadas y probablemente no deberías usarlas directamente (aunque podrían servir como un buen punto de partida para algo):

https://github.com/preet/scratch/tree/master/gis/polytoolkit

https://github.com/preet/scratch/tree/master/gis/shapefiles/shptk

Answer 2

Definitivamente parece que tienes problemas de geometría. Es poco probable que pueda obtener resultados limpios a partir de un archivo de entrada sucio, independientemente del software utilizado, a menos que primero resuelva sus problemas de geometría. Una vez que haya resuelto sus problemas de geometría, podría intentar lo siguiente si sigue teniendo problemas:

1) Asegúrese de que su conjunto de datos de cuadrícula tiene la misma proyección que su conjunto de datos de polígonos del mundo. Si no es así, recréalo en la proyección adecuada.

2) Convertir todas las características en una sola pieza - mucho más fácil de procesar

3) Elimine todos los campos extraños manteniendo sólo el campo id que le permitirá volver a unir sus atributos después de que se haya realizado la intersección - de nuevo mucho más fácil de procesar

4) En lugar de intersecar todo el conjunto de datos de la cuadrícula con todo el conjunto de datos de los polígonos del mundo, intente hacer un bucle sobre sus polígonos de la cuadrícula, seleccionando los polígonos que se intersecan en su conjunto de datos del mundo y realizando un recorte basado en su polígono de la cuadrícula. Esto le permitirá aislar cualquier problema y al final podrá fusionar los resultados para lograr su objetivo original.

5) Pruebe a utilizar polígonos de cuadrícula más grandes.

Answer 3

Otro enfoque podría haber sido intentar una conversión de vector a trama para crear un conjunto de datos de puntos y luego utilizar el conjunto de datos de puntos como base para escribir algún código para crear sus mosaicos.