Menggunakan SRTM Global DEM untuk perhitungan Slope?

17

Saya telah mengunduh SRTM GDEM (resolusi ~ 90 km).

Saya menggunakan ArcGIS 10.

Saya telah mencoba menggunakan analis spasial untuk menghitung kemiringan.

Namun, saya tidak dapat menghitung kemiringan.

Nilai output hanya memiliki dua rentang 0 dan 0,1-90.

Saya tidak begitu yakin apa masalahnya?

arcgis-desktop arcgis-10.0 dem srtm pengguna2543
sumber
Ini tergantung di mana Anda menganalisis di dunia. Ada berbagai proyeksi untuk setiap lokasi. Di mana Anda memeriksa?
DJJ
5
Resolusi sebenarnya ~ 90m, bukan ~ 90km.
Akheloes
Hanya komentar, jika Anda sedang dalam pemeliharaan untuk Desktop, Anda dapat masuk ke ArcGIS Online dan menggunakan layanan elevasi mereka (tanpa perlu ekstensi NA). Lapisan kemiringan bebas untuk digunakan sebagai lapisan referensi. Di Australia, kami memiliki data SRTM 1 detik (~ 30 juta) blogs.esri.com/esri/arcgis/2014/07/11/…
Simon

Jawaban:

29

Ini sepertinya tempat yang bagus untuk menggambarkan cara yang sederhana, cepat, dan lebih akurat untuk menghitung lereng untuk DEM yang luas secara global .

Prinsip

Ingatlah bahwa kemiringan permukaan pada suatu titik pada dasarnya adalah rasio terbesar "naik" ke "lari" yang dihadapi pada semua kemungkinan bantalan dari titik itu. Masalahnya adalah bahwa ketika proyeksi memiliki skala distorsi, nilai "run" akan dihitung secara salah. Lebih buruk lagi, ketika distorsi skala bervariasi dengan bantalan - yang merupakan kasus dengan semua proyeksi yang tidak sesuai - bagaimana kemiringan bervariasi dengan bantalan akan diperkirakan secara keliru, mencegah identifikasi akurat dari kenaikan maksimum: rasio run (dan memiringkan perhitungan aspek).

Kita dapat menyelesaikan ini dengan menggunakan proyeksi konformal untuk memastikan bahwa distorsi skala tidak bervariasi dengan bantalan, dan kemudian mengoreksi perkiraan kemiringan untuk memperhitungkan distorsi skala (yang bervariasi dari titik ke titik di seluruh peta). Caranya adalah dengan menggunakan proyeksi konformal global yang memungkinkan ekspresi sederhana untuk distorsi skalanya.

Proyeksi Mercator sesuai dengan tagihan: dengan asumsi skala tepat di Ekuator, distorsinya sama dengan garis potong garis lintang. Artinya, jarak pada peta tampaknya dikalikan dengan garis potong. Ini menyebabkan setiap perhitungan kemiringan untuk menghitung kenaikan: (dtk (f) * lari) (yang merupakan rasio), di mana f adalah garis lintang. Untuk memperbaikinya, kita perlu mengalikan lereng yang dihitung dengan detik (f); atau, secara setara, membaginya dengan cos (f). Ini memberi kami resep sederhana:

Hitung kemiringan (seperti naik: lari atau persen) menggunakan proyeksi Mercator, lalu bagi hasilnya dengan kosinus garis lintang.

Alur kerja

Untuk melakukan ini dengan kisi yang diberikan dalam derajat desimal (seperti SRTM DEM), lakukan langkah-langkah berikut:

  1. Buat garis lintang. (Ini hanya kotak koordinat y.)

  2. Hitung kosinusnya.

  3. Proyek kedua DEM dan cosinus dari garis lintang menggunakan proyeksi Mercator di mana skala benar di Khatulistiwa.

  4. Jika perlu, konversi satuan ketinggian agar sesuai dengan satuan koordinat yang diproyeksikan (biasanya meter).

  5. Hitung kemiringan proyeksi DEM baik sebagai kemiringan murni atau persen ( bukan sebagai sudut).

  6. Bagilah kemiringan ini dengan kisi proyeksi kosinus (garis lintang).

  7. Jika diinginkan, proyeksi ulang grid slope ke sistem koordinat lain untuk analisis atau pemetaan lebih lanjut.

