Divida um shapefile complexo em uma grade

Eu tenho um shapefile decentemente detalhado com recursos de polígono / multipolígono (o arquivo tem cerca de 500mb). Na verdade, é um arquivo de forma do mundo inteiro, com os recursos representando as costas. Eu preciso dividir esses dados usando uma grade. Para ser claro, não quero 'classificar' os dados, mas, na verdade, cortar os polígonos em blocos. Sei que essa pergunta já foi feita antes, mas as soluções que encontrei não funcionaram para mim.

Eu tentei:

• Usando o QGIS e cruzando o conteúdo do meu shapefile com uma grade vetorial - os resultados são terríveis. A maior parte da massa terrestre principal desaparece magicamente, embora pareça que pedaços menores de terra às vezes o fazem. Devo observar que esse método funciona muito bem com dados muito mais simples (ou seja, menos pontos)

• Usando as ferramentas de interseção do OGR. Eu tentei através do ogr2ogr e até rodando minha própria ferramenta C ++. Ambos têm o mesmo problema que o QGIS. Eles também não apresentam esse problema para arquivos simples, mas falham nos mais complexos. Para referência, estou usando um shapefile da Austrália e da Nova Zelândia, com menos de 20mb de tamanho, e o QGIS e o OGR falham em 'gridificá-lo'.

Alguém sugeriu o uso do PostGIS em um ponto, pois ele possui uma função de interseção - mas o ST_Intersect do PostGIS usa o mesmo back-end do GEOS que o OGR. Na verdade, ambos chamam a mesma função, tanto quanto eu posso dizer, então não acho que o PostGIS traga resultados diferentes.

Eu estava procurando sugestões sobre o que mais eu poderia tentar. Preciso de um aplicativo ou kit de ferramentas robusto que possa dividir arquivos de forma altamente detalhados em blocos.

EDIT: Adicionando mais algumas informações

Em resposta a Simbamangu:

• O shapefile é basicamente o litoral do OpenStreetMap. É uma versão mesclada do arquivo 'processing_p' (para que não seja dividida em blocos) que recebi por e-mail na lista de desenvolvedores. Observe que a divisão de ladrilhos (em pedaços de 100 km x 100 km com sobreposição) não é necessariamente o que eu quero - eu não quero sobreposição e quero a liberdade de escolher o tamanho da grade ou apenas usaria o processado_p padrão.

