Monday 13 May 2013

Membuat Histogram Dengan GMT

Pada kesempatan kali ini, saya akan berbagi cara membuat histogram dengan GMT. Di script berikut, jumlah gempa akan diplot terhadap kedalaman gempa. Berikut adalah penggalan data gempa yang saya dapat dari USGS dan saya simpan dengan nama "gempa.dat" :
 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    ... .......   
...................................................................


Selanjutnya, buatlah file warna untuk topografi seperti berikut :
  -11000 36      38      175     -5000   56      58      195
  -5000  56      58      195     -1000   70      72      214
  -1000  70      72      214     -500    81      102     217
  -500   81      102     217     -150    100     129     223
  -150   100     129     223     -70     131     161     230
  -70    131     161     230     -20     164     192     240
  -20    164     192     240     0       170     200     255
  0      0       97      71      50      16      122     47
  50     16      122     47      500     232     215     125
  500    232     215     125     1200    161     67      0      
  1200   161     67      0       1700    130     30      30
  1700   130     30      30      2800    110     110     110
  2800   110     110     110     4000    255     255     255
  4000   255     255     255     6000    255     255     255
Simpan dengan nama "elev.cpt". 
Kemudian buatlah file warna untuk gempa seperti 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".

Kemudian buatlah baris-baris perintah berikut :
psbasemap -JM12 -X4.5 -R100/115/-15/0 -Ba2g2WSne:."HISTOGRAM": --HEADER_FONT_SIZE=20 -P -K -Y15> map3.ps
grdimage indo.nc -R -JM -Celev.cpt -K -O -Na -P >> map3.ps
gawk "{print $7, $6, $8}" gempa.dat | psxy -JM -R -Sc0.2c -G245/245/0 -W1 -P -O -K -Ctabel.cpt >> map3.ps
psxy -R -JM -W10.0 -Sf0.8i/0.1ilt -Gblack -P -O -K trench.gmt -m>> map3.ps
psscale -Ctabel.cpt -D-1/0.3/4/0.35 -B200:km: -O -X14.2c -Y1c -K>> map3.ps
gawk "{print $8}" gempa.dat | pshistogram -R0/800/0/1000 -JX9c -W10 -L2.0 -G0/0/255 -Ba100g50:"Kedalaman (km)":/a100g50:"Jumlah Gempa":WSne:."Grafik Kedalaman Vs Jumlah Gempa": --HEADER_FONT_SIZE=20 -P -O -X-13 -Y-13.5 >> map3.ps
ps2raster map3.ps -Gc:\programs\gs\gs8.53\bin\gswin32c -Tj -P -Ftectonic

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



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

1 comment:

  1. Halo kak, bisa dibuatkan videonya di youtube? Agar bisa mempermudah pemahaman, terimakasih ����

    ReplyDelete