Kesalahan dalam perhitungan kemiringan akan mencapai 0,3% (karena prosedur ini menggunakan model bumi bulat dan bukan model ellipsoidal, yang diratakan sebesar 0,3%). Kesalahan itu secara substansial lebih kecil dari kesalahan lain yang masuk ke perhitungan kemiringan sehingga dapat diabaikan.

Perhitungan sepenuhnya global

Proyeksi Mercator tidak dapat menangani kedua kutub. Untuk pekerjaan di daerah kutub, pertimbangkan untuk menggunakan proyeksi Stereografis kutub dengan skala sebenarnya di kutub. Distorsi skala sama dengan 2 / (1 + sin (f)). Gunakan ungkapan ini sebagai ganti detik (f) dalam alur kerja. Secara khusus, alih-alih menghitung kisi-kisi cosinus (garis lintang), hitung kisi-kisi yang nilainya (1 + sin (lintang)) / 2 ( sunting : gunakan -ukuran untuk Kutub Selatan, seperti yang dibahas dalam komentar). Kemudian lanjutkan persis seperti sebelumnya.

Untuk solusi global yang lengkap, pertimbangkan untuk memecah grid terestrial menjadi tiga bagian - satu di sekitar setiap kutub dan satu di sekitar khatulistiwa -, melakukan perhitungan kemiringan secara terpisah di setiap bagian menggunakan proyeksi yang sesuai, dan meredam hasilnya. Tempat yang masuk akal untuk membagi dunia adalah di sepanjang lingkaran lintang di lintang 2 * ArcTan (1/3), yaitu sekitar 37 derajat, karena pada lintang ini faktor koreksi Mercator dan Stereografis sama satu sama lain (memiliki nilai umum 5/4) dan akan menyenangkan untuk meminimalkan ukuran koreksi yang dilakukan. Sebagai pemeriksaan perhitungan, kisi-kisi harus dalam persetujuan yang sangat dekat di mana mereka tumpang tindih (sejumlah kecil ketidaktepatan titik apung dan perbedaan karena resampling dari kisi-kisi yang diproyeksikan harus menjadi satu-satunya sumber perbedaan).

Referensi

John P. Snyder, Proyeksi Peta - Manual Kerja . Kertas Profesional USGS 1395, 1987.

whuber
sumber
7
Saya menemukan diri saya pada posisi itu, seperti yang sering saya lakukan, mengucapkan terima kasih kepada Whuber sekali lagi karena menggambarkan solusi serta mengemukakan alasan yang membangunnya. Topiku tidak ada padamu, Pak.
matt wilkie
Terima kasih @matt. Saya tidak bermaksud menyarankan sebelumnya bahwa jawaban Anda (yang sekarang dihapus) harus dibatalkan: Sebenarnya, saya telah membatalkannya karena Anda membagikan tautan ke referensi USGS yang menarik yang dapat berguna bagi banyak pembaca. (Komentar saya adalah penting hanya dari bagian sekunder di kertas itu, bukan kertas itu sendiri.)
whuber
ahh terimakasih atas klarifikasinya. Saya telah memulihkan jawabannya, memercayai orang-orang memiliki cukup informasi di depan mereka sekarang untuk membuat pilihan yang terinformasi :)
matt wilkie
2
Berasal dari latar belakang bahasa Prancis, saya perlu beberapa saat untuk menerjemahkan terminologi yang diperlukan untuk lebih memahami jawaban yang bagus ini, jadi saya pikir menjatuhkan tautan ini adalah bantuan yang baik untuk pemula seperti saya: webhelp.esri.com/arcgisdesktop/9.2/…
Akheloes
Ini pendekatan hebat dan saya sudah menggunakan solusi Anda untuk menghasilkan raster lereng global. Satu petunjuk dari pengalaman praktis: Karena nilai garis lintang selatan khatulistiwa negatif, Anda harus menggunakan nilai garis lintang absolut dalam persamaan berikut: (1 + sin (garis lintang)) / 2
Saleika
18

Jawaban asli

Saya menduga unit horisontal untuk raster Anda berada dalam derajat atau detik busur. Anda perlu memproyeksikan ulang raster ini ke proyeksi spasial di mana unit horizontal dan vertikal Anda sama (yaitu, jika unit vertikal dalam meter, maka saya sarankan menggunakan UTM, yang memiliki unit horizontal meter).

Untuk memproyeksikan ulang raster dengan ArcCatalog / ArcGIS, lihat di:

ArcToolbox> Alat Manajemen Data> Proyeksi dan Transformasi> Raster> Project Raster