• Por padrão, os dados da costa apresentam erros de geometria relatados pelo QGIS. Corrijo esses erros com uma pequena ferramenta que montei usando algum código que encontrei projetado para solucionar especificamente esse problema (reparando erros de geometria nos dados da costa: https://github.com/tudelft-gist/prepair ). A execução dos arquivos com esta ferramenta corrige praticamente todos os erros que o QGIS detecta. Eu só tento fazer o cruzamento depois de limpar os arquivos.

• Exatamente o que fiz usando o QGIS: abra os dados para garantir que eles estejam bem no QGIS. Tente dividi-lo em ladrilhos, criando uma camada de ladrilhos usando a Grade de vetores com um espaçamento especificado e, em seguida, cruzando as duas camadas - não é possível. Tente usar um conjunto de dados menor - selecione os recursos na Oceania (Aus, NZ) para tentar um conjunto de dados menor - esse arquivo de forma tem um tamanho <20mb. Mais uma vez, tente dividi-lo, não funciona.

• O que fiz com o OGR: ogr2ogr diretamente usando as opções '-spat' e '-clipsrc' com spat_extent. Também escrevi uma pequena ferramenta C ++ que funciona no WKT, então eu converto o shapefile em WKT usando ogr2ogr e, em seguida, alimento o arquivo de texto no meu aplicativo. Ele percorre o arquivo e chama o método Intersection () documentado aqui: http://www.gdal.org/ogr/classOGRGeometry.html . Eu acho que acaba fazendo exatamente a mesma coisa que usar ogr2ogr diretamente.

Em resposta a Brent:

1. Faz. Tudo está em WGS84 Lat / Lon
2. Eu pensaria que o oposto é verdadeiro - que, para um determinado conjunto de blocos de grade, levaria muito mais tempo para interceptar um multipolígono gigante em vez de um monte de recursos fragmentados que poderiam ser mais espacialmente localizados em cada bloco, mas isso é uma sugestão interessante - vou tentar e relatar de volta.
3. Nenhum campo de atributo é mantido durante o processo, estou interessado apenas em geometria.
4. Não tenho certeza, mas acho que você está dizendo que devo selecionar os polígonos que se sobrepõem a um determinado bloco de grade e, em seguida, executar a interseção. Isso é muito complicado manualmente com o QGIS. Minha ferramenta já faz isso em certa medida com uma caixa de seleção delimitadora. Há um pouco de velocidade, mas o resultado final ainda é ruim e não é visivelmente diferente.
5. Esta não é uma opção. No momento, estou tentando dividir os dados para que sejam de 1 grau lat x 1 grau lon, e estou procurando uma metodologia geral / robusta que funcione com todos os casos. Tentei aumentar o tamanho da grade (ou seja, 10 x 10) para ver se eu obteria melhores resultados e não vejo correlação entre o tamanho da grade e a qualidade da saída.

Edição # 2:

Eu tentei brincar com isso mais e, em geral, parece que os resultados não são confiáveis, tanto usando o GEOS quanto com o QGIS (que usa fTools, não sei se, por sua vez, usa o GEOS novamente). Eu estava errado ao afirmar que o tamanho da grade não tem nada a ver com os resultados - quanto maior a grade, melhores os resultados (isso é bom saber, mas ainda não é uma solução). Aqui está uma captura de tela de uma grade realmente espaçada que funcionou principalmente, mas falhou parcialmente em um bloco:

A geometria está limpa - o QGIS mostra 0 erros com a ferramenta "Verificar validade". Não pretendo abordar esse problema passo a passo; verificar se certos recursos falharam ou não na interseção em um conjunto de dados tão grande quando não é visualmente aparente (e não será com blocos menores) não é prático.

Acabei de criar minhas próprias ferramentas para fazer isso.

Eu usei a biblioteca Clipper ( http://www.angusj.com/delphi/clipper.php ) junto com o OGR para dividir meus dados configurados. Algo a se notar é que a realização de interseções ingenuamente com essa biblioteca leva muito tempo, então, em vez disso, usei uma abordagem de quadtree ... ou seja, divida em quatro células da grade, divida cada uma delas em mais quatro, etc, até obter a resolução desejada. A lib funciona muito bem, eu anexei uma captura de tela mostrando os resultados no hemisfério oriental:

O resultado acima levou cerca de 4,5 horas em um processador de 1,33 GHz.

Aqui estão as ferramentas no caso de alguém encontrar um problema semelhante no futuro. Observe que eles foram hackeados com provas de conceitos e você provavelmente não deve usá-los diretamente (pode servir como um bom ponto 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 você tem problemas de geometria. É improvável que seja capaz de produzir resultados limpos de um arquivo de entrada sujo, independentemente do software usado, a menos que você resolva primeiro os problemas de geometria. Depois de resolver os problemas de geometria, tente o seguinte se ainda estiver com problemas:

1) Verifique se o conjunto de dados da grade tem a mesma projeção que o conjunto de dados do polígono mundial. Caso contrário, recrie-o na projeção adequada.

2) Converta todos os recursos em peça única - muito mais fácil de processar

3) Remova todos os campos estranhos mantendo apenas o campo id, que permitirá que você junte seus atributos novamente após a interseção ter sido executada - novamente muito mais fácil de processar

4) Em vez de cruzar todo o conjunto de dados da grade com o conjunto de dados do polígono mundial, tente fazer um loop sobre os polígonos da grade, selecionando os polígonos que se cruzam no conjunto de dados do mundo e executando um clipe com base no polígono da grade. Isso permitirá que você isole qualquer problema e, no final, poderá mesclar os resultados para atingir seu objetivo original.

5) Tente usar polígonos de grade maiores.

Outra abordagem pode ter sido tentar uma conversão de vetor para varredura para criar um conjunto de dados de pontos e, em seguida, usar o conjunto de dados de pontos como base para escrever algum código para criar seus blocos.