Tworzenie NMT z chmury punktów LAS w SAGA GIS

By móc pracować w nich z efektami skanowania LiDAR konieczne jest stworzenie na bazie chmur punktów (zazwyczaj zapisywanych obecnie w formacie LAS) rastrowych Numerycznych Modeli Terenu zarówno w postaci modeli powierzchni czy też pokrycia terenu, jak i rzeźby. Dla uzyskania tej drugiej opcji (co przecież jest największą zaletą skanowania LiDAR), konieczne jest odfiltrowanie roślinności i zabudowy.

Programów umożliwiających takie działania jest oczywiście sporo. Wiele z nich to bardzo drogie, komercyjne oprogramowanie. Dlatego skupimy się na darmowym, otwartoźródłowym SAGA GIS.

Krok 1: importujemy chmurę punktów LAS

Modules-File-Shapes-Import-Import LAS files

Istotne jest by zaznaczyć opcję „classification”. Pozostałe opcje wedle uznania – do stworzenia modelu rzeźby terenu nie będą potrzebne.

W ten sposób oprzemy się na istniejącej klasyfikacji chmury punktów. W przypadku plików z programu ISOK nie zawsze może być ona najlepsza dla celów archeologicznych, ale z moich prywatnych doświadczeń wynika, że zazwyczaj wystarcza.

Krok 2 (opcjonalny): podgląd chmury punktów.

Modules – Shapes – Point Cloud – Visualisation – Point Cloud Viewer

SAGA posiada moduł podglądu chmury punktów. Nie daje ona tylu możliwości co specjalistyczne programy obsługujące dane LiDAR, ale wystarczy do podejrzenia tego, co nam się udało zaimportować 😉

Krok 3: Odfiltrowanie roślinności i zabudowy.

Modules – Shapes – Point Clouds – Tools – Point Clouds Reclassifier/Subset Extractor

W oknie dialogowym musimy wybrać:

– w rubryce Point Cloud wybieramy chmurę, nad którą pracujemy.

– w rubryce Attribute: classification

– w rubryce Mode of operation: Extract Subset

– w rubryce old value możemy wpisać „2”

W efekcie otrzymamy chmurę punktów odwzorowujących ukształtowanie powierzchni. W przeglądarce chmury będzie wyglądała tak: Jak widać nieco więcej tutaj dziur, wynikających z odfiltrowania powierzchni.

Krok 4: griding – tworzenie rastra NMT.

Modules – Grid – Gridding – Interpolation from points – Triangulation

W oknie dialogowym:

– w rubryce Points wybieramy chmurę po odfiltrowaniu roślinności i budynków (domyślnie ma rozszerzenie nazwy subset_classification).

– w rubryce Atribute wybieramy „z” (czyli wysokość)

Po zatwierdzeniu pojawi się kolejne okno dialogowe:

Tutaj w rubryce Cellsize wpisujemy „1”. Po zatwierdzeniu rozpocznie się proces triangulacji, który może trwać nawet kilka godzin. Jest to oczywiście uzależnione od ilości punktów (dla terenów zalesionych trwa to krócej, bo wyeliminowanie roślinności zmniejsza liczbę punktów) i jakości komputera. Zwykle można w tym czasie zrobić sobie kawę. Albo nawet dwie… albo obiad…

Alternatywnie można zamiast triangulacji użyć funkcji:

Modules – Shapes – Point Clouds – Conversion – Poin Cloud to Grid

Zwykle trwa to nieco szybciej niż triangulacja.

Krok 5 – zapis uzyskanego rastra

Modules – File – Grid – Export – Export ESRI Arc/Info grid (o ile chcemy mieć plik z rozszerzeniem .asc czytany przez większość programów GIS, w tym Qgis). Opcjonalnie możemy importować do geotiff. Lub do dowolnego rastrowego formatu, jaki nam do głowy przyjdzie. Akurat te dwa pierwsze są najczęściej używane…

Niekiedy zamiast triangulacji warto zastosować inną metodę gridingu. W przypadku niektórych terenów zalesionych, po odfiltrowaniu roślinności, liczba punktów, które są podstawą obliczenia rastra NMT jest znacznie mniejsza od wyjściowej (n.p. najczęściej stosowanej w programie ISOK 4 pkt. na 1m kwadratowy). Po wykonaniu triangulacji efekt końcowy może być daleko niezadowalający, jak w tym przypadku:

