Tengo una trama con la que me gustaría hacer algunas interpolaciones de puntos. Aquí es donde estoy:
from osgeo import gdal
from numpy import array
# Read raster
source = gdal.Open('my_raster.tif')
nx, ny = source.RasterXSize, source.RasterYSize
gt = source.GetGeoTransform()
band_array = source.GetRasterBand(1).ReadAsArray()
# Close raster
source = None
# Compute mid-point grid spacings
ax = array([gt[0] + ix*gt[1] + gt[1]/2.0 for ix in range(nx)])
ay = array([gt[3] + iy*gt[5] + gt[5]/2.0 for iy in range(ny)])
Hasta ahora, he probado la función interp2d de SciPy :
from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
Sin embargo, aparece un error de memoria en mi sistema Windows de 32 bits con un ráster 317 × 301:
Traceback (most recent call last):
File "<interactive input>", line 1, in <module>
File "C:\Python25\Lib\site-packages\scipy\interpolate\interpolate.py", line 125, in __init__
self.tck = fitpack.bisplrep(self.x, self.y, self.z, kx=kx, ky=ky, s=0.)
File "C:\Python25\Lib\site-packages\scipy\interpolate\fitpack.py", line 873, in bisplrep
tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
MemoryError
Admito que tengo una confianza limitada en esta función SciPy, ya que los parámetros bounds_error
o fill_value
no funcionan según lo documentado. No veo por qué debería tener un error de memoria, ya que mi raster es 317 × 301, y el algoritmo bilineal no debería ser difícil.
¿Alguien ha encontrado un buen algoritmo de interpolación bilineal, preferiblemente en Python, posiblemente diseñado con NumPy? ¿Alguna pista o consejo?
(Nota: el algoritmo de interpolación vecino más cercano es fácil:
from numpy import argmin, NAN
def nearest_neighbor(px, py, no_data=NAN):
'''Nearest Neighbor point at (px, py) on band_array
example: nearest_neighbor(2790501.920, 6338905.159)'''
ix = int(round((px - (gt[0] + gt[1]/2.0))/gt[1]))
iy = int(round((py - (gt[3] + gt[5]/2.0))/gt[5]))
if (ix < 0) or (iy < 0) or (ix > nx - 1) or (iy > ny - 1):
return no_data
else:
return band_array[iy, ix]
... pero prefiero los métodos de interpolación bilineal)
fuente
MemoryError
porque NumPy intenta acceder más allá de tuband_array
? Debe verificarax
yay
.gt
.Respuestas:
Traduje la fórmula a continuación (de Wikipedia ) a Python-speak para obtener el siguiente algoritmo, que parece funcionar.
Tenga en cuenta que el resultado se devolverá con una precisión aparentemente mayor que los datos de origen, ya que se clasifica al
dtype('float64')
tipo de datos de NumPy . Puede usar el valor de retorno con.astype(band_array.dtype)
para hacer que el tipo de datos de salida sea el mismo que la matriz de entrada.fuente
Lo probé localmente con resultados similares, pero estoy en una plataforma de 64 bits, por lo que no alcanzó el límite de memoria. Quizás, en lugar de eso, intente interpolar piezas pequeñas de la matriz a la vez, como en este ejemplo .
También podría hacer esto con GDAL, desde la línea de comando sería:
Para hacer la operación equivalente en Python, use ReprojectImage () :
fuente
ReprojectImage
técnica incorporada de GDAL .He tenido el problema exacto en el pasado, y nunca lo resolví usando interpolate.interp2d. He tenido éxito usando scipy.ndimage.map_coordinates . Intenta lo siguiente:
scipy.ndimage.map_coordinates (band_array, [ax, ay]], order = 1)
Esto parece dar el mismo resultado que bilineal.
fuente
scipy.interpolate.interp2d () funciona bien con scipy más moderno. Creo que las versiones anteriores asumen cuadrículas irregulares y no aprovechan las cuadrículas regulares. Me sale el mismo error que con scipy. versión = 0.11.0, pero en scipy. versión = 0.14.0, felizmente funciona en algunos resultados del modelo 1600x1600.
Gracias por las sugerencias en su pregunta.
fuente