Divide un archivo shape complejo en una cuadrícula

Tengo un archivo de forma decentemente detallado con características de polígono / multipolígono (el archivo es de aproximadamente 500 mb). En realidad, es un archivo de forma del mundo entero, con las características que representan las costas. Necesito dividir estos datos usando una cuadrícula. Para ser claros, 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 encontré no funcionaron para mí.

He intentado:

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

• Usando las herramientas de intersección de OGR. Lo probé a través de ogr2ogr e incluso haciendo rodar mi propia herramienta C ++. Ambos tienen el mismo problema que QGIS. Tampoco exhiben este problema para archivos simples, pero fallan los más complejos. Como referencia, estoy usando un archivo shape de Australia y Nueva Zelanda, de menos de 20mb de tamaño, y tanto QGIS como OGR no logran 'cuadricularlo'.

Alguien sugirió usar PostGIS en un punto, ya que tiene una función de intersección, pero ST_Intersect de PostGIS usa el mismo back end GEOS que OGR. De hecho, ambos llaman a la misma función por lo que puedo decir, por lo que no creo que PostGIS arroje resultados diferentes.

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

EDITAR: Agregar más información

En respuesta a Simbamangu:

• El archivo de forma es básicamente información de la costa de OpenStreetMap. Es una versión fusionada del archivo 'procesado_p' (por lo que no está dividido en mosaicos) que obtuve enviando por correo electrónico su lista de desarrolladores. Tenga en cuenta que su división de mosaicos (en trozos de 100 km x 100 km con superposición) no es necesariamente lo que quiero: no quiero superponer y quiero la libertad de elegir el tamaño de la cuadrícula, o simplemente usaría el por defecto procesado_p.

• Por defecto, los datos de la costa tienen errores de geometría informados por QGIS. Arreglo estos errores con una pequeña herramienta que reuní 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 costa: https://github.com/tudelft-gist/prepair ). Repasar los archivos con esta herramienta corrige prácticamente todos los errores que QGIS detecta. Solo intento hacer la intersección después de limpiar los archivos.

• Exactamente lo que hice con QGIS: abra los datos para asegurarse de que se vean bien en QGIS. Intente dividirlo en mosaicos creando una capa de mosaicos usando Vector Grid con un espaciado específico y luego intersectando las dos capas, no se vaya. Intente usar un conjunto de datos más pequeño: seleccione características en Oceanía (Aus, NZ) para probar un conjunto de datos más pequeño; este archivo de forma tiene un tamaño <20mb. 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 de C ++ que funciona en WKT, por lo que convierto el archivo de forma a WKT usando ogr2ogr, luego alimento el archivo de texto a 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 termina haciendo exactamente lo mismo que usar ogr2ogr directamente.

En respuesta a Brent:

1. Lo hace. Todo está en WGS84 Lat / Lon
2. Pensé que lo contrario es cierto: que para un conjunto dado de mosaicos de cuadrícula, tomaría más tiempo intersectar un multipolígono gigante en lugar de un conjunto de características fragmentadas que podrían estar más localizadas espacialmente en cada mosaico, pero esto es Una sugerencia interesante: lo intentaré e informaré.
3. No se mantienen campos de atributos durante el proceso, solo estoy interesado en la geometría.
4. No estoy seguro, pero creo que estás diciendo que debería seleccionar los polígonos que se superponen a un mosaico de cuadrícula dado y luego realizar la intersección. Esto es demasiado engorroso manualmente con QGIS. Mi herramienta ya lo hace hasta cierto punto con una casilla de verificación. Hay un poco de velocidad, pero el resultado final sigue siendo pobre y no notablemente diferente.
5. Esta no es una opción. En este momento estoy tratando de dividir los datos para que sea 1 grado lat x 1 grado lon, y estoy buscando una metodología general / robusta que funcione en todos los casos. Intenté aumentar el tamaño de la cuadrícula (es decir, 10x10) para ver si obtendría mejores resultados y no veo ninguna correlación entre el tamaño de la cuadrícula y la calidad de la salida.

Editar # 2:

He intentado jugar más con esto y, en general, parece que los resultados no son confiables tanto con GEOS como con QGIS (que usa fTools, no sé si eso a su vez usa GEOS nuevamente). 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 aún no es una solución). Aquí hay una captura de pantalla de una cuadrícula realmente espaciada que funcionó principalmente, pero falló parcialmente en un mosaico:

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

Acabo de crear mis propias herramientas para hacer esto.

Usé la biblioteca Clipper ( http://www.angusj.com/delphi/clipper.php ) junto con OGR para dividir mi configuración de datos. Algo a tener en cuenta es que realizar intersecciones ingenuamente con esta lib lleva mucho tiempo, por lo que en su lugar utilicé un enfoque de árbol cuádruple ... es decir, dividir en cuatro celdas de cuadrícula, dividir cada una de ellas en cuatro más, etc., hasta obtener la resolución deseada. Sin embargo, la biblioteca funciona muy bien, adjunto una captura de pantalla que muestra los resultados en el hemisferio oriental:

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

Aquí están las herramientas en caso de que alguien se encuentre con un problema similar en el futuro. Tenga en cuenta que están pirateados juntos como prueba de conceptos y probablemente no debería usarlos directamente (aunque podría 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

Definitivamente parece que tienes problemas de geometría. Es poco probable que pueda obtener resultados limpios de un archivo de entrada sucio, independientemente del software utilizado, a menos que primero aborde sus problemas de geometría. Una vez que solucione los problemas de geometría, puede intentar lo siguiente si todavía tiene problemas:

1) Asegúrese de que su dataset de cuadrícula tenga la misma proyección que su dataset de polígono mundial. Si no, recrearlo en la proyección adecuada.

2) Convierta todas las funciones en una sola parte, mucho más fácil de procesar

3) Elimine todos los campos extraños manteniendo solo el campo de identificación que le permitirá unir sus atributos después de que se haya realizado la intersección, nuevamente mucho más fácil de procesar

4) En lugar de intersectar todo el conjunto de datos de la cuadrícula con todo el conjunto de datos de polígonos del mundo, intente recorrer sus polígonos de cuadrícula, seleccionando los polígonos de intersección en su conjunto de datos del mundo y realizando un clip basado en su polígono de cuadrícula. Esto le permitirá aislar cualquier problema y, al final, puede fusionar los resultados para lograr su objetivo original.

5) Intenta usar polígonos de cuadrícula más grandes.

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