Pilih referensi spasial yang diproyeksikan yang mencakup wilayah yang Anda minati, misalnya, coba zona UTM. Ada banyak opsi lain yang paling baik didokumentasikan dalam manual . Catatan, Anda tidak dapat membuat dataset kemiringan untuk seluruh Earth (jika itu yang ingin Anda lakukan).

Jawaban yang lebih baik, menggunakan GDAL dengan skala

Sekarang data SRTM tersedia secara global, saya benar-benar dapat melihat dan bekerja dengan file. The gdaldemutilitas dari GDAL dapat menghitung kemiringan dan Hillshade menggunakan skala pilihan untuk rasio unit vertikal menjadi horisontal. Manual merekomendasikan 111120 m / ° untuk sesuatu seperti ubin SRTM. Jadi misalnya, dari shell OSGeo4W:

$ gdaldem slope -s 111120 -compute_edges N44E007.hgt N44E007_slope.tif

The -compute_edgespilihan membuat tepi lebih mulus, jika Anda ingin menjahit beberapa ubin bersama-sama. Atau menghitung ubin untuk wilayah yang luas. Kerugian dengan teknik "skala" adalah bahwa jarak dalam arah EW dan NS tidak sama, kecuali di khatulistiwa, jadi untuk ubin yang lebih dekat ke kutub, mungkin ada beberapa kesalahan representasi lereng yang aneh.

Mike T
sumber
Ada baiknya menekankan komentar terakhir Anda: ini adalah solusi yang buruk untuk poin tidak di dekat Khatulistiwa. Ini bukan masalah kecil "keliru penyajian": hasilnya akan sangat keliru, terutama di tempat-tempat yang lebih dekat ke Polandia daripada Khatulistiwa. Dokumentasi untuk gdaldemnegara bagian "Untuk lokasi yang tidak dekat dengan garis khatulistiwa, akan lebih baik untuk memproyeksikan ulang grid Anda menggunakan gdalwarp sebelum menggunakan gdaldem." Sayangnya itu tidak akan berfungsi untuk kumpulan data yang mencakup globe, kecuali jika Anda memecahnya menjadi potongan-potongan kecil (74 zona UTM, mungkin?), Memproyeksikannya, menghitung lereng, dan membuat mosaik hasilnya.
Whuber
7

Sederhananya, tidak ada satu. Menurut definisi, sistem koordinat berdasarkan derajat tidak diproyeksikan. Secara umum, kami mengatakan WGS84 adalah proyeksi "geografis", tetapi itu tidak benar, hanya untuk kenyamanan.

Saya rasa saya ingat pernah membaca tentang perangkat lunak atau proses untuk bekerja secara akurat dengan model elevasi di ruang geografis yang tidak diproyeksikan, tetapi saya tidak dapat menemukannya sekarang. Bagaimanapun itu akan menjadi percobaan atau membangun sendiri dari jenis proses kode.

Ahhh, menemukannya: Pengembangan Dataset Kemiringan Global untuk Estimasi Kejadian Tanah Longsor yang Dihasilkan dari Gempa Bumi (USGS). Halaman 4 menjelaskan masalahnya dengan baik

... panjang satu derajat bervariasi tergantung pada lokasi garis lintangnya. Di garis katulistiwa, blok satu derajat demi satu derajat cukup persegi ketika dikonversi ke satuan meter (111.321 meter di arah x dengan 110.567 meter di arah y ... tetapi lebih dekat ke kutub jarak di arah-x tumbuh lebih kecil sebagai fungsi kosinus lintang, karena konvergensi meridian. Sebagian besar paket GIS, termasuk ArcGIS, beroperasi hanya pada piksel persegi, dan menggunakan faktor untuk menyesuaikan dimensi x, y, atau z untuk unit umum tidak mungkin.

Makalah ini kemudian menjelaskan perhitungan spesifik dan alat perangkat lunak ( gdal , python , numpy ) yang mereka gunakan untuk menyelesaikan masalah mendasar ini. Makalah tidak menyertakan kode, tetapi jika diminta dengan baik mereka mungkin berbagi. Bagaimanapun, saya mungkin hanya bertanya di mana hasilnya, menjadi USGS itu mungkin sudah online di suatu tempat. :)

