Alat smoothing / generalisasi raster apa yang tersedia?

46

Saya memiliki DEM yang ingin saya menghaluskan atau menggeneralisasi untuk menghilangkan ekstrem topografi (memotong puncak dan mengisi lembah). Idealnya, saya juga ingin memiliki kendali atas jari-jari atau tingkat "keburaman". Pada akhirnya, saya akan membutuhkan seperangkat raster yang berkisar dari sedikit buram hingga benar-benar buram. (Secara teoritis, blurriest akan menjadi raster konstan rata-rata aritmatika dari semua nilai).

Apakah ada alat atau metode yang dapat saya gunakan (berdasarkan Esri, GDAL, GRASS)? Apakah saya perlu pulang memanggang rutinitas Gaussian blur saya sendiri ? Bisakah saya menggunakan filter low-pass (mis. Filter ArcGIS ), dan jika demikian, apakah saya perlu menjalankannya beberapa kali untuk mendapatkan efek radius besar?

Mike T
sumber
Bagaimana dengan hanya mengekspor raster ke ukuran sel yang lebih besar? Bukankah ini juga berakibat bungkam ekstrem?
1
Ya, itu juga akan mengurangi ekstrem (dengan asumsi bahwa resampling implisit melibatkan beberapa bentuk rata-rata) tetapi itu adalah cara yang mengerikan untuk menghaluskan DEM: Anda akan membuat sejumlah kecil blok besar. BTW, orang biasanya tidak perlu mengekspor raster untuk melakukan ini; agregasi serta resampling ke selsize yang berbeda adalah operasi dasar yang biasanya ditemukan dalam perangkat lunak berbasis raster.
Whuber

Jawaban:

29

Gaussian blur hanya fokus rata-rata tertimbang. Anda dapat membuatnya kembali dengan akurasi tinggi dengan urutan lingkungan lingkaran jarak pendek (tidak berbobot) berarti: ini adalah aplikasi dari Central Limit Theorem .

Anda punya banyak pilihan. "Filter" terlalu terbatas - hanya untuk lingkungan 3 x 3 - jadi jangan repot-repot. Pilihan terbaik untuk DEM besar adalah dengan mengambil perhitungan di luar ArcGIS ke dalam lingkungan yang menggunakan Fast Fourier Transforms: mereka melakukan perhitungan fokus yang sama tetapi (sebagai perbandingan) mereka melakukannya dengan sangat cepat. (GRASS memiliki modul FFT . Ini dimaksudkan untuk pemrosesan gambar tetapi Anda mungkin dapat menekannya untuk DEM Anda jika Anda dapat mengubah skala dengan presisi yang masuk akal ke dalam kisaran 0..255.) Kecuali itu, dua solusi setidaknya adalah layak dipertimbangkan:

  1. Buat satu set bobot lingkungan untuk memperkirakan keburaman Gaussian untuk lingkungan yang cukup besar. Gunakan lintasan blur berturut-turut ini untuk membuat urutan DEM yang lebih mulus.

    (Bobot dihitung sebagai exp (-d ^ 2 / (2r)) di mana d adalah jarak (dalam sel jika Anda suka) dan r adalah jari-jari efektif (juga dalam sel). Mereka harus dihitung dalam lingkaran yang memanjang setidaknya untuk 3r . Setelah melakukannya, bagilah masing-masing berat dengan jumlah mereka semua sehingga pada akhirnya jumlah mereka menjadi 1.)

  2. Atau, lupakan bebannya; jalankan rata-rata fokus lingkaran berulang kali. Saya telah melakukan ini untuk mempelajari bagaimana grid yang diturunkan (seperti kemiringan dan aspek) berubah dengan resolusi DEM.

