Saya memiliki raster yang ingin saya lakukan interpolasi titik dengan. Di sinilah tempat saya:
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)])
Hingga kini, saya sudah mencoba fungsi interp2d SciPy :
from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
namun saya mendapatkan kesalahan memori pada sistem Windows 32-bit saya dengan raster 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
Saya akui, saya memiliki keyakinan terbatas pada fungsi SciPy ini, karena parameter bounds_error
atau fill_value
tidak berfungsi seperti yang didokumentasikan. Saya tidak melihat mengapa saya harus memiliki kesalahan memori, karena raster saya adalah 317 × 301, dan algoritma bilinear seharusnya tidak sulit.
Adakah yang menemukan algoritma interpolasi bilinear yang baik, lebih disukai dengan Python, mungkin disesuaikan dengan NumPy? Ada petunjuk atau saran?
(Catatan: algoritma interpolasi tetangga terdekat adalah kue mudah:
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]
... tapi saya lebih suka metode interpolasi bilinear)
sumber
MemoryError
karena NumPy mencoba mengakses di luar Andaband_array
? Anda harus memeriksaax
danay
.gt
.Jawaban:
Saya telah menerjemahkan rumus di bawah ini (dari Wikipedia ) ke dalam Python-speak untuk menghasilkan algoritma berikut, yang tampaknya berfungsi.
Perhatikan bahwa hasilnya akan dikembalikan dengan ketepatan yang lebih tinggi daripada data sumber, karena ini diklasifikasi ke
dtype('float64')
tipe data NumPy . Anda dapat menggunakan nilai kembali dengan.astype(band_array.dtype)
untuk membuat tipe data output sama dengan array input.sumber
Saya mencobanya secara lokal untuk hasil yang serupa, tetapi saya menggunakan platform 64-bit sehingga tidak mencapai batas memori. Mungkin alih-alih coba interpolasi potongan-potongan kecil array pada suatu waktu, seperti dalam contoh ini .
Anda juga bisa melakukan ini dengan GDAL, dari baris perintah:
Untuk melakukan operasi yang setara di Python, gunakan ReprojectImage () :
sumber
ReprojectImage
teknik bawaan GDAL .Saya memiliki masalah yang tepat di masa lalu, dan tidak pernah menyelesaikannya menggunakan interpolate.interp2d. Saya telah berhasil menggunakan scipy.ndimage.map_coordinates . Coba yang berikut ini:
scipy.ndimage.map_coordinates (band_array, [ax, ay]], order = 1)
Ini tampaknya memberikan hasil yang sama dengan bilinear.
sumber
scipy.interpolate.interp2d () berfungsi baik dengan scipy yang lebih modern. Saya pikir versi yang lebih lama menganggap grid tidak teratur dan tidak mengambil keuntungan dari grid biasa. Saya mendapatkan kesalahan yang sama seperti yang Anda lakukan dengan Scipy. versi = 0.11.0, tetapi pada scipy. versi = 0.14.0, itu dengan senang hati bekerja pada beberapa keluaran model 1600x1600.
Terima kasih atas petunjuknya dalam pertanyaan Anda.
sumber