Fragment NMT z zastosowaniem cieniowania (hillshading). Widać wyraźne, duże trójkąty, wynikające z triangulacji pomiędzy zbyt „luźno” rozłożonymi punktami.

W takiej sytuacji można zastosować np. metodę krigingu (ordinary kriging), która dobrze się sprawdza w podobnych sytuacjach (zob. dawny wpis o tworzeniu NMT z planów warstwicowych).

Efekt końcowy jest wyraźnie lepszy.

W skrajnych przypadkach może okazać się, że oparcie się na istniejącej klasyfikacji chmury nie będzie wystarczająco korzystne. Konieczna będzie reklasyfikacja chmury. Ale o tym może innym razem…

Reklama

SAGA GIS – łączenie kilku rastrowych NMT

Czasami potrzebujemy połączyć (zmozaikować) kilka Numerycznych Modeli Terenu i stworzyć z nich jeden większy. Połączony NMT ułatwia tworzenie wizualizacji, analiz czy nawet ujednoliconego sposobu wyświetlania (kolor, przezroczystość etc). Dla doświadczonych użytkowników programów GIS to nie jest wielki problem, dla początkujących może być. Dla nich właśnie krótki tutorial dla SAGA GIS, który ma tę zaletę, że doskonale radzi sobie z rastrami, a do tego jest nieco łatwiejszy w obsłudze niż GRASS.

1. Otwieramy pliki rastrowe: Modules – File – GDAL/OGR – GDAL: Import Raster:

W oknie dialogowym wybieramy raster, który chcemy otworzyć. Pamiętajmy, że w SAGA każdorazowo musimy kliknąć na ikonę otwartej warstwy w zakładce „Data”, by zaimportowana warstwa wyświetliła się w oknie mapy i znalazła swoje miejsce w drzewie warstw w zakładce „Maps”.

Tak mniej więcej obraz uzyskamy po zaimportowaniu wszystkich arkuszy NMT.

2. Drugim krokiem jest po prostu wejście w opcję: Modules – Grid – Grid System – Mosaicking

Pojawi się okno dialogowe, w którym klikamy na „Input Grids”

3. W drugim okienku dialogowym, które pojawi się po wybraniu opcji „Input Grids” wybieramy warstwy, które mają zostać połączone:

Po zatwierdzeniu w oknie dialogowym mozaikowania wybieramy następujące opcje:

Interpolation: B-Spline Interpolation

Overlapping Areas: mean

4. Po zatwierdzeniu powinniśmy uzyskać jednolity, połączony raster, który oczywiście musimy wywołać kliknięciem do okna mapy:

5. Gotowy raster możemy zapisać np. w postaci pliku Arc Grid (.asc): Modules – File – Grids – Export – Export ESRI Arc/Info Grid

W okienku dialogowym wybieramy nowo stworzoną mozaikę (zazwyczaj na końcu listy w zakładce Grid systems) oraz wybieramy nazwę i lokalizację dla gotowego pliku:

I gotowe.

NMT w 3D w SAGA GIS

Darek Bobak zauważył przy wpisie o Rozprzy, że GIS to nie tylko ładne mapy. To prawda. Ale nie zaszkodzi by przy okazji jakąś ładną mapę spłodzić. W niniejszym wpisie kilka słów poświęcę spłodzeniu estetycznie wyglądającego Numerycznego Modelu Terenu w ujęciu 3D (czy też pseudo 3D), czyli modnego wśród archeologów rzutu izometrycznego.

W tym zbożnym celu wykorzystuje się zwykle program Surfer. Jest to rzeczywiście wspaniałe oprogramowanie, niestety ma jedną podstawową wadę – jest drogie. No i przeznaczony jest wyłącznie na system Windows, co zniechęca takich nawiedzonych linuksiarzy jak ja (no dobra – tak na prawdę to nie mam kasy). Oczywiście wielu osób nie zniechęca cena Surfera. Właściwie nie ma ona żadnego znaczenia wobec dostępności wersji Pirate Edition. Tutaj jednak w ramach propagowania legalnego i otwartego oprogramowania spróbujemy zastąpić Surfera przez SAGA GIS.

