HEC-RAS 2D i QGIS

Pierwsza część mojego przyspieszonego kursu modelowania w HEC-RAS, którą znajdziecie pod tym linkiem, została poświęcona w całości modelowaniu 2D na poziomie (bardzo) podstawowym. W wyniku naszych analiz otrzymaliśmy dane mówiące nam o głębokościach zalewu, rzędnych zwierciadła wody oraz prędkościach przepływu. Są to 3 podstawowe i domyślnie tworzone warstwy wynikowe. HEC-RAS umożliwia utworzenie dodatkowych warstw, takich jak na przykład iloczyn prędkości i głębokości (DV), dający nam bardzo dobre rozeznanie o obszarach szczególnie niebezpiecznych dla ludzi, albo coś przydatne głównie z hydraulicznego punktu widzenia jak warstwy z wartościami liczb Frouda lub Couranta.

HEC-RAS na początek

W tym celu, musimy otworzyć Ras Mappera, rozwinąć zakładkę Results, kliknąć prawym przyciskiem myszy (PPM) na nazwę naszego planu i wybrać Create a New Results Map Layer. Rozwijamy interesujący nas typ mapy, wybieramy czy chcemy uzyskać wartości minimalne, maksymalne czy dla konkretnego kroku czasowego i potwierdzamy Add map. Warstwa zostanie dodana, ale musimy ją jeszcze przeliczyć klikając na Compute/Update Stored Map. W nawiasie może się również pojawić krótka informacja taka jak Map not created, Map files up to date lub inna. Warto również zmienić w Layer Properties legendę i dla DV ustawić górny zakres na 2, przy czym warto wiedzieć, że już dla wartości 0,5 istnieje poważne ryzyko dla osób dorosłych.

Kiedyś wydrukuję sobie taką mapę i zawieszę na ścianie, bo jest piękna.

Przechodzimy do QGISa

Jak już zostało wcześniej wielokrotnie wspomniane, Ras Mapper, delikatnie mówiąc, nie jest najlepszym narzędziem gisowym, w jakim dane mi było pracować, więc zapiszemy nasze wyniki do pliku tiff i dalej będziemy pracować w QGIS. Twórcy oprogramowania chyba to przewidzieli i bardzo ułatwili nam eksport map. W przypadku głębokości, rzędnych i prędkości wystarczy kliknąć ppm na warstwę, rozwinąć Export Layer i wybrać Export Raster, a pozostałe stworzone przez nas warstwy od razu zapisują się do odpowiedniego formatu. Wszystkie te dane znajdziemy w folderze z naszym projektem, a jeszcze dokładniej – w folderze z nazwą naszego planu.

Warto korzystać z otwartych źródeł danych

Po dodaniu tychże rastrów do QGISa warto zmienić symbologię, coś niebieskiego dla głębokości, a coś z czerwienią dla DV. My w tym przykładzie ograniczymy się tylko do tych dwóch warstw, przyjrzymy się bowiem poziomem zalania zabudowy, ale rzucimy również okiem na obszary szczególnie niebezpieczne dla przebywających w pobliżu osób. Zapytacie, moi drodzy, skąd wziąć warstwę z budynkami. Śpieszę z odpowiedzią! Źródła są co najmniej 2, OpenStreetMap oraz Baza Danych Obiektów Topograficznych (BDOT10k). W pierwszym przypadku dane możemy pobrać ze strony Geofabrik.de, a w drugim z Geoportalu. Tylko zwróćcie uwagę na układy współrzędnych!

Usuwamy to co zbędna

Nie potrzebujemy obiektów dla całego województwa czy powiatu, dlatego od razu po wczytaniu pliku shapefile do QGISa ograniczymy je wyłącznie do zakresu naszego modelu. Nie chcemy przecinać budynków, więc darujemy sobie narzędzie Clip i zamiast niego najpierw wybierzemy poprzez selekcję i dalej ekstrakcję te poligony, które przecinają się ze stworzony przez nas już wcześniej zasięgiem obszaru 2D. W Processing Toolbox wyszukujemy Extract by Location. Extract features from będzie naszą warstwą z budynkami, warunkiem, który musi zostać spełniony będzie Intersect, a By comparing to the feature from to nasz obszar 2D, po czym zapiszemy wyselekcjonowane poligony do nowego pliku.

Analiza wyników

