Mapy hipsometryczne

Jak wykorzystać dane o numerycznym modelu terenu (NMT)? Skąd je wziąć i jak przygotować do wygodniej pracy?

Dzisiaj wpis mocno techniczny, bez szczególnych analiz. Gratka dla osób szukających sposobów na rozdzielenie danych punktowych na gminy, powiaty i województwa oraz wskazówki jak pokazywać takie dane na mapie (i agregować je zgodnie z podziałem terytorialnym). Oprzemy się na danych o wysokości terenu, ale to tylko przykład.

Potrzebne nam będą następujące biblioteki:

Kolejne fragmenty kodu to kilka skryptów złączonych w jedno. Stąd mogą trafić się powtórzenia w kodzie albo te same dane pod różnymi nazwami zmiennych. Nie przejmujcie się tym – czytając kod i rozumiejąc go znajdziecie sposób na wybranie tego, co najbardziej Was interesuje.

Zaczniemy od zgromadzenia danych o wysokości terenu. Centralny Ośrodek Dokumentacji Geodezyjnej i Kartograficznej (CODGiK) udostępnia takie dane w postaci plików tekstowych (jeden plik to jedno województwo). Struktura plików jest bardzo prosta – w każdym wierszu mamy informację dla jednego punktu: współrzędne punktu i jego wysokość nad poziomem morza. Pliki na serwerze CODGiK są spakowane – pobierzmy je i rozpakujmy.

Aby nie bawić się w ręczne pobieranie każdego pliku zeskanujmy automatycznie stronę i wyciągnijmy linki do plików. Nie raz już scrappowaliśmy strony, więc pójdzie łatwo:

Po tej operacji mamy 16 plików tekstowych w folderze dane/. Weźmy na warsztat Małopolskę.

Pierwszym krokiem będzie wczytanie danych i przygotowanie ich w odpowiednim układzie współrzędnych (takim, jaki znamy z geografii). Niestety pliki zapisane są w układzie PUWG 1992, musimy je przemapować. Oczywiście istnieją odpowiednie do tego funkcje, które potrafią to zrobić na obiektach typu Shape:

Dane zaokrągliliśmy, aby szybciej się rysowały. Oczywiście nie musicie tego robić. Pierwsza nasza mapka to wysokość w wybranym województwie:

Pięknie to wygląda, prawdziwa mapa hipsometryczna. Możemy zmienić paletę z ciągłej na dyskretną – przyczyna to pokazanie niebieskim kolorem tego co jest poniżej poziomu morza.

Nałóżmy nasze kolorki na rzeczywistą mapę – pobraną chociażby z Google Maps. Łatwiej będzie rozpoznać poszczególne górki i doliny. Niech będzie Kraków:

Ale to ciągle tylko dane z obszaru Krakowa. Co jest w którym miejscu?

Teraz ładnie widać na przykład Wzgórze Wawelskie. Oraz przyczynę krakowskiego smogu – położenie miasta w dolinie.

Jak wygląda rozkład wysokości na pokazanym obszarze?

Sprawdźmy co się stanie jak podniesie się woda w Krakowie – na przykład niech poziom morza będzie teraz na wysokości 150 metrów:

Widzimy, że całe koryto Wisły zostaje zapełnione wodą. Wawel ostaje się na lądzie. Między innymi do tego celu mogą służyć mapy z wysokością terenu.

No dobrze – wybraliśmy Kraków na postawie jakiegoś otoczenia współrzędnych jego środka. Ale na początku napisałem, że dowiemy się jak dane punktowe przypisać do gmin (i co za tym idzie powiatów i województw). Poniższy kod to właśnie robi dla jednego (małopolskiego) województwa.

Algorytm działania jest następujący (i stosunkowo prosty, chociaż obliczenia są długotrwałe):

  • przygotuj punkty w takim samym układzie jak dane o granicach gmin
  • dla wszystkich punktów z danymi przygotuj macierz binarną mówiącą o tym czy dany punkt należy do obszarów kolejnych gmin
  • na podstawie tej macierzy przypisz do punktu kod TERYT gminy (i powiatu, i województwa – bo te są fragmentami kodu gminy)
  • aby nie liczyć wszystkiego na raz dla wszystkich danych (potrzeba dużo pamięci) dzielimy punkty z danymi na mniejsze paczki

Dla jednego województwa proces trwał u mnie kilkanaście minut. Ale dzięki temu możemy w łatwy sposób wybierać dane o wysokości dla konkretnej gminy. Zobaczmy gminę w której leżą Rysy. Ale na początek – całe województwo małopolskie, na które nałożymy granice gmin:

Teraz wybierzmy jedną gminę – po jej kodzie TERYT. Do tego wybierzmy granice powiatu w którym leży gmina. I narysujmy taką mapkę:

