10. Mapy i poziomice

Zadanie Jubileuszowe…

…bo to już 10 zajęcia.

Na wykładzie była mowa, o różnego rodzaju wykresach. W tym o obrazowaniu powierzchni i poziomicach. Szukałem jakiegoś fajnego zadania. Wyszło mi dosyć zagmatwane. W związku z tym zadanie jest nieobowiązkowe (spóźniony prezent mikołajkowy), ale zachęcam do pobawienia się. Uwzględnię to w ocenach 😄.

Obrazowanie powierzchni

Geoportal

Jest taki serwis geoportal.gov.pl zawierajacy różniaste informacje na temat naszego kraju. Na prawdę różniaste. Pod adresem mapy.geoportal.gov.pl jest mapa z której można pobierać różne dane.

1. Wybieramy obszar, który nas interesuje i dobieramy odpowiedni zoom (albo kółkiem myszy, albo suwaczkiem po lewej stronie).

   

2. Usługa pozwala na pobieranie danych w różnych formatach

   • WCS (Usługa WCS (Web Coverage Service) – jest usługą pobierania danych przestrzennych, zapisanych w modelu rastrowym jak np. ortofotomapa czy dane numerycznego modelu terenu. )
   • WFS (Usługa WFS (Web Feature Service) – jest usługą służącą do pobierania danych w postaci wektorowej, na podstawie kryteriów zdefiniowanych przez użytkownika. Formatem służącym do przekazywania danych jest Geography Markup Language (GML).)

   Nas będzie interesował WCS., a zwłaszcza numeryczny model terenu..

3. Geoportal oferuje dwa zestawy danych:

   • Numeryczny model terenu (NMT)
   • Numeryczny model pokrycia terenu (NMPT)

   Pierwszy z nich opisuje kształt terenu, a drugi stara się odwzorować wszystko to co „wystaje nad powierzchnią”.

GIS

GIS System informacji geograficznych i oprogramowanie, które służy do robienia różnych rzeczy związanych z tworzeniem map. Co nie jest tematem naszych zajęć.

Mapping Toolbox

Z drugiej strony MATLAB posiada Mapping Toolbox, który…

Nie uważam, że jest to narzędzie absolutnie niezbędne w CV inżyniera specjalności BMI (Biomechanika Inżynierska??), ale kto wie co się w życiu przyda…

Do zrobienia

Korzystając z mapy dostępnej na Geoportalu wybrać jakiś (niezbyt wielki obszar, chyba nie można pobrać danych z obszaru większego niż 2 km^2^)

Będziemy zainteresowani Numerycznym Modelu Terenu w formacie Arc/Info ASCII Grid.

A podstawowym zadaniem będzie zobrazowanie powierzchni i narysowanie poziomic dla wybranego obszaru.

Format Arc/Info ASCII Grid

jest bardzo prosty. Jego opis zawiera na przykład Wikipedia na stronie Esri grid.

W naszym przypadku interpolowane¹ pomiary pochodzą z siatki o rozmiarach ca 1 m $\times$ 1 m. Nagłówek przedstawia podstawowe informacje na temat liczby wierszy (nrows) oraz kolumn (ncols), współrzędne lewego dolnego rogu obszary (xllcorner oraz yllcorner) w używanym ukłądzie współrzędnych. oraz rozmiar komórki (dx, dy). Całą reszta to zapisane „ciurkiem” interpolowane wartości wysokości.

Przykłądowy nagłówek wygląda tak:

ncols        462
nrows        147
xllcorner    364010.629344662302
yllcorner    361718.207231069566
dx           1.001068528529
dy           1.000744024319
 118.90000152587890625 118.8899993896484375 …

Układ współrzędnych

Ukłąd współrzędnych nie jest specjalnie istotny z naszego punktu widenia. Geoportal korzysta ze współrzędnych mapy 1992 (EPSG 2180).

Z naszego punktu widzenia nie jest to istotne. Możemy założyć, że startujemy z punktu o współrzędnych $(0,0)$ i współrzędne zmieniają się ze skokiem $dx,dy$.

Czytanie danych

Format ten rozumie funkcja MATLABa readgeoraster() z Mapping Toolbox. Natomisat funkcja mapshow() z tego samego pakietu pozwala wyrysowanie ładnej mapy. Nic prostszego.

>> [Z,R] = readgeoraster('result.asc',"OutputType","double");
>> mapshow(Z,R,'DisplayType','surface')

W tablicy Z zapisane zostaną przeczytane z pliku wysokości natomiast struktura R zawiera wszystkie wartości niezbędne do wyrysowania mapy.

Skąd wziąć Mapping Toolbox

Nie wiem, czy w laboratorium jest on zainstalowany, można to zprawdzić wydając polecenie MATLABa ver. Na moim laptopie:

>> ver
-----------------------------------------------------------------------------------------------------------
MATLAB Version: 24.2.0.2773142 (R2024b) Update 2
MATLAB License Number: 99999999
Operating System: Linux 6.8.0-49-generic #49-Ubuntu SMP PREEMPT_DYNAMIC Mon Nov  4 02:06:24 UTC 2024 x86_64
Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------------
MATLAB                                                Version 24.2        (R2024b)
Simulink                                              Version 24.2        (R2024b)
Curve Fitting Toolbox                                 Version 24.2        (R2024b)
Image Acquisition Toolbox                             Version 24.2        (R2024b)
Image Processing Toolbox                              Version 24.2        (R2024b)
MATLAB Compiler                                       Version 24.2        (R2024b)
Signal Processing Toolbox                             Version 24.2        (R2024b)
Statistics and Machine Learning Toolbox               Version 24.2        (R2024b)
Symbolic Math Toolbox                                 Version 24.2        (R2024b)

I widać, że nie jest zainstalowany. Nie warto go instalować na potrzeby wydania dwu komend.

1. Skorzystamy z MATLABa on-line (trzeba się zalogować)
2. Plik z danymi można przesłać na dysk w chmurze wchodząc na stronę drive.mathworks.com (może być potrzebne logowanie).

Odczyt danych z mapy

• wybieram pobieranie danych
• numeryczny model terenu
• wartstwa: DTM_PL-KRON86-NH (nieistotna)
• naciskam ołóweczek (Rysuj…) i zaznaczam prostokąt na mapie
• a następnie Pobierz plik

Dane zapisane są w formacie tekstowym nagłówek (6 linii) i kolejne linie z danymi, po jednej linii na każdy wiersz.

Plik przesyłam do chmury

Wydaję polecenia

>> [Z,R] = readgeoraster('dane.asc',"OutputType","double");
>> mapshow(Z,R,'DisplayType','surface')
>> demcmap(Z)

1. Pierwsza linia czyta dane
2. Druga rysuje ją
3. Trzecia dobiera kolory na podstawie wysokości w tablicy Z

Efekt

Żeby dostać poziomice

mapshow(Z,R,'DisplayType','contour')

lub lepiej

mapshow(Z,R,'DisplayType','contour','LineColor','k')

Można jeszcze dodać na poziomicach ich wysokości (dodajac parametr 'ShowText','on'), ale uczyni to obrazek nieczytelnym…

Sprawdzić https://www.mathworks.com/help/map/ref/contourm.html

O czym jeszcze warto pomyśleć?

1. Można zaimportować dane w formacie GeoTIFF. Nie sprawdzałem tego, ale funkcja readgeoraster() chyba też potrafi takie dane czytać i interpretować…
2. Numeryczny Model Pokrycia Terenu powinien (teoretycznie) pokazywać kształt budynków…

A co jeżeli nie Mapping Tool?

1. Plik otrzymany z geoportalu można zaimportować używając narzędzi Import Tool
2. Trzeba pamiętać, żeby w odpowiednim miejscu zaznaczyć, że poszczególne kolumny są oddzielane (Delimiter) odstępami (a nie przecinkami), i poprosić, aby wielokrotne odstępy traktował jako jeden (w Delimiter options), a później zaimportować tylko dane dotyczące wysokości (nie uwzględniając nagłówka)
3. Można również podejrzeć wartości dx i dy.
4. Do rysowania poziomic można użyć funkcji contour() a do narysowania mapy funkcji surf()

Funkcja surf()

Klasyczne użycie funkcji surf wymaga przygotowania siatki. Służy do tego funkcja meshgrid() generująca pary punktów układające się tak jak na obrazku

W tych punktach będą nam potrzebne wartości funkcji $z=f(x,y)$.

Klasyczne użycie funkcji surf() wygląda tak

[X,Y] = meshgrid(1:0.5:10,1:20); % generujemy siatkę
Z = sin(X) + cos(Y);             % wyliczamy wartości funkcji
surf(X,Y,Z)                      % rysujemy wykres

Ale można też tak:

x = 1:3;
y = 1:5;
[X,Y] = meshgrid(x,y)
X = 5×3
     1     2     3
     1     2     3
     1     2     3
     1     2     3
     1     2     3
Y = 5×3
     1     1     1
     2     2     2
     3     3     3
     4     4     4
     5     5     5

Wyliczamy $x^2 +y^2$ w punktach siatki

X.^2 + Y.^2
ans = 5×3
     2     5    10
     5     8    13
    10    13    18
    17    20    25
    26    29    34

Odczytana z geoportalu tablica zawiera właśnie funkcji (wysokość nad poziom morze) w węzłach siatki. Trzeba przygotować tylko wartośći X i Y. Potrzebne nam będą do tego wartości dx i dy oraz ncols i nrows odczytane z pliku ASCI Grid. Startować możemy albo z punktu o współrzędnych (xllcorner, yllcorner) albo od zera (nie budujemy mapy!)

W naszym przypadku x zmienia się w zakresie od 0 do (ncols - 1) * dx. Analogicznie y. Funkcja meshgrid() zrobi resztę. Spoglądając na wykres „z góry” dostaniemy pokolorowaną mapę… (Warto też zerknąć do slajdów z wykładu)

Funkcja contour()

Użycie analogiczne jak funkcji surf()

Ta sama siatka (X, Y) oraz wartości odczytane z mapy pozwolą na stworzenie poziomic.