Kedua metode akan bekerja dengan baik, dan setelah beberapa lintasan pertama akan ada sedikit untuk memilih antara keduanya, tetapi ada pengembalian yang semakin berkurang: jari-jari efektif n sarana fokus berturut-turut (semua menggunakan ukuran lingkungan yang sama) hanya (kira-kira) yang akar kuadrat dari n kali radius rata-rata fokus. Dengan demikian, untuk pengaburan dalam jumlah besar, Anda harus memulai lagi dengan lingkungan radius besar. Jika Anda menggunakan rata-rata fokus tertimbang, jalankan 5-6 melewati DEM. Jika Anda menggunakan bobot kira-kira Gaussian, Anda hanya perlu satu lintasan: tetapi Anda harus membuat matriks bobot.

Pendekatan ini memang memiliki rata-rata aritmatika DEM sebagai nilai pembatas.

whuber
sumber
1
Jika data Anda memiliki lonjakan, Anda dapat mencoba filter median ( en.wikipedia.org/wiki/Median_filter ) terlebih dahulu sebelum menerapkan blur yang lebih umum seperti yang disarankan oleh whuber.
MerseyViking
@ Hastings Itu saran yang bagus. Saya belum pernah melihat DEM dengan outlier lokal, tetapi sekali lagi saya tidak pernah harus memproses DEM mentah (seperti hasil LIDAR mentah) juga. Anda tidak dapat melakukan filter median dengan FFT, tetapi Anda hanya (biasanya) membutuhkan lingkungan 3 x 3 sehingga tetap merupakan operasi yang cepat.
whuber
Terimakasih Saya harus mengakui bahwa saya hanya pernah menggunakan data LiDAR yang sudah diproses sebelumnya, tetapi ada beberapa lonjakan signifikan dalam data SRTM yang akan mendapat manfaat dari filter median. Mereka cenderung menjadi 2 atau 3 sampel lebar, jadi filter median yang lebih besar akan diperlukan.
MerseyViking
@Mersey Anda masih ok dengan filter median yang lebih besar dari 5 x 5 atau 7 x 7. Jika Anda merenungkan (katakanlah) filter 101 x 101, bersiaplah untuk menunggu! Anda juga menyarankan poin penting yang perlu dijabarkan: adalah ide yang sangat bagus untuk melakukan analisis eksplorasi DEM sebelum melakukan apa pun. Ini akan termasuk mengidentifikasi lonjakan (outlier lokal) dan mengkarakterisasi ukuran dan luasannya. Anda ingin memastikan itu benar-benar artefak (dan bukan fenomena nyata) sebelum Anda memusnahkannya dengan filter!
whuber
1
+1 untuk FFT pada data ketinggian. Saya sudah benar-benar membuat itu bekerja di rumput untuk data NED 32bit untuk menghapus striping bi-directional. Pada akhirnya, ini juga bermasalah karena memperkenalkan kembali efek terasering yang mewabahi banyak DEM turunan kontur lainnya.
Jay Guarneri
43

Saya telah menjelajahi pendekatan signal.convolusi SciPy (berdasarkan buku masak ini ), dan saya mengalami beberapa keberhasilan yang sangat bagus dengan cuplikan berikut:

import numpy as np
from scipy.signal import fftconvolve

def gaussian_blur(in_array, size):
    # expand in_array to fit edge of kernel
    padded_array = np.pad(in_array, size, 'symmetric')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    return fftconvolve(padded_array, g, mode='valid')

Saya menggunakan ini dalam fungsi lain yang membaca / menulis float32 GeoTIFFs melalui GDAL (tidak perlu skala ulang ke 0-255 byte untuk pemrosesan gambar), dan saya telah menggunakan ukuran piksel yang mencoba (misalnya, 2, 5, 20) dan memiliki output yang sangat bagus (divisualisasikan dalam ArcGIS dengan piksel 1: 1 dan kisaran min / maks konstan):

DTM Gaussian

Catatan: jawaban ini telah diperbarui untuk menggunakan fungsi pemrosesan signal.fftconvoltve berbasis FFT jauh lebih cepat .

