¿Suavizar/interpolar raster en Python utilizando GDAL?

Question

Estoy desarrollando en Python y utilizando GDAL de OSGEO para manipular e interactuar con rásters y shapefiles.

Quiero tomar un shapefile que tiene características puntuales e interpolarlas en un raster de superficie. En este momento estoy utilizando el método 'RasterizeLayer' que quema un valor de la característica de punto en el raster (que se establece con todos los valores nodata), pero deja todos los píxeles sin tocar como un valor 'nodata'. Por lo tanto, me queda un raster de tipo tablero de ajedrez.

Lo que tengo después de usar RasterizeLayer:

Lo que quiero para un producto final:

Creo que la función que estoy buscando se conoce como 'Spline_sa()' de la importación de arcgisscripting.

¿Dispone GDAL de una función similar o existe algún otro método para obtener el resultado deseado?

Answer 1

Yo le echaría un vistazo a NumPy y Scipy - hay un buen ejemplo de interpolación de datos puntuales en el archivo Libro de cocina SciPy utilizando el scipy.interpolate.griddata función. Obviamente esto requiere que tengas los datos en un array numpy;

• Utilizando los enlaces python de GDAL puede leer sus datos en Python utilizando gdal.Dataset.ReadAsArray() para una trama.
• Con OGR, se haría un bucle a través de la capa de características y se extraerían los datos de puntos del shapefile (o mejor aún, se escribiría el shapefile en un CSV utilizando GEOMETRY=AS_XYZ [véase el [Formato de fichero OGR CSV].](http://www.gdal.org/ogr/drv_csv.html) y leer el csv en Python).

Una vez que tengas una salida cuadriculada, puedes usar GDAL para escribir el array numpy resultante en un raster.

Por último, si no tienes suerte con la biblioteca de interpolación Scipy, siempre puedes probar con scipy.ndimage también.

Answer 2

Eche un vistazo al API de cuadrícula GDAL . No sé si eso está expuesto en los enlaces de Python, pero si no es así, puedes llamar a la aplicación gdal_grid mediante la función subproceso módulo.

La API de rejilla de GDAL sólo utiliza la ponderación de distancia inversa, la media móvil y el vecino más próximo, no implementa splines. Otra opción es utilizar Scipy .

Answer 3

Un poco viejo a este hilo pero he escrito un módulo simple que utiliza el algoritmo KNN de sklearn llamado skspatial.

https://github.com/rosskush/skspatial

Puedes importar un shapefile usando geopandas y seleccionar una columna e interpolará una superficie que puede ser exportada a un raster. Es muy básico y probablemente no es la mejor manera de hacerlo, pero mantiene todo python puro al menos.