1. Budujemy Numeryczny Model Terenu. W tym celu importujemy do programu punkty pomiarowe w formacie shapefile: File – GDAL/OGR – Import vector data.

Aby teraz wyświetlić zaimportowaną warstwę musimy w lewym oknie wejść w zakładkę „Data”.

I kliknąć dwukrotnie na miniaturę zaimportowanej warstwy. Od tej pory najlepiej będzie oglądać ją w zakładce „Maps” i podzakładce „Tree”, która wyświetla drzewo warstw GIS. Natomiast w oknie podglądu mapy będzie to wyglądało tak:

Do stworzenia modelu możemy użyć opcji triangulacji (Interpolation from Points – Triangulation) albo krigingu (Spatial and geostatistics – Kriging – Ordinary kriging). Jeśli wszystko pójdzie dobrze uzyskujemy mniej więcej taki efekt (tyle, że w kolorystyce red-blue):

2. Teraz należy dobrać odpowiedni sposób prezentacji naszego modelu. W domyślnie skonfigurowanym oknie programu SAGA, w prawym okienku możemy ustawić sobie sposób prezentacji. Wybieramy więc opcję: Colors – Scaling i w ramce klikamy na Colors. Pojawi się okno dialogowe, w którym wybieramy Presets. W Presets settings wybieramy schemat kolorystyczny dla naszej mapy. Niech będzie to „Topography”.

Należy pamiętać, że każdą zmianę w oknie ustawień musimy zatwierdzić klikając Apply.

W efekcie uzyskamy taki mniej więcej obraz:

3. Teraz zbudujemy sobie warstwę cieniowaną, która nada nieco głębi naszej mapie. Wybieramy Terrain Analysis – Lighting – Analytical Hillshading. W rubryczce Exaggeneration możemy ustawić przewyższenie, np. 3.

Azymut i kąt padania światła możemy zostawić bez zmian (ale z tą opcją warto poeksperymentować).

W efekcie powinniśmy otrzymać taki obrazek:

4. Kolejnym etapem jest ustawienie przezroczystości dla warstwy z cieniowaniem. Robimy to w prawym okienku (czyli oknie ustawień). Możemy ustawić np. 40. W efekcie nasz obrazek będzie wyglądał tak:

5. Teraz możemy wykonać projekcję 3D klikając na stosowną ikonkę na pasku górnym. W oknie ustawień widoku 3D (dostępnym w każdym momencie z 3D View – Properties) możemy ustawić przewyższenie (Exaggeneration) na wartość np. 3 (co podkreśli formę terenową). Warto też ustawić nieco większą niż domyślna wartość rozdzielczości (Resolution), przynajmniej na 800 (to w dużej mierze zależeć będzie od jakości waszej karty graficznej – im wyższa wartość tym lepszej karty będziecie potrzebowali; 2000 jest już całkiem nieźle i więcej nie trzeba zazwyczaj). W oknie 3D będzie to wyglądało tak:

Jak widać, coś już widać, ale dalekie to jest od jakości, której byśmy się spodziewali. Oczywiście widoczna na obrazku pikseloza w znacznej mierze zależeć będzie od jakości waszego modelu. Jeżeli jest to duży NMT, w dodatku oparty na pomiarach LiDAR, to piksele nie powinny przeszkadzać. Jeśli jest mniejszy, będzie bolało. Na szczęście jest na to rada.

6. Żeby pozbyć się wielkich pikseli wystarczy dalej pogrzebać w opcjach wyświetlania map, czyli w naszym prawym okienku programu. W ustawieniach Display wybieramy Interpolation i ustawiamy opcję B-Spline:

Jeśli zastosujemy ten zabieg dla obu naszych warstw rastrowych (NMT i warstwy cieniowanej), to nasz widok 3D powinien wyglądać tak:

Co jest już znacznie ładniejsze od wersji zamieszczonej wyżej.

