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.
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%
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.
Data 27 September 2019??
ReplyDeleteApa tdk terlalu cepat om?
2018 Pak Bambs, salah ketik di narasinya :D
ReplyDeletecara menambahkan logo/gambar gmn ya
ReplyDeletepake data dengan format .ras (SUN RASTER)
Delete