Teraz podobnie jak wcześniej z Krakowem – dla łatwiejszej orientacji w położeniu naszych kolorowych punktów – nałóżmy je na mapę z Google Maps:

Fajnie jest mieć dane podzielone po kodach TERYT, prawda? Mamy ogrom danych – zamiast trzymać je w plikach wrzućmy je do bazy. Idea jest podobna jak poprzednio dla województwa, z tą różnicą że tutaj musimy wykonać to samo dla wszystkich 16 plików. I zamiast zapisywać plik lokalnie – dodajemy kolejne wiersze do stosownej tabeli. W dużej części jest to powtórzenie powyższego kodu (ale dla ułatwienia daję całość).

Gotowe! Tylko kilkadziesiąt godzin i sprawa załatwiona. Ale za to jak uprości to dalszą pracę! Zobaczmy sami.

Najpierw parametry dostępu do bazy danych:

A teraz się z nią łączymy i zapytaniami w SQL wyciągamy średnią wysokość na poziomie gminy, powiatu i województwa:

Możemy już zamknąć połączenie z bazą danych. Jeśli spojrzycie w wielkość danych jakie zajmują w pamięci będziecie zachwyceni ;) Całą matematykę (grupowanie i liczenie średnich) wykonał serwer bazodanowy – my mamy gotowy wynik.

Zobaczmy jak wyglądają średnie wysokości dla każdego z województw. Uwaga – poniżej do wczytania mapy (plików shp) i połączenia jej z danymi wykorzystuję pakiet sf – o wiele bardziej przyjemny w użyciu niż cała reszta readOGR() i tak dalej:

To samo możemy zrobić na poziomie powiatów:

oraz gmin:

W samych mapach oczywiście nie ma zaskoczenia – na południu mamy góry, a reszta to mniej więcej równina na jakiejś średniej wysokości… Właśnie – jakiej? Skorzystajmy z danych uśrednionych na poziom gmin:

Z rozkładu widzimy, że 3/4 kraju jest na wysokości nie większej niż około 225 metrów. Gdyby poziom wód podniósł się o te 225 metrów byłby prawdziwy potop…

Na koniec możemy przygotować dwie funkcje, które na mapie z Google Maps narysują nam wysokość terenu – na podstawie kodu TERYT gminy lub powiatu (korzystając z danych zapisanych w bazie):

Zobaczmy kilka przykładów:

Warszawa:

Kraków:

to czerwone po lewej to miejsce gdzie znajduje się klasztor w Tyńcu (Bielańsko-Tyniecki Park Krajobrazowy)

powiat Nowy Dwór Gdański

czyli Żuławy Wiślane (ich fragment) – zwróćcie uwagę, że to co jest mocno zielone jest poniżej poziomu morza.

Bieszczady:

pięknie widać połoniny (dwa czerwone pasma nad Wetliną) oraz masyw z Tarnicą i Krzemieniem (ta większa czerwona plama)

okolice Karpacza:

tam gdzie najbardziej czerwono leży Wielki Szyszak

Z dzisiejszego postu nauczyliśmy się (mam nadzieję) przede wszystkim:

    • jak sprawdzić czy dany punkt leży w określonym obszarze – to ćwiczyliśmy w dopisywaniu kodu TERYT do punktów
    • jak połączyć dane punktowe z mapą
    • że pakiet sf jest bardzo wygodny jeśli chodzi o pliki SHP
    • jak dane zapisać do bazy danych oraz je z niej pobrać

Jeśli Ci się podobało – udostępniaj, ślij znajomym. I wracaj tutaj w przyszłości.

2 komentarze do wpisu „Mapy hipsometryczne

  1. Super, zawsze wiele się uczę z tych przykładów!
    Przy okazji, zainteresowałem się tematem przekształceń, odwzorowania: kod:
    CRS_puwg1992 <- crs("+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs ")
    # układ WGS84 (współrzędne "normalne")
    CRS_wgs84 <- crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs;")

    W pakiecie simple features (SF) można posługiwać się czterocyfrowym kodem (EPSG) odwzorowań już przy ładowaniu lub przy przekształcaniu danych, np:

    dane_nmt %
    sample_n(size = 1000, replace=F) %>%
    st_as_sf( coords = c(„x”, „y”), crs = 2180, agr = „constant”, precision = 0.1) %>%
    st_transform(crs = 4326) %>%
    as_Spatial(.)

  2. korekta, ucieło pierwszą linię:

    dane_nmt =
    read_delim(dane_plik, ” „, col_names = c(„x”, „y”, „z”), col_types = „ddd”) %>%
    sample_n(size = 1000, replace=F) %>%
    st_as_sf( coords = c(„x”, „y”), crs = 2180, agr = „constant”, precision = 0.1) %>%
    st_transform(crs = 4326) %>%
    as_Spatial(.)

Dodaj komentarz