Chcę przypisać określoną pozycję długą / krótką na mapie do rzędnej z plików danych SRTM3, ale nie mam pojęcia, jak znaleźć konkretną wartość. Więc chcę jakiś przykład, jak mogę znaleźć w podniesieniu N50E14.hgt do 50 ° 24'58.888 "N, 14 ° 55'11.377" E.

Format danych

Zajmę się tym jako małe ćwiczenie w programowaniu czytnika danych. Zajrzyj do dokumentacji :

Dane SRTM są podzielone na dwa poziomy: SRTM1 (dla Stanów Zjednoczonych i ich terytoriów oraz posiadłości) z danymi próbkowanymi w odstępach jednosekundowych długości i szerokości geograficznej oraz SRTM3 (dla świata) próbkowanymi z trzema sekundami łukowymi.

Dane są dzielone na kafelki szerokości i długości geograficznej jeden na jeden stopień w projekcji „geograficznej”, co oznacza prezentację rastrową o równych odstępach szerokości i długości geograficznej bez projekcji, ale łatwą do manipulowania i mozaikowania.

Nazwy plików odnoszą się do szerokości i długości geograficznej lewego dolnego rogu płytki - np. N37W105 ma lewy dolny róg na 37 stopniach szerokości geograficznej północnej i 105 stopniach długości geograficznej zachodniej. Mówiąc ściślej, współrzędne te odnoszą się do geometrycznego środka lewego dolnego piksela, który w przypadku danych SRTM3 będzie wynosił około 90 metrów.

Pliki wysokości mają rozszerzenie .HGT i są podpisane dwubajtowymi liczbami całkowitymi. Bajty są w kolejności „big-endian” firmy Motorola, z najważniejszym bajtem jako pierwszym, bezpośrednio odczytywanym przez systemy takie jak Sun SPARC, Silicon Graphics i komputery Macintosh wykorzystujące procesory Power PC. DEC Alpha, większość komputerów PC i Macintosh zbudowanych po 2006 roku korzysta z kolejności Intel („little-endian”), więc może być konieczna zamiana bajtów. Wysokości podano w metrach w odniesieniu do geoidy WGS84 / EGM96. Pustki danych mają przypisaną wartość -32768.

Jak postępować

Dla swojej pozycji, 50 ° 24'58.888 "N 14 ° 55'11.377" E, już znalazłeś prawidłową płytkę, N50E14.hgt. Sprawdźmy, który piksel Cię interesuje. Pierwsza szerokość geograficzna, 50 ° 24'58.888 "N:

24'58.888" = (24 * 60)" + 58.888" = 1498.888"

sekundy łukowe. Podzielony przez trzy i zaokrąglony do najbliższej liczby całkowitej daje wiersz siatki 500. To samo obliczenie długości geograficznej daje kolumnę 1104 siatki.

W dokumentacji szybkiego startu brakuje informacji o tym, jak wiersze i kolumny są zorganizowane w pliku, ale w pełnej dokumentacji jest to zaznaczone

Dane są przechowywane w głównej kolejności wiersza (wszystkie dane dla wiersza 1, a następnie wszystkie dane dla wiersza 2 itd.)

Pierwszy wiersz w pliku jest najprawdopodobniej tym najbardziej wysuniętym na północ, tzn. Jeśli interesuje nas wiersz 500 od dolnej krawędzi , musimy spojrzeć na wiersz

1201 - 500 = 701

od początku, jeśli plik . Nasza komórka siatki to liczba

(1201 * 700) + 1104 = 841804

od początku pliku (tj. pomiń 700 wierszy, a w 701-szym weź próbkę 1104). Dwa bajty na próbkę oznaczają, że musimy pominąć pierwsze 1683606 bajtów w pliku, a następnie odczytać dwa bajty, aby uzyskać komórkę siatki. Dane są typu big-endian, co oznacza, że ​​musisz zamienić dwa bajty np. Na platformach Intela.

