Jawa Barat merupakan salah satu wilayah yang perlu mendapat perhatian dalam kajian kegempaan karena memiliki riwayat kejadian gempa bumi yang tersebar pada beberapa lokasi. Informasi gempa bumi yang dipublikasikan BMKG (2026) memuat parameter seperti waktu kejadian, magnitudo, kedalaman, koordinat, dan wilayah terdampak, sehingga data tersebut dapat digunakan untuk melihat sebaran kejadian gempa secara spasial. Dalam kajian bahaya gempa, dampak guncangan tidak hanya berkaitan dengan besarnya magnitudo atau kedalaman pusat gempa, tetapi juga dipengaruhi oleh kondisi tanah di lokasi tertentu. Oleh karena itu, analisis kondisi tanah menjadi penting untuk memahami wilayah mana yang secara karakteristik tapak lebih rentan terhadap penguatan guncangan.
Dalam penelitian ini, data gempa bumi dibatasi pada kejadian dengan magnitudo minimal 4 mb. Pembatasan ini dilakukan agar analisis tidak terlalu dipengaruhi oleh kejadian gempa kecil yang dampaknya relatif rendah, tetapi lebih difokuskan pada gempa yang lebih berpotensi menimbulkan guncangan terasa dan perlu diperhatikan dalam kajian kerentanan kondisi tanah. BMKG Stasiun Geofisika Kelas II Pasuruan (2026) menjelaskan bahwa pada skala IV MMI, getaran dapat dirasakan oleh banyak orang di dalam rumah, sedangkan pada skala V MMI getaran dapat dirasakan hampir semua penduduk dan barang-barang dapat terpelanting. Dengan demikian, penggunaan batas magnitudo ≥ 4 mb bertujuan untuk menyeleksi kejadian gempa yang lebih relevan terhadap kajian guncangan, sementara data gempa tetap digunakan sebagai overlay dan bukan sebagai variabel prediksi.
Salah satu indikator yang umum digunakan untuk menggambarkan karakteristik tanah adalah Vs30, yaitu rata-rata kecepatan gelombang geser hingga kedalaman 30 meter. USGS (2025) menjelaskan bahwa Vs30 digunakan sebagai indikator kondisi tapak seismik karena dapat membedakan karakteristik tanah lunak, tanah sedang, tanah keras, hingga batuan. Nilai Vs30 yang lebih rendah umumnya menunjukkan tanah yang lebih lunak, sedangkan nilai Vs30 yang lebih tinggi menunjukkan tanah yang lebih keras. Dalam penelitian ini, kerentanan yang dimaksud bukan kerentanan sosial atau kerusakan bangunan, melainkan kerentanan kondisi tanah terhadap guncangan gempa berdasarkan karakteristik Vs30. Dengan demikian, data gempa bumi digunakan sebagai informasi overlay, sedangkan variabel yang dimodelkan secara spasial adalah nilai Vs30.
Pendekatan spasial diperlukan karena nilai Vs30 tidak berdiri sendiri sebagai angka pada satu titik, tetapi memiliki pola kedekatan antarwilayah. Lokasi yang berdekatan cenderung memiliki karakteristik tanah yang lebih mirip dibandingkan lokasi yang berjauhan, sehingga metode interpolasi spasial dapat digunakan untuk memperkirakan nilai Vs30 pada wilayah yang tidak memiliki titik pengamatan langsung. Esri (2025) menjelaskan bahwa Kriging merupakan metode geostatistik yang menghasilkan permukaan prediksi dari sekumpulan titik pengamatan dengan mempertimbangkan struktur spasial data melalui semivariogram. Penelitian ini menggunakan Universal Kriging karena hasil pemeriksaan awal menunjukkan adanya kecenderungan tren spasial pada nilai Vs30 berdasarkan koordinat lokasi. Dengan pendekatan tersebut, penelitian ini diharapkan dapat menghasilkan peta prediksi Vs30 dan peta ketidakpastian prediksi sebagai dasar untuk memahami variasi spasial kerentanan kondisi tanah terhadap guncangan gempa bumi di Jawa Barat.
Berdasarkan latar belakang, dapat diidentifikasi beberapa permasalahan utama dalam penelitian ini, yaitu:
Berdasarkan tujuan penelitian, pertanyaan prediksi pada penelitian ini ialah:
Data penelitian ini terdiri dari data Vs30, data kejadian gempa bumi, dan data batas administrasi wilayah Jawa Barat. Data Vs30 diperoleh dari USGS Global Vs30 Mosaic Map Viewer. USGS menjelaskan bahwa Global Vs30 Mosaic menyediakan informasi Vs30, yaitu rata-rata kecepatan gelombang geser pada kedalaman 30 meter, dan data tersebut dapat diakses melalui aplikasi peta Vs30 USGS. Data gempa bumi diperoleh dari katalog United States Geological Survey (USGS) melalui pustaka ObsPy dengan layanan FDSN Event Web Service. USGS menjelaskan bahwa layanan FDSN Event Web Service dapat digunakan untuk melakukan pencarian informasi gempa berdasarkan parameter tertentu, seperti waktu, lokasi, magnitudo, dan kedalaman kejadian. Data gempa pada penelitian ini diambil untuk wilayah sekitar Jawa Barat dengan batas latitude -8 sampai -5 dan longitude 105 sampai 109 pada periode 2000 sampai 2026. Data batas administrasi diperoleh dari Geoportal Badan Informasi Geospasial (BIG) berupa batas administrasi desa/kelurahan Jawa Barat, yang digunakan sebagai pembatas wilayah analisis dan dasar visualisasi peta.
Unit observasi utama dalam penelitian ini adalah titik lokasi Vs30 di
Jawa Barat. Setiap titik memiliki informasi koordinat longitude dan
latitude serta nilai vs30_mosaic. Variabel
vs30_mosaic digunakan sebagai variabel utama yang
diprediksi menggunakan Universal Kriging. Koordinat longitude dan
latitude digunakan untuk membentuk objek spasial, kemudian
ditransformasikan ke sistem koordinat UTM sehingga proses analisis jarak
dapat dilakukan dalam satuan meter. Pada Universal Kriging, koordinat
hasil transformasi digunakan sebagai variabel x_utm dan
y_utm untuk mewakili tren spasial. Data gempa bumi
digunakan sebagai data pendukung dengan variabel longitude, latitude,
magnitude, depth_km, tahun, bulan, tanggal, dan kota. Data gempa tidak
digunakan sebagai variabel prediksi, tetapi ditampilkan sebagai overlay
pada peta hasil prediksi Vs30 untuk melihat posisi kejadian gempa
terhadap karakteristik kondisi tanah.
Metode yang digunakan dalam penelitian ini adalah Universal Kriging.
Metode ini dipilih karena nilai Vs30 memiliki unsur lokasi dan
menunjukkan indikasi tren spasial berdasarkan koordinat. Esri (2025)
menjelaskan bahwa Universal Kriging digunakan ketika data memiliki
komponen tren yang berubah menurut lokasi. Dalam penelitian ini,
variabel yang diprediksi adalah vs30_mosaic, sedangkan
koordinat hasil transformasi UTM, yaitu x_utm dan
y_utm, digunakan sebagai komponen tren spasial. Dengan
mengunakan parameter \(\beta\), Model
Universal Kriging dituliskan sebagai:
\[ Z(s) = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon(s) \]
Keterangan:
\(Z(s)\) = nilai Vs30 pada lokasi (s), (x) dan (y) adalah koordinat UTM
\(\varepsilon(s)\) = residual spasial
Pebesma (2025) menjelaskan bahwa dalam package gstat, Universal
Kriging dapat dijalankan dengan formula seperti z ~ x + y,
sedangkan Ordinary Kriging menggunakan formula z ~ 1. Oleh
karena itu, model pada penelitian ini dijalankan dengan formula
vs30_mosaic ~ x_utm + y_utm.
Struktur spasial residual dimodelkan menggunakan semivariogram. SciKit-GStat (2025) menjelaskan bahwa variogram menggambarkan hubungan antara jarak antar titik dengan tingkat perbedaan nilai observasi. Semivariogram empiris pada penelitian ini dicocokkan dengan beberapa model teoritis, yaitu Spherical, Exponential, dan Gaussian. Model variogram terbaik dipilih berdasarkan nilai error fitting yang paling kecil dan bentuk kurva yang paling mengikuti semivariogram empiris. Dalam Kriging, bobot prediksi tidak ditentukan secara manual, tetapi dihitung dari model variogram. Esri (2025) menjelaskan bahwa Kriging menggunakan bobot dari semivariogram untuk memprediksi nilai pada lokasi yang tidak memiliki pengamatan.
Prediksi nilai Vs30 dilakukan pada grid wilayah Jawa Barat. Hasil utama dari model adalah nilai prediksi Vs30 dan nilai ketidakpastian prediksi. Esri (2025) menjelaskan bahwa nilai variansi prediksi yang lebih rendah menunjukkan tingkat keyakinan prediksi yang lebih tinggi, sedangkan nilai variansi yang lebih tinggi menunjukkan ketidakpastian yang lebih besar. Evaluasi model dilakukan menggunakan cross-validation dengan ukuran ME, MAE, RMSE, R² CV, dan korelasi antara nilai observasi dan prediksi. Dokumentasi gstat menjelaskan bahwa fungsi krige.cv dapat digunakan untuk validasi silang pada simple, ordinary, maupun universal kriging (Pebesma, 2025). Oada penelitian ini, informasi waktu pada data gempa hanya digunakan sebagai informasi deskriptif, sedangkan model utama difokuskan pada pemetaan nilai Vs30 berdasarkan lokasi.
Data Vs30 yang digunakan dalam penelitian ini berjumlah 42.727 titik setelah proses pembersihan data. Nilai Vs30 minimum sebesar 180 m/s dan maksimum sebesar 900 m/s, dengan rata-rata sebesar 450,01 m/s dan median sebesar 437 m/s. Nilai standar deviasi sebesar 192,68 menunjukkan bahwa nilai Vs30 di Jawa Barat memiliki variasi yang cukup besar antarwilayah. Berdasarkan rentang nilainya, wilayah kajian tidak hanya didominasi oleh satu jenis karakteristik tanah, tetapi mencakup kondisi tanah yang relatif lebih lunak hingga lebih keras.
| Data | n | Minimum | Rata-rata | Median | Maksimum | Standar Deviasi |
|---|---|---|---|---|---|---|
| Vs30 | 42.727 | 180 | 450,01 | 437 | 900 | 192,68 |
| Magnitudo gempa | 224 | 4 | 4,49 | 4,4 | 6,5 | - |
| Kedalaman gempa | 224 | 1,41 | 85,67 | 78,1 | 329,6 | - |
Data kejadian gempa bumi yang lebih dari 4 magnitudo berjumlah 224 kejadian. Magnitudo gempa berada pada rentang 4 sampai 6,5, dengan rata-rata magnitudo sebesar 4,49 dan median sebesar 4,4. Kedalaman gempa berkisar antara 1,41 km sampai 329,6 km, dengan rata-rata kedalaman sebesar 85,67 km. Hal ini menunjukkan bahwa data gempa yang digunakan mencakup kejadian dengan kedalaman dangkal hingga dalam. Dalam penelitian ini, data gempa tidak digunakan sebagai variabel yang diprediksi, tetapi digunakan sebagai overlay untuk melihat posisi kejadian gempa terhadap hasil pemetaan Vs30.
Gambar 3.1. Peta Sebaran Titik Vs30 di Jawa Barat
Peta sebaran nilai Vs30 menunjukkan bahwa kondisi tanah di Jawa Barat bervariasi secara spasial. Warna yang lebih terang menunjukkan nilai Vs30 lebih rendah, sedangkan warna yang lebih gelap menunjukkan nilai Vs30 lebih tinggi. Hal ini menunjukkan bahwa wilayah Jawa Barat memiliki karakteristik tanah yang berbeda-beda, dari tanah yang relatif lebih lunak hingga lebih keras.
Gambar 3.2. Peta Kontur Density Gempa Bumi Dengan Minimal 4 mb di Jawa Barat
Peta kontur density gempa menunjukkan bahwa kejadian gempa bumi yang lebih dari sama dengan 4 mb tidak tersebar merata. Kepadatan kejadian gempa terlihat lebih tinggi pada beberapa area, terutama bagian barat dan sebagian selatan Jawa Barat. Kontur yang lebih rapat menunjukkan konsentrasi kejadian gempa yang lebih tinggi.
Pemeriksaan pola spasial dilakukan untuk melihat apakah data Vs30 memiliki ketergantungan spasial dan pola tren berdasarkan lokasi. Hasil Moran’s I menunjukkan nilai sebesar 0,7808. Nilai tersebut menunjukkan adanya autokorelasi spasial positif yang kuat, sehingga titik-titik yang saling berdekatan cenderung memiliki nilai Vs30 yang relatif mirip. Kondisi ini mendukung penggunaan metode Kriging karena metode tersebut mempertimbangkan hubungan spasial antar titik pengamatan dalam proses prediksi.
Gambar 3.3. Grafik Tren Vs30 terhadap Koordinat X
Gambar 3.4. Grafik Tren Vs30 terhadap Koordinat Y
Selain autokorelasi spasial, pola tren nilai Vs30 juga diperiksa
secara visual melalui grafik hubungan Vs30 terhadap koordinat X dan Y
UTM. Pada grafik terhadap koordinat X, garis tren tidak sepenuhnya datar
dan menunjukkan perubahan nilai Vs30 dari arah barat ke timur. Sementara
itu, pada grafik terhadap koordinat Y, pola tren terlihat lebih jelas,
yaitu nilai Vs30 cenderung meningkat pada bagian tengah wilayah kajian
kemudian menurun pada koordinat Y yang lebih tinggi. Pola tersebut
menunjukkan bahwa nilai rata-rata Vs30 tidak sepenuhnya konstan di
seluruh wilayah kajian. Dengan adanya kecenderungan tren spasial
tersebut menjadi dasar penggunaan Universal Kriging. Metode ini sesuai
karena dapat mempertimbangkan komponen tren berdasarkan koordinat
lokasi, yaitu x_utm dan y_utm, sekaligus tetap
memodelkan struktur spasial residual melalui variogram.
Gambar 3.5. Histogram Jarak Nearest Neighbor Titik Vs30
Kerapatan titik pengamatan juga diperiksa melalui histogram jarak nearest neighbor. Hasil pemeriksaan menunjukkan bahwa sebagian besar titik Vs30 memiliki jarak yang relatif dekat dengan titik terdekatnya, terutama pada kisaran kurang dari 2.000 meter. Hal ini menunjukkan bahwa data pengamatan cukup rapat untuk mendukung proses interpolasi spasial. Namun, keberadaan beberapa titik dengan jarak yang lebih jauh tetap perlu diperhatikan, sehingga hasil prediksi sebaiknya dibaca bersama peta ketidakpastian.
Pemodelan variogram dilakukan untuk menggambarkan struktur spasial residual pada Universal Kriging. Pada penelitian ini, tiga model variogram teoritis dibandingkan, yaitu Spherical, Exponential, dan Gaussian.
| Model | SSErr | Nugget | Partial Sill | Total Sill | Range |
|---|---|---|---|---|---|
| Exponential | 5.376.743,52 | 6.557,63 | 26.344,10 | 32.901,73 | 21.164,40 |
| Spherical | 11.624.913,28 | 10.849,90 | 21.313,42 | 32.163,32 | 58.438,94 |
| Gaussian | 19.803.820,04 | 12.542,58 | 19.228,07 | 31.770,65 | 24.661,09 |
Berdasarkan hasil fitting, model Exponential memiliki nilai SSErr paling kecil, yaitu 5.376.744. Nilai tersebut lebih kecil dibandingkan model Spherical sebesar 11.624.913 dan Gaussian sebesar 19.803.820. Oleh karena itu, model Exponential dipilih sebagai model variogram terbaik untuk proses Universal Kriging.
Gambar 3.6. Perbandingan Model Variogram Spherical, Exponential, dan Gaussian
Model Exponential memiliki nugget sebesar 6.557,6, partial sill sebesar 26.344,1, total sill sebesar 32.901,73, dan range sebesar 21.164,4 meter. Nugget menunjukkan variasi pada jarak yang sangat dekat atau variasi yang tidak tertangkap oleh model. Partial sill menunjukkan variasi spasial yang dapat dijelaskan oleh struktur spasial, sedangkan total sill menunjukkan keseluruhan variasi yang terbentuk. Range sebesar sekitar 21,16 km menunjukkan bahwa hubungan spasial antar titik Vs30 masih cukup kuat sampai jarak tersebut. Setelah melewati jarak tersebut, keterkaitan spasial antar titik cenderung semakin melemah.
Prediksi nilai Vs30 dilakukan pada 1.485 titik grid di wilayah Jawa Barat menggunakan Universal Kriging.
| N | ME | MAE | RMSE | R² CV | Korelasi Observasi-Prediksi |
|---|---|---|---|---|---|
| 2.000 | 0,965 | 72,19 | 102,20 | 0,724 | 0,851 |
Berdasarkan tabel tersebut, terlihat bahwa nilai ME yang mendekati nol menunjukkan bahwa model tidak memiliki bias prediksi yang besar. Nilai korelasi sebesar 0,851 menunjukkan hubungan yang kuat antara nilai Vs30 hasil observasi dan hasil prediksi. Dengan demikian, model Universal Kriging dapat dikatakan cukup baik dalam memetakan pola spasial Vs30 di wilayah kajian.
Gambar 3.7. Peta Prediksi Vs30 di Jawa Barat
Peta prediksi Vs30 menunjukkan bahwa bagian utara Jawa Barat cenderung memiliki nilai Vs30 yang lebih rendah dibandingkan beberapa wilayah tengah dan selatan. Nilai Vs30 yang lebih rendah menunjukkan karakteristik tanah yang relatif lebih lunak, sedangkan nilai Vs30 yang lebih tinggi menunjukkan tanah yang relatif lebih keras. Berdasarkan hasil klasifikasi grid prediksi, sebagian besar wilayah berada pada kelas SC atau tanah keras, yaitu sebanyak 958 grid atau 64,51%. Kelas SD atau tanah sedang/lunak mencakup 470 grid atau 31,65%. Sementara itu, kelas SB atau batuan mencakup 53 grid atau 3,57%, dan kelas SE atau tanah sangat lunak hanya mencakup 4 grid atau 0,27%.
Gambar 3.8. Peta Ketidakpastian Prediksi Vs30
Peta ketidakpastian prediksi menunjukkan bahwa sebagian besar wilayah memiliki nilai ketidakpastian yang relatif rendah. Namun, terdapat area dengan ketidakpastian lebih tinggi, terutama pada bagian tepi wilayah. Hal ini wajar karena area tepi biasanya memiliki jumlah titik pengamatan sekitar yang lebih terbatas dibandingkan area tengah.
Overlay kejadian gempa dilakukan untuk melihat posisi titik gempa terhadap peta prediksi Vs30.
Gambar 3.9. Overlay Kejadian Gempa Bumi Minimal 4 mb terhadap Peta Prediksi Vs30
Hasil overlay menunjukkan bahwa sebagian besar titik gempa berada pada wilayah dengan kelas SC atau tanah keras, yaitu sebanyak 177 kejadian atau 79,1%. Sebanyak 37 kejadian atau 16,5% berada pada kelas SD atau tanah sedang/lunak. terakhir, terdapat 9 kejadian atau 4.02% berada pada kelas SB atau batuan dan 1 kejadian atau 0,45% berada pada kelas SE atau tanah sangat lunak.
Hasil ini menunjukkan bahwa lokasi kejadian gempa pada data penelitian tidak hanya berada pada wilayah dengan Vs30 rendah, tetapi tersebar pada beberapa kelas kondisi tanah. Namun, dari sisi kerentanan kondisi tanah, wilayah dengan Vs30 yang lebih rendah tetap perlu diperhatikan karena tanah yang lebih lunak secara konsep lebih berpotensi mengalami penguatan guncangan dibandingkan tanah yang lebih keras.
Secara substantif, hasil penelitian menunjukkan bahwa kerentanan kondisi tanah terhadap guncangan gempa bumi di Jawa Barat bervariasi antarwilayah. Wilayah dengan nilai Vs30 lebih rendah menunjukkan karakteristik tanah yang relatif lebih lunak, sehingga lebih perlu diperhatikan karena berpotensi mengalami penguatan guncangan. Sebaliknya, wilayah dengan nilai Vs30 lebih tinggi menunjukkan kondisi tanah yang relatif lebih keras.
Penggunaan Universal Kriging didukung oleh adanya autokorelasi spasial positif pada nilai Vs30 serta indikasi tren berdasarkan koordinat lokasi. Model Exponential menjadi model variogram terbaik karena memiliki nilai SSErr paling kecil. Hasil cross-validation juga menunjukkan performa model yang cukup baik, sehingga peta prediksi Vs30 dapat digunakan sebagai gambaran awal variasi kondisi tanah di Jawa Barat.
Namun, penelitian ini memiliki beberapa keterbatasan. Variabel yang dimodelkan hanya Vs30, sehingga hasil analisis hanya menggambarkan kerentanan dari sisi kondisi tanah. Data gempa magnitudo ≥ 4 mb hanya digunakan sebagai overlay dan kontur density, bukan sebagai variabel prediksi. Oleh karena itu, penelitian ini tidak memprediksi kejadian gempa, magnitudo, PGA, maupun tingkat kerusakan. Selain itu, hasil prediksi tetap dipengaruhi oleh kerapatan titik Vs30, sehingga wilayah dengan ketidakpastian tinggi perlu diinterpretasikan secara hati-hati.
Berdasarkan hasil analisis, nilai Vs30 di Jawa Barat menunjukkan variasi spasial yang cukup jelas. Hasil Moran’s I sebesar 0,7808 menunjukkan adanya autokorelasi spasial positif yang kuat, sedangkan grafik tren terhadap koordinat X dan Y menunjukkan adanya kecenderungan perubahan nilai Vs30 berdasarkan lokasi. Kondisi ini mendukung penggunaan Universal Kriging dalam pemetaan nilai Vs30. Model variogram terbaik adalah Exponential karena memiliki nilai SSErr paling kecil dibandingkan model Spherical dan Gaussian. Hasil cross-validation menunjukkan performa model yang cukup baik, dengan RMSE sebesar 102,20, R² CV sebesar 0,724, dan korelasi observasi-prediksi sebesar 0,851. Peta prediksi menunjukkan bahwa sebagian besar wilayah berada pada kelas SC atau tanah keras, tetapi masih terdapat wilayah kelas SD atau tanah sedang/lunak yang perlu diperhatikan karena secara karakteristik lebih rentan terhadap penguatan guncangan.
Overlay kejadian gempa magnitudo ≥ 4 mb menunjukkan bahwa titik gempa tersebar pada beberapa kelas Vs30, dengan sebagian besar berada pada kelas SC. Dengan demikian, hasil penelitian ini dapat digunakan sebagai gambaran awal kerentanan kondisi tanah terhadap guncangan gempa bumi berdasarkan nilai Vs30. Namun, hasil ini tidak digunakan untuk memprediksi kejadian gempa, magnitudo, PGA, atau tingkat kerusakan.
Badan Meteorologi, Klimatologi, dan Geofisika. (2026). Gempa bumi dirasakan. Diambil dari https://www.bmkg.go.id/gempabumi/gempabumi-dirasakan
BMKG Stasiun Geofisika Kelas II Pasuruan. (2026). Skala MMI (Modified Mercalli Intensity). Diambil dari https://stageof-tretes.bmkg.go.id/skala-mmi
Esri. (2025a). How kriging works. ArcGIS Pro Documentation. https://doc.esri.com/en/arcgis-pro/latest/tool-reference/spatial-analyst/how-kriging-works.html
Esri. (2025b). Kriging (Spatial Analyst tools). ArcGIS Pro Documentation. Diambil dari https://doc.esri.com/en/arcgis-pro/latest/tool-reference/spatial-analyst/kriging.html?tabs=dialog
Esri. (2025c). Kriging in Geostatistical Analyst. ArcGIS Pro Documentation. Diambil dari https://doc.esri.com/en/arcgis-pro/latest/help/analysis/geostatistical-analyst/kriging-in-geostatistical-analyst.html
Esri. (2025d). Understanding universal kriging. ArcGIS Pro Documentation. Diambil dari https://doc.esri.com/en/arcgis-pro/latest/help/analysis/geostatistical-analyst/understanding-universal-kriging.html
Pebesma, E. (2025). gstat: Spatial and spatio-temporal geostatistical modelling, prediction and simulation. CRAN. Diambil dari https://cran.r-project.org/web/packages/gstat/gstat.pdf
SciKit-GStat. (2025). Variography. SciKit-GStat Documentation. Diambil dari https://scikit-gstat.readthedocs.io/en/latest/userguide/variogram.html
U.S. Geological Survey. (2025). Vs30 models and data. USGS Earthquake Hazards Program. Diambil dari https://earthquake.usgs.gov/data/vs30/