Senin, 06 Januari 2020

Cara Plot dan Ekstrak Data GSMaP

Data GSMaP yang sudah didownload masih dalam format *.dat.gz, contoh file tanggal 01 Januari 2020 :  gsmap_nrt.20200101.0.1d.daily.00Z-23Z.dat.gz. Untuk itu file tersebut perlu diekstrak terlebih dahulu dengan menggunakan unzip sehingga nama filenya menjadi gsmap_nrt.20200101.0.1d.daily.00Z-23Z.dat.

Untuk plotting data diatas, dibutuhkan file deskriptor yang lazim juga disebut dengan kontrol file sehingga ekstensinya adalah *.ctl sebagai berikut :

DSET ^gsmap_nrt.%y4%m2%d2.0.1d.daily.00Z-23Z.dat
TITLE GSMaP_NRT 0.1deg Daily (00:00Z-23:59Z)
UNDEF -999.9
OPTIONS YREV LITTLE_ENDIAN TEMPLATE
XDEF   3600 LINEAR  0.05 0.1
YDEF   1200  LINEAR -59.95 0.1
zdef 1 levels 1013
tdef 365 linear 00:00z1jan2020 1dy
VARS 1
precip          0 99 daily averaged precip(mm/hr)
ENDVARS


untuk penjelasan dari tiap item dalam file deskriptor (*.ctl) diatas dapat dibaca : disini
 

Copy semua baris file *.ctl diatas dan Paste di Notepad, kemudian simpan dengan nama gsmap.ctl pada folder yang sama dengan folder penyimpanan data. 

Selanjutnya copy dan paste baris perintah GrADS berikut ke notepad dan beri nama filenya dengan nama plotgsmap.gs :

'reinit'

tanggal = 01 Januari 2020
'open PATH/gsmapctl.ctl'
'set t 1'
'set lon 95 141'
'set lat -15 15'
'set mpdset hires'
'set csmooth on'
'set gxout shaded'
'set clevs 0 5 20 50 100 150'
'd precip*24'
'cbarn'
'draw title Distribusi Curah Hujan GSMaP 01 Januari 2020'
'printim
PATH/01/01Januari2020.png white'

'set gxout print'
'set prnopts %6.2f 1 1'

