Wygładzanie / interpolacja rastra w Pythonie przy użyciu GDAL?

20

Rozwijam się w Pythonie i używam GDAL z OSGEO do manipulacji i interakcji z rastrami i plikami kształtów.

Chcę wziąć plik kształtu, który ma cechy punktowe i interpolować go do rastra powierzchniowego. Obecnie używam metody „RasterizeLayer”, która wypala wartość z funkcji punktowej do rastra (która jest ustawiona dla wszystkich wartości nodata), ale pozostawia wszystkie nietknięte piksele jako wartość „nodata”. Pozostaje mi więc raster typu szachownica.

Co mam po użyciu RasterizeLayer:

[Raster z używania gdal.rasterizelayer]

Co chcę dla produktu końcowego:

Uważam, że funkcja, której szukam, jest znana jako „Spline_sa ()” z importu arcgisscripting.

Czy GDAL ma podobną funkcję, czy może istnieje inna metoda uzyskania pożądanego wyniku?

— Doug
źródło

18

Rzuciłbym okiem na NumPy i Scipy - istnieje dobry przykład interpolacji danych punktów w SciPy Cookbook przy użyciu funkcji scipy.interpolate.griddata . Oczywiście wymaga to posiadania danych w tablicy numpy;

Za pomocą powiązań GDAL w Pythonie możesz wczytywać dane do Pythona za gdal.Dataset.ReadAsArray()pomocą rastra.
Za pomocą OGR można przechodzić przez warstwę funkcji i wyodrębniać dane punktowe z pliku shapefile (lub jeszcze lepiej, zapisać plik shapefile do pliku CSV za pomocą GEOMETRY=AS_XYZ[patrz format pliku OGR CSV] i wczytać plik csv do Pythona).

Gdy uzyskasz wynik w postaci siatki, możesz następnie użyć GDAL do zapisania wynikowej tablicy numpy do rastra.

Na koniec, jeśli nie masz szczęścia z biblioteką interpolacji Scipy, zawsze możesz spróbować również scipy.ndimage .

— om_henners
źródło

Dzięki za pomoc! Podaję Scipy.interpolate.griddata podejście wirowe. Prześlę z powrotem moje wyniki.

— Doug

1

Przepraszam, że tyle czasu zajęło mi powrót do tego postu. Powyższa odpowiedź jest w zasadzie tym, co zrobiłem, aby rozwiązać mój problem. Użyłem biblioteki interpolacji Scipy, aby wypełnić te przestrzenie nodata, a następnie zapisałem ją z powrotem do pasma rastrowego. Dzięki za pomoc chłopaki!

— Doug

@Doug Bez obaw - chętnie leczy!

— om_henners

1

Jak szybkie jest to rozwiązanie? Czy można go stosować do siatki 10k x 10k, gdzie znana jest tylko każda wartość 100x100? Próbowałem gdal_fillnodata, który jest niesamowicie szybki w porównaniu do dowolnej interpolacji, ale nie działa dobrze w przypadku zbyt rzadkich punktów. W tej chwili używam triangulacji z Sagi, ale jest ona bardzo wolna dla średnich tablic i kończy się niepowodzeniem dla dużych.

— Miro

12

Zobacz interfejs API gridowania GDAL . Nie wiem, czy jest to ujawnione w powiązaniach Pythona, ale jeśli nie, wywołujesz narzędzie gdal_grid za pośrednictwem modułu podprocesu .

Interfejs API siatki GDAL używa tylko odwrotnego ważenia odległości, średniej ruchomej i najbliższego sąsiada, nie implementuje splajnów. Inną opcją jest użycie Scipy .

— użytkownik2856
źródło

1

Trochę stary w tym wątku, ale napisałem prosty moduł, który wykorzystuje algorytm KNN ze sklearn o nazwie skspatial.

https://github.com/rosskush/skspatial

Możesz zaimportować plik kształtu za pomocą geopand i wybrać kolumnę, która interpoluje powierzchnię, którą można wyeksportować do rastra. Jest to bardzo prosty i prawdopodobnie nie najlepszy sposób, aby to zrobić, ale utrzymuje wszystko w czystości co najmniej python.

— rosskush
źródło