Tuesday, 14 May 2013

Membuat Cross Section Gempa Dengan GMT

Cross section atau irisan melintang kadang diperlukan untuk melihat sebaran gempa di suatu wilayah dan untuk melihat kemiringan penunjaman lempeng. Kali ini saya akan berbagi cara membuat cross section tersebut. Sebelumnya buatlah file text berisi data gempa. Berikut adalah contoh data gempa yang saya miliki :
 PDE    1973  01 02 005320.30  -9.85  117.43  66  5.5 mbGS    ... .......     
 PDE    1973  01 02 022709.20   1.03  126.21  61  5.4 mbGS    ... .......     
 PDE    1973  01 07 221702.60   5.68  127.30  79  5.1 mbGS    ... .......     
 PDE    1973  01 09 061425.30   6.98  126.15  51  5.3 mbGS    ... .......     
 PDE    1973  01 13 111020.40  -2.70  101.27 105  5.2 mbGS    ... .......     
 PDE    1973  01 16 054257.90   0.54  125.96  15  5.3 mbGS    ... .......     
..............................................................................
Simpan dengan nama "gempa.dat".
Kemudian berikut adalah gempa utama yang terjadi :
PDE    1973  01 02 005320.30  -8.308  109    92  5.5 mbGS    ... .......  
Simpan dengan nama "utama.dat".

Kemudian  buatlah file warna berikut :
0    255 0 0        100    255 0 0
100    0 255 0        300    0 255 0
300    0 0 255        800    0 0 255
simpan dengan nama "tabel.cpt". 
Bila pada postingan sebelumnya saya membuat manual file warna untuk topografi, kali ini saya akan memasukkan script untuk membuat file warna secara otomatis.

 Selanjutnya buatlah script berikut :

makecpt -Cglobe -Z > elev.cpt
psbasemap -JM12 -X4.5 -R100/115/-15/0 -Ba2g2WSne:."CROSS SECTION": --HEADER_FONT_SIZE=20 -P -K -Y15 > map2.ps
pscoast -JM -R -Ggreen -Sblue -Dh -Wthin -K -P -O >> map2.ps
grdimage indo.nc -R -JM -Celev.cpt -K -O -Na -P >> map2.ps
gawk "{print $7, $6, $8}" gempa.dat | psxy -JM -R -Sc0.2c -W1 -Gred -P -O -K -Ctabel.cpt>> map2.ps
gawk "{print $7, $6}" utama.dat | psxy -JM -R -Sa0.9c -W1 -Gyellow -P -O -K >> map2.ps
psxy -R -JM -W10.0 -Sf0.8i/0.1ilt -Gblack -O -K trench.gmt>> map2.ps
echo 106.0 -10.0 14 0 1 LT A> line.dat
echo 114.0 -5.2 14 0 1 LT B>> line.dat
psxy line.dat -JM -R -W10 -P -O -K >> map2.ps
pstext line.dat -JM -R -P -O -K >> map2.ps
gawk "{print $7, $6, $8}" gempa.dat | project -Q -C106.0/-10.0 -E114.0/-5.2 -Fpz -W-100/100 > cross.dat
psbasemap -JX15/-10 -X-1 -R0/700/0/350 -Ba100:"Distance (km)":/:"Depth (km)":50WSne --LABEL_FONT_SIZE=18p -P -O -K -Y-12.5 >> map2.ps
echo 0 0 14 0 1 BC A> teks.dat
echo 700 0 14 0 1 BC B>> teks.dat
pstext teks.dat -JX -R -P -O -Gblue -N -K >> map2.ps
psxy cross.dat -JX -R -Sc0.2 -W1 -Gred -P -O -K >> map2.ps
gawk "{print $7, $6, $8}" utama.dat | project -Q -C106.0/-10.0 -E114.0/-5.2 -Fpz -W-100/100 > maincross.dat
psxy maincross.dat -JX -R -Sa0.9 -W1 -Gyellow -P -O >> map2.ps
ps2raster map2.ps -Gc:\programs\gs\gs8.53\bin\gswin32c -Tj -P -Fmap2

Simpan dengan ekstensi ".bat". Setelah dieksekusi hasilnya akan seperti berikut :



Perintah makecpt pada script diatas adalah untuk membuat file warna topografi dan bathymetri. Untuk pilihan warna lainnya bisa dilihat di command prompt dan ketik makecpt. Sedangkan perintah project disini digunakan untuk memproyeksikan data pada sebuah lintasan.

Selamat mencoba dan semoga bermanfaat. Salam Orang Indonesia...

