Mam raster, z którym chciałbym wykonać interpolację punktową. Oto gdzie jestem:
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)])
Do tej pory wypróbowałem funkcję interp2d SciPy :
from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
jednak pojawia się błąd pamięci w moim 32-bitowym systemie Windows z rastrem 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
Przyznam, że mam ograniczone zaufanie do tej funkcji SciPy, ponieważ parametry bounds_error
lub fill_value
nie działają zgodnie z dokumentacją. Nie rozumiem, dlaczego powinienem mieć błąd pamięci, ponieważ mój raster ma wymiary 317 × 301, a algorytm dwuliniowy nie powinien być trudny.
Czy ktoś napotkał dobry algorytm interpolacji dwuliniowej, najlepiej w Pythonie, być może dostosowany do NumPy? Wszelkie wskazówki lub porady?
(Uwaga: algorytm interpolacji najbliższego sąsiada to łatwe ciasto:
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]
... ale zdecydowanie wolę metody interpolacji dwuliniowej)
źródło
MemoryError
ponieważ NumPy próbuje uzyskać dostęp poza twoimband_array
? Powinieneś sprawdzićax
iay
.gt
.Odpowiedzi:
Przetłumaczyłem poniższy wzór (z Wikipedii ) na Python-speak, aby uzyskać następujący algorytm, który wydaje się działać.
Zauważ, że wynik zostanie zwrócony z pozornie większą precyzją niż dane źródłowe, ponieważ jest on klasyfikowany w górę do
dtype('float64')
typu danych NumPy . Można użyć wartości zwracanej za pomocą,.astype(band_array.dtype)
aby typ danych wyjściowych był taki sam jak tablica wejściowa.źródło
Próbowałem lokalnie, aby uzyskać podobne wyniki, ale jestem na platformie 64-bitowej, więc nie osiągnęło limitu pamięci. Być może zamiast tego spróbuj interpolować małe fragmenty tablicy na raz, tak jak w tym przykładzie .
Możesz to również zrobić za pomocą GDAL, z wiersza poleceń byłoby to:
Aby wykonać równoważną operację w Pythonie, użyj ReprojectImage () :
źródło
ReprojectImage
techniki GDAL .W przeszłości miałem dokładny problem i nigdy nie rozwiązałem go za pomocą interpolate.interp2d. Odniosłem sukces, używając scipy.ndimage.map_coordinates . Spróbuj wykonać następujące czynności:
scipy.ndimage.map_coordinates (tablica_pasów, [ax, ay]], kolejność = 1)
Wydaje się, że daje to taką samą wydajność jak dwuliniowa.
źródło
scipy.interpolate.interp2d () działa dobrze z bardziej nowoczesnym scipy. Myślę, że starsze wersje zakładają nieregularne siatki i nie wykorzystują zwykłych siatek. Otrzymuję ten sam błąd, co ty z Scipy. wersja = 0.11.0, ale na scipy. wersja = 0.14.0, na szczęście działa na niektórych modelach o rozdzielczości 1600 x 1600.
Dziękuję za wskazówki w twoim pytaniu.
źródło