Membuat Peta Seismisitas Berbasis Waktu

Kali ini saya kan membagikan cara membuat peta seismisitas. Namun peta seismisitas yang saya buat, pewarnaannya berdasarkan waktu selisih antara waktu referensi dengan origin time gempa. Sehingga akan terlihat pola sebaran gempa yang terjadi dalam waktu yang saling berdekatan. Disini saya menggunakan data gempa Palu yang saya dapatkan dari repogempa.bmkg.go.id. Anda dapat menggunakan data saya disini.

Selanjutnya buatlah skrip berikut :

set base=2018 09 27 00 00 00

gawk -F"/|:| " "NR>12 {t2=mktime(\"%base%\"); tmst=sprintf(\"%%s %%s %%s %%s %%s %%s\", substr($1,1,4), substr($2,1,2), substr($3,1,2), substr($5,1,2), substr($6,1,2), substr($7,1,6)); t1=mktime(tmst); print (t1-t2)/3600, $0}" datagempa.csv | gawk "{if ($5==\"S\") {$4=$4*-1}; print $1, $4, %$6, $8, $9}" > dataolah

Di skrip tersebut, saya menggunakan tanggal referensi 27 September 2018. Kemudian saya menghitung selisih waktu dalam satuan jam antara tanggal referensi dengan origin time gempa dan menyimpannya dalam file dataolah. NR>12 artinya saya menggunakan data setelah baris ke 12 karena terdapat header pada baris 1 hingga 12 yang ingin saya abaikan.

Selanjutnya, kita harus membuat file palette warna sesuai dengan data selisih waktu yang didapat sebelumnya.Tetapi kita harus tahu nilai dari data yang akan kita buat palette warnanya terlebih dahulu. Caranya, ketikkan perintah berikut :

gmtinfo dataolah

Maka akan muncul informasi seperti ini :

dataolah: N = 704 <0.405278/812.447> <-4.58/0.69>  <118.71/123.45> <10/255>   <2.2/7.5> 

Artinya, terdapat 704 data gempa. Dimana nilai didalam tanda kurung <..../....> menunjukkan nilai maksimum dan minimum di tiap kolom. Data selisih waktu terdapat di kolom 1. Sehingga kita menggunakan nilai 0.405278 sebagai nilai minimum selisih waktu dalam setuan jam dan 812.447 sebagai nilai maksimum selisih waktu dalam satuan jam.

Kemudian, untuk membuat file palette warna, ketikkan perintah berikut :

makecpt -Cseis -T0/813 > warna.cpt

Nilai 0 dan 813 kita dapatkan dengan membulatkan data minimum dan maksimum dari selisih waktu.

Selanjutnya, kita bisa membuat peta seismisitas sebagai berikut :

set F=peta.ps
set R=118/123.5/-4.7/0.8

REM $1 = Selisih waktu
REM $2 = Lintang
REM $3 = Bujur
REM $4 = Kedalaman
REM $5 = Magnitudo

REM Peta Utama
grdimage indo.nc -R%R% -JM5i -Celev.cpt -K -Ba1f1WesN -Y7 > %F%
gawk "{if ($5<=3) print $3, $2, $1, $5}" dataolah | psxy -R -Cwarna.cpt -JM -O -K -Sc0.1 -W0.1 >> %F%
gawk "{if ($5>3 && $5<=5) print $3, $2, $1, $5}" dataolah | psxy -R -Cwarna.cpt -JM -O -K -Sc0.2 -W0.1 >> %F%
gawk "{if ($5>5 && $5<6) print $3, $2, $1, $5}" dataolah | psxy -R -Cwarna.cpt -JM -O -K -Sc0.3 -W0.1 >> %F%
gawk "{if ($5>=6) print $3, $2, $1, $5}" dataolah | psxy -R -Cwarna.cpt -JM -O -K -Sc0.5 -W0.1 >> %F%

REM Cross section bawah
project -C90/0 -E101/0  -Q -G1 > track
gawk "{print $3, $2, $4, $5, $1}" dataolah | project -C118/0 -E123.5/0 -Fxyzpqrs > cross.dat

REM $1 = Bujur
REM $2 = Lintang
REM $3 = Kedalaman
REM $4 = Magnitudo
REM $5 = Selisih waktu

psbasemap -JX5i/1.5i -R118/123.5/-260/0 -Ba1:"Bujur":/:"Kedalaman (Km)":100WSne --FONT_LABEL=14p -K -O -Y-5 >> %F%
gawk "{if ($4<=3) print $1, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.1 -W0.1 -Cwarna.cpt -O -K >> %F%
gawk "{if ($4>3 && $4<=5) print $1, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.2 -W0.1 -Cwarna.cpt -O -K >> %F%
gawk "{if ($4>5 && $4<6) print $1, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.3 -W0.1 -Cwarna.cpt -O -K >> %F%
gawk "{if ($4>=6) print $1, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.5 -W0.1 -Cwarna.cpt -O -K >> %F%

REM Cross section samping
project -C120/-4.7 -E120/0.8  -Q -G1 > track
gawk "{print $3, $2, $4, $5, $1}" dataolah | project -C120/-4.7 -E120/0.8 -Fxyzpqrs > cross.dat

REM $1 = Bujur
REM $2 = Lintang
REM $3 = Kedalaman
REM $4 = Magnitudo
REM $5 = Selisih waktu

psbasemap -JX5i/1.5i -R-4.7/0.8/-260/0 -Ba1:"Lintang":/:"Kedalaman (Km)":100WSne --FONT_LABEL=14p -K -O -X14 -Y5 -p270/90 >> %F%
gawk "{if ($4<=3) print $2, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.1 -W0.1 -Cwarna.cpt -O -K -p >> %F%
gawk "{if ($4>3 && $4<=5) print $2, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.2 -W0.1 -Cwarna.cpt -O -K -p >> %F%
gawk "{if ($4>5 && $4<6) print $2, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.3 -W0.1 -Cwarna.cpt -O -K -p >> %F%
gawk "{if ($4>=6) print $2, $3*-1, $5, $4}" cross.dat | psxy -JX -R -Sc0.5 -W0.1 -Cwarna.cpt -O -K -p >> %F%

REM Legenda warna
psscale -D6/5/10/0.5 -Cwarna.cpt -O --FONT_ANNOT=14p -B120:"Hours": >> %F%

REM Cetak peta
psconvert -A -TG -P %F%

Selamat mencoba dan semoga sukses.

Salam.