7. Teraz możemy wzbogacić naszą mapę o warstwę izolinii (czyli warstwę warstwic). Wybieramy: Modules – Shapes – Grid – Vectorization – Contour Lines from Grid. W oknie dialogowym najważniejszym parametrem będzie wybranie odpowiedniego cięcia (Equidistance). Jeżeli potrzebujemy wartości dziesiętnej to musimy pamiętać by wpisać ją z kropką, a nie przecinkiem. Cięcie co 25 cm będziemy więc wpisywać: „0.25”. W efekcie powinniśmy uzyskać wektorową mapę warstwicową:

Możemy teraz zapisać nasz widok 3D w postaci obrazka. Wybieramy 3D View – Save as Image. Domyślne rozmiary zapisywanego obrazka mogą być nieco zbyt małe, więc możemy je nieco powiększyć, by uzyskać lepszą rozdzielczość. Niestety program nie zrobi tego automatycznie, więc musimy sobie samodzielnie wpisać wartości np. dwukrotnie większe niż domyślnie zaproponowane. Efekt końcowy wygląda tak:

Całość możemy wzbogacić o legendę. W menu „Map” klikamy „Copy Legend to Clipboard”. Legendę możemy wstawić na obrazek posługując się dowolnym programem graficznym, np. Gimpem. Otwieramy obrazek i wklejamy na niego legendę pozostającą w schowku systemowym. I gotowe. Nie będzie to Surfer, ale też będzie nieźle…

Grodzisko

W ostatnich tygodniach stosunkowo niewiele pojawia się tutaj nowych wpisów. Czuję się zobowiązany wytłumaczyć czytelnikom z tych zaniedbań:

Po pierwsze – mam sporo pracy wynikającej z przygotowywania kilku wystąpień konferencyjnych. Pracuję równolegle nad kilkoma tekstami dotyczącymi bardzo różnych tematów i obszarów od Drohiczyna na wschodzie, po zachodnie skraje Pomorza.

Po drugie, w międzyczasie próbowałem (bez powodzenia) ratować zabytki archeologiczne w Ostrowitem, co wymagało poświęcenia kilku dni na dość ciężką i bezowocną pracę.

Po trzecie zaangażowany byłem w rozmaite bieżące sprawy. Jedną z nich jest grodzisko.

Jakie to jest grodzisko – na razie pozostanie tajemnicą. Dość powiedzieć, że grodzisko to było badane przez archeologów w latach 70. i 80. XX wieku. Badania te  nigdy nie zostały w pełni opublikowane. Grodzisko ulega stałemu niszczeniu w wyniku prac polowych i obecnie jest nieomal nieczytelne w terenie.

Zacząłem od stworzenia GIS/AIS (Systemu Informacji Geograficznej/Archeologicznej) dla stanowiska oraz Numerycznego Modelu Terenu. Wracam więc do spraw, które były już poruszane wcześniej dla innego stanowiska. Tamten wpis oparłem na tutorialu Piotra Pociaska, który pragnę przypomnieć, choć nowsze wersje SAGA GIS (2.0.8) spowodowały, że pewne szczegóły interfejsu prezentowane wówczas się zmieniły.

Pierwszym krokiem była wektoryzacja zachowanego w archiwach planu warstwicowego w Qgis.

W tym celu wielką płachtę mapy warstwicowej, wykreślonej wówczas jeszcze za pomocą ołówka, kalki technicznej, linijki i kątomierza należało zeskanować i poddać georeferencji (na bazie usługi WMS z geoportalu).

Dalej należało ręcznie zwektoryzować poziomice tworząc tabelę atrybutów zawierającą wartość wysokości bezwzględnej.

Kolejnym było stworzenie Numerycznego Modelu Terenu w SAGA GIS:

Tutaj wersja stworzona metodą interpolacji (Interpolation from Points – Triangulation).

Tutaj wersja z nałożoną ortofotomapą.

Tutaj wersja stworzona metodą krigingu (kriging – ordinary kriging).

A tutaj otwarte okna podglądu 3D na trzy modele wykonane różnymi metodami.

SAGA i Qgis – tworzenie numerycznego modelu terenu dla stanowisk archiwalnych