Chcielibyśmy teraz sprawdzić dla każdego każdego budynku 4 rzeczy: jaka część jego powierzchni może być zalana (liczba komórek), jaka jest średnia, maksymalna oraz minimalna głębokość wody. Do tego celu wykorzystamy Zonal Statistics. Naszą Input layer jest właśnie stworzony nowy plik shp, Raster layer to raster z głębokościami wody i interesującymi nas Statistics to calculate są wspomniane Count, Mean, Minimum oraz Maximum. Po ich zaznaczeniu koniecznie potwierdzamy wybór przyciskiem OK, w innym przypadku nie zostanie to zapamiętane. Tym razem jednak darujemy sobie tworzenie nowego pliku i zapiszemy całość jako warstwę tymczasową.

Łączenie danych

Po wyłączeniu programu stracilibyśmy efekt naszej pracy, więc żeby tego uniknąć posłużymy się funkcją Joins, do której dostajemy się wchodzą w Properties warstwy z analizowanymi budynkami. Zielonym krzyżykiem dodajemy nowe połączenie. Join layer to warstwa tymczasowa Zonal Statistics, Join field oraz Target field zostawimy fid i dołączymy tylko 4 obliczone pola zaznaczając najpierw Joined fields, a następnie to co stworzyliśmy. Nowe kolumny charakteryzują się tym, że mają podkreślenie jako prefix. Potwierdzamy OK, Apply i znowu OK.

Modyfikujemy tabelę z wynikami

Nasza tabela atrybutów jest pełna niepotrzebnych informacji, usuniemy więc to co zbędne wchodząc w tryb edycji i wybierając Delete Fields. Wybieramy wszystko poza fid. Prawdopodobnie zauważyliście, że część jest wyszarzona – to te kolumny, które dołączamy. Będą one tam tak długo, do póki istniej warstwa, do której mogą się odnieść, a nasza jest tymczasowa. Jesteśmy więc zmuszeni jakoś sobie z tym poradzić, najlepiej więc będzie, jeżeli wykorzystamy kalkulator. Klikamy więc Open field calculator i dla każdej wartości obliczamy nowe pole na podstawie dołączenia. Zostawiamy domyślny wybór Create a new field, Output filed name to nazwa naszej nowej kolumny, Output field type będzie liczną zmiennoprzecinkową Decimal number (real), z dokładnością (Precision) do 2 miejsc. W środkowym okienku rozwijamy Fields and Values i wyszukujemy aktualnie interesującą nas rzecz.

Pojawia się jednak problem, w momencie, w którym mamy zalaną tylko jedną komórkę na obszarze budynku, możliwe jest, że kalkulator zaokrągli naszą wartość w dół do 0. Możemy temu zapobiec zmieniając nasze wyrażenie. Wracamy więc do kalkulatora, ale tym razem wybieramy Update existing field, którym jest w moim przypadku Count. W wyrażeniu nie dodajemy zwyczajnie naszej warstwy, lecz wpiszemy coś odrobinę bardziej skomplikowanego.

CASE
WHEN “Zonal Statistics__count” > 0 AND “Zonal Statistics__count” < 1 THEN 1
ELSE “Zonal Statistics__count”
END

Wizualizacja wyników

Teraz powinno już wszystko grać. Kolejną rzeczą, którą warto zrobić, jest zmiana symbologii, żebyśmy od razu widzieli, które budynki są zagrożone. Wchodzimy więc we właściwości warstwy z budynkami poprzez Properties, następnie przechodzimy do Symbology, po czym zmieniamy z Single Symbol na Rule-based. Sklasyfikujemy budynki na 3 kategorie, silnie zagrożone, gdzie głębokość wody będzie wyższa niż 0,3 m (“MAX” >= 0.3), mało zagrożone, dla których ten warunek będzie spełniony w przedziale pomiędzy 0,1 a 0,3 m (“MAX” > 0.1 AND “MAX” <0.3) oraz bezpieczne (ELSE), czyli głębokość wody to nie więcej niż 10 cm.

Warto dodać też opis

Analogicznie do symbologii formuły wykorzystanej przy aktualizacji pól, możemy również przygotować etykiety, czyli Labels. W Properties naszej warstwy wyszukujemy Labels, po czym wpisujemy poniższą formułę:

CASE
WHEN “COUNT”> 1 THEN’Zalana powierzchnia budynku: ‘|| “COUNT”||’ m. kw.’
END

Dzięki czemu wszystkie zalane budynki zostaną opatrzone stosowną informacją. Żeby jednak to wyświetlić, musimy kliknąć ppm na naszą warstwę i wybrać Show labels. Ostatnią rzeczą jaką zrobimy jest usunięcie połączenia do warstwy tymczasowej ze statystykami głębokości i powierzchni zalewu.

Zastanawiacie się pewnie co z warstwą DV? Ano nic. Przecież tam gdzie są budynki prędkość powinna wynosić 0, więc nie ma sensu pokazywać tam niczego poza samą mapą 🙂