21 comments:

  1. Yos, itu cross sectionnya gak terbalik kah ? koq yg di wilayah A datar, tapi ke arah B malah menujam #jamroni

    ReplyDelete
    Replies
    1. bukannya semakin ke B emg gempanya semakin dalam ya pak? keliatan datar krn gempanya blm direlokasi pak...

      Delete
  2. kq saya coba script di atas gak keluar cross section nya ya??

    ReplyDelete
    Replies
    1. format data gempanya sama kayak di sini bang?

      Delete
    2. klo peta seismisitasnya muncul (gambar atas), tapi cross sectionnya gak muncul (gambar bawah),, kenapa ya??

      Delete
    3. kedalamannya harus bernilai positif bang, trus di psbasemap bawah yang -R0/700/0/350 harus disesuaikan ama datanya, 0-700 itu jaraknya, yg 0-350 itu kedalamannya

      Delete
    4. iya datanya ud bner kq kedalamannya udah positif,, trus yg dibagian itunya uda gw ganti jarak sama kedalamannya,, tapi ttp gak muncul,, trus kan contoh yg lw kasih ad gwmpa utamanya,, yg gw bikin gak pake gempa utama,,ap itu ngaruh ya??

      Delete
    5. coba di dua baris terakhir sebelum ps2raster dihapus bang, yg gawk ama psxy itu

      Delete
  3. This comment has been removed by the author.

    ReplyDelete
  4. Halo mas Yosi, saya Thomas, saya ingin belajar GMT, mau tanya apa bedanya GMT dan iGMT, apa mas Yosi pernah pake iGMT, kayaknya iGMT itu GUI nya utk GMT ya? Tks

    ReplyDelete
  5. Halo juga pak Thomas, saya sendiri belum pernah pakai iGMT, mungkin nanti akan saya coba dulu

    ReplyDelete
  6. Mas, kalo mau buat gambar yg bagian bawah aja (tetep overlay) tapi sumbu x-nya bujur atau lintang gimana ?
    makasih..

    ReplyDelete
    Replies
    1. di bagian
      gawk "{print $7, $6, $8}" utama.dat | project -Q -C106.0/-10.0 -E114.0/-5.2 -Fpz -W-100/100 > maincross.dat
      psxy maincross.dat -JX -R -Sa0.9 -W1 -Gyellow -P -O >> map2.ps

      di bagian "project" yg -Fpz ganti saja dengan -Fxyzpqrs, trus ntar di psxy nya diganti jd
      gawk "{print $(no kolom lintang/bujur), $(no kolom kedalaman)}" maincross.dat | psxy -JX -R -Sa0.9 -W1 -Gyellow -P -O >> map2.ps
      no kolom bisa dilihat di file "maincross.dat"

      Delete
  7. mas yosi, tanya mas... misalkan kita pakai data dari global cmt, barisan yg

    gawk "{print $7, $6, $8}" gempa.dat | psxy -JM -R -Sc0.2c -W1 -Gred -P -O -K -Ctabel.cpt>> map2.ps
    gawk "{print $7, $6}" utama.dat | psxy -JM -R -Sa0.9c -W1 -Gyellow -P -O -K >> map2.ps

    gawk "{print $7, $6, $8}" gempa.dat | project -Q -C106.0/-10.0 -E114.0/-5.2 -Fpz -W-100/100 > cross.dat

    gawk "{print $7, $6, $8}" utama.dat | project -Q -C106.0/-10.0 -E114.0/-5.2 -Fpz -W-100/100 > maincross.dat

    perlu atau bisa dihilangkan mas??

    terimakasih mas,
    regards

    Yogaswara

    ReplyDelete
    Replies
    1. dihapus gpp, itu perintah untuk membuat cross section gempa, klo mau cross section mekanisme fokal pakai pscoupe aja

      Delete
  8. kakak mohon pencerahannya, aku mau coba bikin peta ini, pas semua code nya tak copas persis kok ada error di psbasemap ya?? gini "Warning: Axis sub-item a set more than once (typo?)

    ReplyDelete
  9. bisa discreenshot errornya, by email gpp

    ReplyDelete
  10. isian kolom 1,2,3 dst itu apa ya mas?

    ReplyDelete
  11. Mas Yosi, gmn caranya biar gempa utama yg berbentuk bintang itu muncul lebih dari 1 di peta cross section?
    saya udh coba membuat garis cross section yg melalui 2 gempa utama, tetapi hanya 1 yg muncul di peta cross section nya.
    Mohon pencerahannya ya..
    Terima kasih

    ReplyDelete
  12. This comment has been removed by the author.

    ReplyDelete
  13. maaf kak, data trench.gmt nya dimana di donlotnya ya?

    ReplyDelete