Zdarza się niekiedy, że archeolog chciałby stworzyć numeryczny model terenu dla stanowiska archiwalnego. Mamy plan warstwicowy, zrobiony lata temu, przy pomocy tradycyjnych metod pomiarowych, wykreślony elegancko, zgodnie z geodezyjną sztuką. Chcielibyśmy uzyskać to, co w przypadku współczesnych metod pomiarowych nie jest specjalnie skomplikowane – trójwymiarowy obraz, czyli właśnie numeryczny obraz terenu.

Wykorzystamy w tym celu Qgis i SAGA GIS.

Dzisiaj modne są takie wizualizacje. Jeśli więc chcemy być modni to:

1. Musimy zdigitalizować plan warstwicowy, który będzie podstawą naszych działań. Najlepiej udać się do specjalistycznego punktu, gdzie dysponują wielkoformatowym skanerem (chyba, że sami takim dysponujemy).

2. Wprowadzamy skan mapy do Qgis korzystając z opcji georeferencji. Proces ten został już opisany tutaj.

3. Musimy teraz wykonać wektorowy plik shp dla warstwic. W tym celu klikamy „Nowa warstwa wektorowa”:

wybieramy opcję „linia”; w rubrykę nazwa wpisujemy np. „wartosc” (lub „wysokosc”, lub co tam chcemy); z opcji „typ” wybieramy REAL i wprowadzamy pole (ikonka z żółtą gwiazdką) oraz zatwierdzamy. Pliki nazywamy jak chcemy (najlepiej np. „warstwice”).

4. Odrysowujemy pracowicie wszystkie warstwice:

Przy edycji każdej izolinii wpisujemy jej wartość (wysokość nad poziom morza, lub względną), w tabeli atrybutów.

Po zakończonej pracy zapisujemy plik. Teraz będziemy mogli już wyjść z programu Qgis i otworzyć SAGA.

5. Otwieramy SAGA i wprowadzamy do niego plik shp z warstwicami: „modules – GDAL/OGR – Import Vector Data”. Jak już plik znajdzie się w programie trzeba wprowadzić go do mapy klikając dwukrotnie na jego warstwę w zakładce „Data”. Pojawi się okno dialogowe:

w której zaznaczacie opcję „Map” i:

6. Istnieje kilka metod interpolacji, które możemy wykorzystać. Jedną z nich jest Ordinary Kriging opisany na blogu Piotra Pociaska. Metoda ta daje dobre rezultaty, jednak o ile korzystacie z wersji SAGA dostarczanej w postaci binarek przez dewelopera (a nie kompilowanej samodzielnie), to okaże się, że nie będziecie mieli modułu Kriging. Jest on dostępny w wersji na Windows i w źródłach do samodzielnej kompilacji. Jeśli nie chcecie kompilować to:

• Stwórzcie warstwę z punktami: shapes – Points – Points from Lines, tak jak opisuje to Piotr Pociask:

• Teraz wybierzemy metodę triangulacji do stworzenia modelu:

Grid – Gridding – Interpolation from Points – Triangulation

• W oknie dialogowym będziecie musieli wprowadzić wszystkie potrzebne wartości. Później pojawi się okno definiowania „grida”, w którym wpisać trzeba rozmiar piksela. Jeśli model dotyczy niewielkiego obszaru (np. stanowiska), można wpisać „1”:

• Teraz wystarczy dodać uzyskany obrazek do mapy i wybrać z paska ikon tę z napisem „3D”. Uwaga: efekt uzyskany metodą Krigingu zazwyczaj wygląda nieco lepiej, więc albo użyjcie wersji Windows, albo samodzielnie skompilujcie program.

6. O ile chcecie by model wyglądał nieco inaczej możecie poeksperymentować. Wprowadzając do SAGA np. warstwy rastrowe z mapą lub ortofotomapą można uzyskać model z teksturą (ważne, żeby w zakładce „Maps” warstwy, które mają być wyświetlane na modelu znalazły się wyżej):

7. Możliwe, że będziecie też chcieli obrazek nieco uprościć. W tym celu zwróćcie uwagę na panel właściwości, wybierając zakładkę „Colors” (1), a w otwartym okienku dialogowym zakładkę „Presets”:

W otwartym kolejnym okienku dialogowym wybierzcie:

na przykład black green; wtedy efekt końcowy będzie wyglądał tak:

Nie jest to aż takie trudne…