Suavizando / interpolando raster em Python usando GDAL?

20

Estou desenvolvendo em Python e usando o GDAL da OSGEO para manipular e interagir com rasters e shapefiles.

Eu quero pegar um shapefile com recursos de ponto e interpolar em um raster de superfície. No momento, estou usando o método 'RasterizeLayer', que grava um valor do recurso de ponto na varredura (que é definida com todos os valores de nodata), mas deixa todos os pixels intocados como valor de 'nodata'. Portanto, sou deixado com uma varredura do tipo quadriculado.

O que tenho depois de usar o RasterizeLayer:

[Varredura do uso de gdal.rasterizelayer]

O que eu quero para um produto final:

insira a descrição da imagem aqui

Acredito que a função que estou procurando seja conhecida como 'Spline_sa ()' na importação de scripts de arco.

O GDAL tem uma função semelhante ou existe um método diferente para obter a saída desejada?

python raster gdal interpolation Doug
fonte

Respostas:

18

Eu daria uma olhada no NumPy e Scipy - há um bom exemplo de interpolação de dados de pontos no SciPy Cookbook usando a função scipy.interpolate.griddata . Obviamente, isso requer que você tenha os dados em uma matriz numpy;

  • Usando as ligações python GDAL, você pode ler seus dados no Python usando gdal.Dataset.ReadAsArray()uma varredura.
  • Com o OGR, você percorreria a camada de recursos e extrairia os dados do ponto do shapefile (ou, melhor ainda, grave o shapefile em um CSV usando GEOMETRY=AS_XYZ[consulte o formato de arquivo OGR CSV] e leia o csv no Python).

Depois de obter uma saída em grade, você pode usar o GDAL para gravar a matriz numpy resultante em uma varredura.

Por fim, se você não tiver sorte com a biblioteca de interpolação Scipy, sempre poderá tentar o scipy.ndimage .

om_henners
fonte
Obrigado pela ajuda! Estou dando uma olhada na abordagem Scipy.interpolate.griddata. Vou postar de volta meus resultados.
Doug
1
Peço desculpas por demorar tanto para voltar a este post. A resposta acima é basicamente o que eu fiz para resolver meu problema. Usei a biblioteca de interpolação Scipy para preencher esses espaços nodata e depois escrevi de volta para a banda raster. Obrigado pela ajuda pessoal!
Doug
@ Doug Não se preocupe - feliz em healp!
Om_henners 13/03/12
1
Qual é a velocidade desta solução? Pode ser usado para uma grade de 10k x 10k, onde apenas cada valor de 100x100 é conhecido? Eu tentei o gdal_fillnodata, que é incrivelmente rápido em comparação com qualquer interpolação, mas não funciona bem para pontos muito esparsos. No momento, estou usando a triangulação da Saga, mas é muito lento para matrizes médias e falha com as grandes.
Miro
12

Veja a API de grade da GDAL . Não sei se isso está exposto nas ligações do Python, mas se não, você chama o utilitário gdal_grid por meio do módulo de subprocesso .

A API da grade GDAL usa apenas ponderação de distância inversa, média móvel e vizinho mais próximo; ela não implementa splines. Outra opção é usar o Scipy .

user2856
fonte
1

Um pouco antigo para este segmento, mas escrevi um módulo simples que usa o algoritmo KNN do sklearn chamado skspatial.

https://github.com/rosskush/skspatial

Você pode importar um arquivo de forma usando geopandas e selecionar uma coluna, que interpola uma superfície que pode ser exportada para uma varredura. É muito básico e provavelmente não é a melhor maneira de fazê-lo, mas mantém pelo menos tudo puro python.

rosskush
fonte