matt wilkie
sumber
1
Saran makalah itu bahwa proyeksi sama rata azimut dapat digunakan untuk menghitung lereng adalah salah arah dan salah. Memang akan memberikan lereng yang benar di dekat asal proyeksi, tetapi mereka juga akan semakin kurang akurat karena jarak ke asal meningkat.
whuber
terima kasih telah menunjukkan itu. Pembaca, pastikan untuk membaca gis.stackexchange.com/a/40464/108 juga, untuk keseimbangan
matt wilkie
2

Parameter DEM global (di mana sebagian besar rumus didasarkan pada asumsi ruang Euclidean) dapat diturunkan secara efisien menggunakan sistem EQUI7 GRID (Bauer-Marschallinger et al. 2014). EQUI7 GRID membagi dunia menjadi 7 wilayah daratan, semuanya diproyeksikan dalam sistem proyeksi yang sama dengan kehilangan presisi minimum. Lihat contoh DEM global pada resolusi 250 m di EQUI7 GRID. Di sini Anda dapat menemukan beberapa kode sampel yang menunjukkan cara memperoleh parameter DEM global menggunakan SAGA GIS. Setelah Anda menyelesaikan penurunan parameter DEM dalam sistem EQUI7 GRID, Anda dapat mengubah kembali semua peta menjadi longlatkoordinat WGS84 dan kemudian membuat mosaik global menggunakan GDAL.

Tom Hengl
sumber
Bisakah Anda menjelaskan bagaimana ini menjawab pertanyaan? Jika Anda mengusulkan menggunakan proyeksi yang sama untuk perhitungan kemiringan, harap perhatikan bahwa itu adalah solusi yang buruk karena distorsi metrik relatif besar yang terjadi ketika seseorang bergerak keluar dari pusat proyeksi. Meskipun memfokuskan tujuh proyeksi seperti itu pada massa tanah membantu mengatasi masalah itu, itu masih bukan pilihan terbaik.
whuber
Makalah oleh Bauer-Marschallinger et al. (2014) menjelaskan mengapa proyeksi ini dipilih untuk mewakili massa lahan global (mereka diasumsikan memiliki kehilangan minimum dalam presisi). Saya setuju bahwa setiap proyeksi 2D pada akhirnya akan mengarah pada deformasi, tetapi sejauh yang saya tahu, EQUI7 adalah kompromi yang baik antara kehilangan dalam presisi dan komoditas (aljabar 2D). Karena itu, heksagon lagi digunakan untuk mewakili permukaan tanah global (meskipun analisis DEM dengan heksagon 3D masih rumit).
Tom Hengl
Terima kasih untuk referensi Abstraknya menunjukkan bahwa ini memecahkan masalah yang sangat berbeda, yaitu "meminimalkan oversampling data lokal yang muncul selama proyeksi gambar satelit umum ke grid raster biasa." Itu tidak berarti itu akan berfungsi dengan baik untuk tujuan lain, seperti memperkirakan lereng.
whuber
Tentu saja EQUI7 tidak benar-benar menyelesaikan masalah estimasi lereng lokal secara akurat, tetapi ini mungkin solusi yang lebih elegan daripada menggunakan proyeksi Mercator yang disarankan di atas. Pada akhirnya, jika seseorang ingin memperkirakan kemiringan dengan ketepatan sempurna, maka satu-satunya pilihan mungkin adalah (1) menggunakan proyeksi lokal (berjarak sama) per ubin ukuran lebih kecil (misalnya 100 kali 100 km) dengan tumpang tindih 10-20%, seperti disebutkan juga dalam Verdin et al. (2007) , atau (2) untuk menggunakan kisi segi enam ( paket dggridR ).
Tom Hengl
Masalahnya bukan presisi - itu terletak pada pembuatan lereng dan aspek yang bias secara sistematis. Karena proyeksi berjarak sama secara berbeda mendistorsi arah ortogonal ke geodesik yang berasal dari pusatnya, aspeknya akan selalu salah (walaupun cukup akurat dekat dengan pusat di mana semua distorsi rendah) dan kesalahan pada lereng akan tumbuh dengan cepat dengan jarak. Penggunaan banyak proyeksi lokal tentu saja akan berhasil, tetapi sangat bertolak belakang dengan keanggunan yang Anda hargai.
whuber
-2

Kemiringan naik / lari. Hitung kenaikan dan hitung lari dan Anda memiliki jawaban Anda. Sederhana untuk menghitung jarak antara koordinat geografis. Ini akan menyebabkan lebih sedikit kesalahan pemasangan kembali dibandingkan dengan konversi ke UTM, dll.

THK
sumber