Przykładowy program

Uproszczony program Pythona do pobierania odpowiednich danych wyglądałby tak (patrz dokumentacja dotycząca korzystania z modułu struct):

import struct

def get_sample(filename, n, e):
    i = 1201 - int(round(n / 3, 0))
    j = int(round(e / 3, 0))
    with open(filename, "rb") as f:
        f.seek(((i - 1) * 1201 + (j - 1)) * 2)  # go to the right spot,
        buf = f.read(2)  # read two bytes and convert them:
        val = struct.unpack('>h', buf)  # ">h" is a signed two byte integer
        if not val == -32768:  # the not-a-valid-sample value
            return val
        else:
            return None

if __name__ == "__main__":
    n = 24 * 60 + 58.888
    e = 55 * 60 + 11.377
    tile = "N50E14.hgt"  # Or some magic to figure it out from position
    print get_sample(tile, n, e)

Zauważ, że wydajne pobieranie danych musiałoby wyglądać na bardziej wyrafinowane (np. Nie otwieranie pliku dla każdej próbki).

Alternatywy

Możesz także użyć programu, który może odczytywać pliki .hgt po wyjęciu z pudełka. Ale to jest nudne.

GDAL może odczytywać / zapisywać te formaty rastrowe za pomocą sterownika SRTMHGT . Oznacza to, że możesz wyświetlić raster za pomocą QGIS, ArcGIS lub użyć narzędzi GDAL, takich jak gdallocationinfo, aby uzyskać wartości z punktu, np .:

Konwertuj DMS na DD:

• Lat: 50 ° 24'58.888 "N = 50 + (24/60) + (58.888 / 3600) = 50.4163577778
• Długi: 14 ° 55'11.377 "E = 14 + (55/60) + (11.377 / 3600) = 14.9198269444

Następnie z powłoki użyj gdallocationinfo file.hgt -wgs84 long lat:

$ gdallocationinfo N50E14.hgt -wgs84 14.9198269444 50.4163577778
Report:
  Location: (1104P,700L)
  Band 1:
    Value: 216

Wysokość wynosi 216 m.

Jeśli korzystasz z QGIS, sprawdź, czy wtyczka Pythona „Point Sampling Tool” jest zainstalowana. Znajdziesz go w -> Ulepszenia (Python) -> Analizuj.

Wybierz warstwę punktową wymaganych pozycji, a następnie rozpocznij PST, wybierz hgt (lub dowolny plik rastrowy / wielokątny) i wybierz nowy kształt punktowy do wydruku.

To wszystko :-)

  Chris

Odpowiedź Chrisa wskazuje, że próbkowanie punktów z warstwy w QGIS jest proste.

Ponieważ jednak odpowiedź na mój komentarz wyjaśnia, że ​​piszesz własny program do odczytu wartości wysokości z .hgtplików, zapoznaj się z dokumentem Quickstart PDF w dokumentach SRTM. Wyjaśnia, w jaki sposób przechowywane są dane wysokości. Podsumowując:

• Pliki SRTM3 zawierają sekwencję wartości całkowitych big-endian.
• Każda wartość całkowita reprezentuje wzniesienie „w metrach odniesione do geoidy WGS84 / EGM96”, z wyjątkiem wartości -32768, które wskazują piksele bez danych.
• Jest 1201 linii 1201 próbek, więc powinny być w sumie 1442401 wartości całkowite.

Mówisz, że możesz konwertować między współrzędnymi lon / lat i pikselami, więc uzyskanie wysokości jest kwestią odczytu wartości całkowitej z odpowiedniego przesunięcia w pliku. Biorąc pod uwagę współrzędne pikseli xi ywzględem lewego górnego rogu sceny, to w zasadzie offset = (y * 1201) + x. Piksel 0,0to pierwsza liczba całkowita w pliku, a piksel 1200,1200to ostatnia liczba całkowita w pliku.