write('
PATH/'tanggal'.txt', 'X    Y    LON    LAT    Precip')

'q dims'
xline=sublin(result,2)   ;* 2nd line
yline=sublin(result,3)   ;* 3rd line
xmax=subwrd(xline,13)    ;*13th word on xline
ymax=subwrd(yline,13)    ;*13th word on yline

say 'X grid-points: 'xmax
say 'Y grid-points: 'ymax
y=570.5
ymax=595.5
while(y<=ymax)

  x=1006.5
  xmax=1045.5
  while(x<=xmax)

    'set x 'x
    'set y 'y
    'd precip*24'

*    NOTE: It may be useful to test this to find out where the data is contained with in the result
*    It just so happens that in this case, the data is the 1st word of the 2nd line, this is not always true

     precip=sublin(result,2)
     precip=subwrd(precip,1)

*    Get Lat/Lon Data

     'q dims'
     lons=sublin(result,2)
     lats=sublin(result,3)
     lon=subwrd(lons,6)
     lat=subwrd(lats,6)

*    Save data to file
*    Note the "append", so to add to the file instead of overwriting it

     write('
PATH/'tanggal'.txt', x'    'y'    'lon'    'lat'    'precip,append)

     x=x+1
   endwhile
  y=y+1
endwhile 


####
Keterangan : PATH adalah alamat folder target penyimpanan file data dan file .ctl serta file .gs yang dibatasi oleh garis miring.


untuk menjalankan script GrADS diatas, masuk dulu ke aplikasi GrADS dengan klik dua kali pada icon GrADS. Maka akan tampil sebagai berikut :

Tekan ENTER maka akan tampil 

selanjutnya pada GrADS prompt, ketik run PATH/plotgsmap.gs maka akan tampil plot curah hujan tanggal 01 Januari 2020 wilayah Indonesia sekaligus menyimpan data tiap grid dalam format .txt kedalam folder target.  

Anda bisa juga melakukan plotting GSMaP yang telah tersimpan dalam format .txt dengan menggunakan aplikasi ArcGIS... https://arifmarufi.blogspot.com/2020/01/cara-plot-data-curah-hujan-harian-gsmap.html

37 komentar:

  1. Mau tanya, GSMaP kan datanya adanya harian, jika mau memetakan curah hujan bulanan pakai GSMaP sebaiknya bagaimana? Terima kasih

    BalasHapus
    Balasan
    1. Ada 2 cara jika ingin mendapatkan data bulanan :
      1. Bisa dengan menjumlahkan data harian dalam bulan yang sama
      2. Bisa juga dengan melalui alamat : ftp://hokusai.eorc.jaxa.jp/realtime/monthly/

      salam

      semoga berhasil

      Hapus
  2. Pak bagaimana ya jika saat di plot ke grads, hasilnya all undefined value?

    BalasHapus
    Balasan
    1. oh harusnya bisa.. mungkin bisa lebih detil bisa didiskusikan via email saya marufi.arif@yahoo.com

      Hapus
    2. sama saya juga seperti ini pak, izin email pak

      Hapus
  3. Ijin bertanya Pak, kenapa saat display itu precip*24? kenapa tidak d precip saja? apakah pengaruhnya karena plotting akumulasi curah hujan harian ya Pak? terima kasih Pak.

    BalasHapus
    Balasan
    1. itu saya kalikan 24 karena di metadata gsmap nya ditulis bahwa datanya average daily yang menurut saya adalah rata-rata harian. sehingga ketika kita akan melihat data total harian harus dikalikan 24 jam...
      pernah saya coba langsung d precip saja tapi nilainya sangat kecil sekali... terima kasih diskusinya... coba cek lagi metadata gsmap nya... untuk memastikan salam..

      Hapus
    2. bisa dibaca di file README.first.txt di alamat : ftp://hokusai.eorc.jaxa.jp/realtime/

      Hapus
  4. Sore mas Arif.
    saya jalankan plotgsmap nya, ada error
    "Unable to locate ENDWHILE statement for the WHILE statement at line 32"
    knapa ya??

    BalasHapus
    Balasan
    1. perlu diperhatikan identasi dari WHILE dan ENDWHILE nya, identasi adalah posisi pengetikan yang menjorok kedalam.. jika salah biasanya akan error...

      Hapus
    2. kalau masih error, boleh tinggalkan alamat emailnya nanti saya kirimkan file gs yang sudah benar identasinya... salam

      Hapus
    3. Terima Kasih pak Arif Ma'rufi.
      Setelah mentoring dengan bapak kurang lebih 2 jam secara online, tugas saya jadi lebih cepat. Sekaligus mendapatkan ilmu baru, terutama untuk alternatif mendapatkan data hujan di Indonesia.
      Salam hormat dan tetap jaga kesehatan.

      Hapus
    4. sama-sama pak, semoga bermanfaat. dan semoga kita semua selalu dilindungi Allah SWT.

      Salam Sehat

      Hapus
  5. kalau saya mau filter daerah daratn aja gimana ya caranya?

    BalasHapus
    Balasan
    1. ada dua alternatif :
      1. Menggunakan perintah wrseries, baca : https://arifmarufi.blogspot.com/2019/11/cara-download-dan-ekstrak-data-chirps.html

      2. Hasil ekstrak diatas dibuka di ArcGIS dan dilakukan proses clip pada geoprocessing wizard

      semoga bermanfaat

      salam

      Hapus
  6. Izin bertanya pak, ketika saya mencoba menjalankan script plotgsmap di atas, pada tampilan grads tertera 'unknown command : 1' dan setelah saya cek file txt nya ternyata kosong pak. Ini bagaimana pak? Mohon arahannya pak

    BalasHapus
    Balasan
    1. kemungkinan domain wilayahnya salah sehingga tidak berhasil diekstrak...

      kalau boleh tau, itu mau ekstrak wilayah mana ?

      Hapus
    2. Izin pak, saya rencananya akan mengekstrak wilayah Kabupaten Lebak, banten pak

      Hapus
    3. mungkin bisa dijelaskan di email saja mbak...
      email saya marufi.arif@yahoo.com

      salam

      Hapus
  7. Mohon maaf pak, untuk download GrADS nya dimana ya? saya mencoba download di internet tidak ada

    BalasHapus
  8. Izin bertanya Kak, kalau untuk per jamnya curah hujannya
    scritnya seperti apa Kak?

    BalasHapus
    Balasan
    1. prinsipnya sama saja sih.. tinggal download file nya yang jam2an.. dan di bagian d precip*24 diubah d precip saja

      Hapus
  9. Terima kasih mas Arif Makrufi.. blognya sangat bermanfaat.. terima kasih sudah memberi inspirasi..

    BalasHapus
    Balasan
    1. terima kasih supportnya pak, semoga bermanfaat..

      salam sehat

      Hapus
  10. Terima kasih tutorialnya. Jika ingin menghitung curah hujan selama 10 hari bagaimana? Apakah cukup dengan d ave(precip,t=1,t=10)*24?

    BalasHapus
    Balasan
    1. Kalau total 10 hari seharusnya menggunakan fungsi sum (jumlah) bukan ave (rata2)

      Hapus
  11. mohon izin bertanya pak, kok saya ada tulisan Open Error: Missing or invalid dimension starting value
    --> The invalid description file record is:
    --> TDEF 365 LINEAR 00:00Z9okt2022 1dy
    mksudnya gmana ya pak?

    BalasHapus
    Balasan
    1. TDEF 365 LINEAR 00:00Z9okt2022 1dy artinya ada 365 t (time), dimulai dari tanggal 09 oktober 2022.. tp menurut saya yang salah adalah penulisan singkatan bulan.. bukan okt tapi oct..

      jadi penulisan singkatan bulan harus dalam bahasa inggris

      semoga bermanfaat

      Hapus
  12. Komentar ini telah dihapus oleh pengarang.

    BalasHapus
  13. Pak, untuk nilai
    y=570.5
    ymax=595.5
    while(y<=ymax)

    x=1006.5
    xmax=1045.5
    while(x<=xmax)
    didapatkan dari mana ya? apakah area? klo saya cuma mau plot pulau sumatra saya, brp nilai yang harus sy inputkan?
    Terima kasih

    BalasHapus
    Balasan
    1. sepertinya itu hanya untuk area Jambi. itu adalah nomor urut grid X dan Y. Memang membatasi wilayah atau area. untuk menghitung suatu wilayah/area kita lihat dulu file ctl nya, disitu terlihat bahwa XDEF 0.05 0.1 yang berarti start di 0.05 derajat dan ukuran grid nya 0.1 derajat sehingga total gridnya menjadi 3600 grid (karena bujur data yang tersedia 360 derajat). sedangkan YDEF -59.95 0.1 yang berarti start di -59.95 derajat dengan ukuran grid 0.1 derajat sehingga total grid nya menjadi 1200 (karena total lintang yang tersedia dalam data tsb adalah 120 derajat).

      misalnya kita ingin mengambil wilayah dengan koordinat 90 derajat sd 100 derajat, maka kita hitung mulai dari 0.05 derajat dengan penambahan 0.1 derajat sampai ketemu 90 derajat, nomor barisnya adalah Xawal dan Xakhir adalah nomor baris pada 100 derajat.

      demikian pula cara menghitung nomor grid Y nya.. tp start dari -59.95 + 0.1 sesuai dengan panduan file ctl.

      Hapus
  14. cara mendapatkan x dan y itu bagaimana ya pak

    BalasHapus
  15. siang pak, izin bertanya saya juga mengolah data gsmap mau diubah dulu ke nc dan mau coba olah pakai python, tapi saat udat dat ke nc saya mengalami kendala, tertulis file notfound bagaimana cara memperbaikinya ya pak?

    BalasHapus
    Balasan
    1. wah kalau itu perlu diurut dari awal lagi.. saran saya didownload ulang saja file nya

      Hapus