Saya ingin menetapkan posisi panjang / lat tertentu pada peta ke ketinggian dari file data SRTM3, tetapi tidak tahu bagaimana menemukan nilai spesifik. Jadi saya ingin beberapa contoh bagaimana saya dapat menemukan di ketinggian N50E14.hgt ke 50 ° 24'58.888 "N, 14 ° 55'11.377" E.

Format data

Saya akan menganggapnya sebagai latihan kecil tentang cara memprogram pembaca data. Lihatlah dokumentasi :

Data SRTM didistribusikan dalam dua tingkat: SRTM1 (untuk AS dan wilayah dan harta miliknya) dengan data sampel pada interval satu busur-detik dalam lintang dan bujur, dan SRTM3 (untuk dunia) disampel pada tiga busur-detik.

Data dibagi menjadi ubin lintang dan bujur satu per satu derajat dalam proyeksi "geografis", yang berarti presentasi raster dengan interval lintang dan bujur yang sama dalam tanpa proyeksi sama sekali, tetapi mudah untuk memanipulasi dan membuat mosaik.

Nama file merujuk pada lintang dan bujur sudut kiri bawah ubin - misalnya N37W105 memiliki sudut kiri bawah pada 37 derajat lintang utara dan 105 derajat bujur barat. Untuk lebih tepatnya, koordinat ini merujuk ke pusat geometris dari piksel kiri bawah, yang dalam hal data SRTM3 akan mencapai sekitar 90 meter.

File tinggi memiliki ekstensi .HGT dan ditandatangani bilangan bulat dua byte. Byte berada dalam urutan "big-endian" Motorola dengan byte pertama yang paling signifikan, langsung dapat dibaca oleh sistem seperti Sun SPARC, Silicon Graphics dan komputer Macintosh menggunakan prosesor Power PC. DEC Alpha, kebanyakan komputer PC dan Macintosh yang dibangun setelah 2006 menggunakan Intel ("little-endian") sehingga beberapa byte-swapping mungkin diperlukan. Ketinggian dalam meter direferensikan ke geoid WGS84 / EGM96. Void data diberi nilai -32768.

Bagaimana cara melanjutkan

Untuk posisi Anda, 50 ° 24'58.888 "N 14 ° 55'11.377" E, Anda sudah menemukan ubin yang benar, N50E14.hgt. Mari kita cari tahu piksel mana yang Anda minati. Lintang pertama, 50 ° 24'58.888 "N:

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

detik busur. Dibagi dengan tiga dan dibulatkan ke bilangan bulat terdekat memberikan baris grid 500. Perhitungan yang sama untuk hasil garis bujur dalam kolom grid 1104.

Dokumentasi quickstart tidak memiliki informasi tentang bagaimana baris dan kolom disusun dalam file, tetapi dalam dokumentasi lengkap dinyatakan bahwa

Data disimpan dalam urutan utama baris (semua data untuk baris 1, diikuti oleh semua data untuk baris 2, dll.)

Baris pertama dalam file sangat mungkin yang paling utara, yaitu jika kita tertarik pada baris 500 dari tepi bawah , kita sebenarnya harus melihat baris

1201 - 500 = 701

dari awal jika file . Sel kisi kami adalah angka

(1201 * 700) + 1104 = 841804

dari awal file (yaitu lewati 700 baris, dan yang 701 ambil sampel 1104). Dua byte per sampel berarti kita harus melewati 1683606 byte pertama dalam file dan kemudian membaca dua byte untuk mendapatkan sel grid kita. Data adalah big-endian, yang berarti Anda harus menukar dua byte pada platform Intel misalnya.

Program sampel

Program Python sederhana untuk mengambil data yang tepat akan terlihat seperti ini (lihat dokumen untuk penggunaan modul 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)

Perhatikan bahwa pengambilan data yang efisien harus terlihat sedikit lebih canggih (mis. Tidak membuka file untuk masing-masing dan setiap sampel).

Alternatif

Anda juga bisa menggunakan program yang bisa membaca file .hgt di luar kotak. Tapi itu membosankan.

GDAL dapat membaca / menulis format raster ini dengan driver SRTMHGT . Ini berarti Anda dapat melihat raster dengan QGIS, ArcGIS, atau menggunakan utilitas GDAL seperti gdallocationinfo untuk mendapatkan nilai dari suatu titik, misalnya:

Ubah DMS ke DD:

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

Kemudian dari shell, gunakan gdallocationinfo file.hgt -wgs84 long lat:

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

Ketinggian 216 m.

Jika Anda menggunakan QGIS, periksa apakah plugin python "Point Sampling Tool" diinstal. Anda akan menemukannya di -> Enhancements (Python) -> Analisis.

Pilih layer titik Anda dari posisi yang diperlukan, kemudian mulai PST, pilih hgt (atau file raster / poligon apa pun) dan pilih bentuk titik baru untuk output.

Itu saja :-)

  Chris

Jawaban Chris menunjukkan langsung untuk mengambil sampel poin dari sebuah layer di QGIS.

Namun, karena balasan Anda untuk komentar saya menjelaskan Anda sedang menulis program Anda sendiri untuk membaca nilai ketinggian dari .hgtfile, lihat lagi Quickstart PDF di dokumen SRTM. Ini menjelaskan bagaimana data ketinggian disimpan. Untuk meringkas:

• File SRTM3 berisi urutan nilai integer big-endian.
• Setiap nilai integer mewakili ketinggian "dalam meter yang direferensikan ke geoid WGS84 / EGM96", kecuali untuk nilai -32768, yang menunjukkan piksel tanpa data.
• Ada 1201 baris dari 1201 sampel, sehingga harus ada 1442401 nilai integer secara keseluruhan.

Anda mengatakan bahwa Anda dapat mengkonversi antara koordinat lon / lat dan piksel, sehingga mendapatkan ketinggian adalah masalah membaca nilai integer dari offset yang sesuai dalam file. Diberikan koordinat piksel xdan yrelatif ke sudut kiri atas tempat kejadian, itu pada dasarnya offset = (y * 1201) + x. Pixel 0,0adalah integer pertama dalam file dan pixel 1200,1200adalah integer terakhir dalam file.