Mike T
sumber
1
+1 Solusi bagus! Saya tidak tahu pasti, tapi itu pertanda bagus bahwa signal.convolve menggunakan FFT.
whuber
Saya sedang mencari beberapa kode kabur untuk alat menjahit otomatis yang saya tulis dan menemukan ini. Kerja bagus @MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum Ingin mendengar lebih banyak tentang alat Anda. MikeToews Jawaban yang bagus dan cuplikan kode yang sangat dihargai.
Jay Laura
@JayLaura Tidak ada yang istimewa, hanya menulis alat untuk melakukan autostitch pada beberapa gambar yang saya ambil dengan beberapa teman dengan balon. Menggunakan kelas Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/…
Ragi Yaser Burhum
2
@whuber saat merevisi rutin ini, ia tidak menggunakan FFT, tapi sekarang, dan jauh lebih cepat.
Mike T
4

Ini bisa menjadi komentar untuk jawaban MikeT yang luar biasa , jika itu tidak terlalu panjang dan terlalu rumit. Saya sudah sering bermain dengannya dan membuat plugin QGIS bernama FFT Convolution Filters (dalam tahap "eksperimental") berdasarkan fungsinya. Selain menghaluskan, plugin juga dapat mempertajam tepi dengan mengurangi raster yang dihaluskan dari yang asli.

Saya sedikit meningkatkan fungsi Mike dalam proses:

def __gaussian_blur1d(self, in_array, size):
        #check validity
        try:
            if 0 in in_array.shape:
                raise Exception("Null array can't be processed!")
        except TypeError:
            raise Exception("Null array can't be processed!")
        # expand in_array to fit edge of kernel
        padded_array = np.pad(in_array, size, 'symmetric').astype(float)
        # build kernel
        x, y = np.mgrid[-size:size + 1, -size:size + 1]
        g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
        g = (g / g.sum()).astype(float)
        # do the Gaussian blur
        out_array = fftconvolve(padded_array, g, mode='valid')
        return out_array.astype(in_array.dtype)

Pemeriksaan validitas cukup jelas, tetapi yang penting adalah casting untuk mengapung dan kembali. Sebelum ini, fungsi membuat array integer hitam (nol saja), karena pembagian dengan jumlah nilai ( g / g.sum()).

Pavel V.
sumber
3

Di QGIS, saya mendapat hasil yang baik dengan mudah menggunakan penyaringan Gambar Orfeo Toolbox . Masuk akal dan mode batch berfungsi dengan baik. Difusi gaussian, rata-rata, atau anisotropik tersedia.

Catatan yang Radiusmengacu pada jumlah sel, bukan jarak.

Berikut adalah contoh penggunaan Smoothing (gaussian) :

  • Mentah:

    Tanpa filter

  • Tersaring:

    Saring

Tactopoda
sumber
1

Solusi yang bagus untuk animasi Gaussian blur dan keren. Mengenai alat Filter Esri yang disebutkan di atas, yang pada dasarnya hanya alat "Statistik Fokus" Esri yang dikodekan dengan ukuran 3x3. Alat Statistik Fokus memberi Anda lebih banyak opsi pada bentuk filter bergerak Anda, ukuran, dan statistik yang ingin Anda jalankan. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Anda juga dapat membuat filter "tidak teratur" tempat Anda memasukkan file teks Anda sendiri dengan bobot yang digunakan untuk setiap sel. File teks memiliki baris sebanyak yang Anda inginkan di area filter Anda, dengan nilai-nilai dibatasi spasi untuk kolom. Saya kira Anda harus selalu menggunakan jumlah baris dan kolom ganjil, sehingga sel target Anda ada di tengah.

Saya membuat spreadsheet excel untuk bermain dengan bobot berbeda yang baru saja saya salin / tempel ke file ini. Seharusnya mencapai hasil yang sama seperti di atas jika Anda menyesuaikan formula.

David A
sumber