Interpolasi Data Netcdf dengan Function interp (mpl_toolkits.basemap)

Jadi...apa motif dari pembuatan artikel ini? Having fun solving a problem, beberapa orang punya masalah karena kesulitan menemukan program untuk menginterpolasi data netcdf. Ya, ini sebenarnya selinganku belajar R sih. Metode belajarku sepertinya harus punya selingan lain supaya tidak bosan atau ini cuma karena procrastinationku yang sudah parah.

Mari kita mulai saja. Sebenarnya untuk mengecek hasil dari interpolasi ini sudah benar atau belum itu cukup sulit karena jika datanya di plot pada peta maka akan terlihat sama saja karena pada dasarnya interpolasi itu mengisi grid yang kosong saja bukan mengubah datanya.

Inti dari interpolasi ini adalah function interp() yang ada pada matplotlib (mpl_toolkits.basemap) a.k.a toolkits dari matplotlib...function-function tambahan untuk menunjang matplotlib sebagai library plotter. Hasil keluaran atau return value dari function ini adalah matriks 2 dimensi hasil interpolasi data input yang mana kerapatan interpolasinya itu di tentukan parameter "matriks koordinat outputnya". Kelihatannya runyam tapi akan lebih mudah jika langsung di beri contoh nyatanya dan yup, dokumentasi library yang baik harus punya banyak contoh programnya IMO dan matplotlib melakukan itu.

Seperti biasa, scriptnya sudah tersedia disini. Parameter-parameter yang dibutuhkan function interp() yaitu:
- matriks data input yang akan di interpolasi
- matriks dimensi input 1 dimensi (x dan y)
- matriks dimensi output 2 dimensi (x dan y, dibuat dengan bantuan np.meshgrid())

Algoritma dari program ini punya 2 tahap besar yaitu: interpolasi dan re-create nc file untuk wadah hasil interpolasi.

Untuk fase interpolasi, algoritmanya yaitu:

Open data input netcdf dengan mode read => ekstrak variabel input beserta dimensi-dimensinya => Menyiapkan matriks koordinat input (dapatkan dari object Basemap) => Menyiapkan matrik output => interpolasi => menyiapkan matriks lintang bujur yang sudah sesuai kerapatan interpolasi.

Lalu masuk fase re-create nc file-nya yaitu:

Open data netcdf dengan mode write (boleh merujuk ke data yang belum pernah ada karena ini mode write) => create variabel untuk data netcdf baru => set attribute (optional) => write data (dari matriks hasil interpolasi dan matrik-matriks lainnya dari algoritma fase interpolasi sebelumnya).

Di source code tersebut, bagian yang terlihat sepertinya arbitary yaitu

 #extract variable in data  
 plotvar = fh.variables[var][0, 710:1096, 1395:1772]     #composition: time, latitude, longitude  
 lat = fh.variables['latitude'][710:1096]     #grid 0,0 is bottom left and max,max is upper right  
 lon = fh.variables['longitude'][1395:1772]     #grid resolution is 0.02 degree

Ini dimaksudkan untuk mengambil dan menginterpolasi bagian pulau sulawesi saja. Data input awalnya mencakup wilayah Indonesia dan sekitarnya yang cukup luas. Angka-angka arbitary itu (710:1096, 1395:1772 dsb) didapatkan dengan memperhitungkan kalau jarak antar gridnya disini itu 0.02 dan pada koordinat 0,0 atau 0,0,0 (karena waktu pada data netcdf disini cuma punya 1 nilai maka bisa dianggap array atau matriks 2 dimensi) itu merujuk pada lintang dan bujur batas bawah pojok kiri dan sebaliknya (max, max, max merujuk pada batas atas pojok kanan). Ya, disini sampelnya adalah wilayah pulau sulawesi. Kenapa sulawesi? Saya lahir di sulawesi, nuff said.

Tampilan output dari data yang telah di interpolasi. Tidak jauh beda (atau mungkin sama) dengan output pada postingan sebelumnya pada wilayah yang sama.


Jadi, semisal ingin menginterpolasi semua wilayah pada data netcdf, tinggal mengubah kode di atas tadi menjadi

 #extract variable in data  
 plotvar = fh.variables[var][:]     #composition: time, latitude, longitude  
 lat = fh.variables['latitude'][:]     #grid 0,0 is bottom left and max,max is upper right  
 lon = fh.variables['longitude'][:]     #grid resolution is 0.02 degree

Dan tidak lupa juga mengubah koordinat pada function basemap() agar mencakup semua wilayah dari data netcdf dimana untuk itu kita perlu tahu mana batas bawah pojok kiri (lintang dan bujur) dan batas atas pojok kanannya (lintang dan bujur). Ini bisa diketahui dengan mengecek isi matriks dari variabel data netcdf di terminal python ataupun menjalankan kodenya.

 ~$ python  
 Python 2.7.12 (default, Nov 19 2016, 06:48:10)   
 [GCC 5.4.0 20160609] on linux2  
 Type "help", "copyright", "credits" or "license" for more information.  
 >>> import numpy as np  
 >>> from netCDF4 import Dataset  
 >>> dset = Dataset('/home/genomexyz/prjsatandro/tes1.nc', mode = 'r')  
 >>> dset.variables['latitude'][:]  
 array([-19.97999954, -19.95999908, -19.93999863, ..., 19.95999908,  
     19.97999954, 20.    ], dtype=float32)  
 >>> dset.variables['longitude'][:]  
 array([ 90.    ,  90.01999664,  90.04000092, ..., 149.94000244,  
     149.95999146, 149.97999573], dtype=float32)

Terlihat untuk batas bawah pojok kirinya, lintang dan bujurnya adalah -19.98o lintang dan 90o bujur dan batas atas pojok kanannya yaitu 20o lintang dan 149.98o bujur.

 m = Basemap(resolution='l', projection='merc', \  
 llcrnrlon=90, llcrnrlat=-19.98, urcrnrlon=149.98, urcrnrlat=20)


Bagian terpentingnya, pada script ini kerapatan hasil interpolasi bisa di atur pada function buatan yaitu create_grid_output(). Kerapatan di sini maksudnya, misal data netcdf awalnya jarak antar gridnya 10 km, Jika kita ingin jarak antara grid hasil interpolasi nanti 2 km, maka kerapatannya yang diinginkan adalah 5 kali lipat. Jadinya pada scriptnya function create_grid_output() menjadi

 #create grid output  
 x1, y1 = create_grid_output(x, y, 5)

Dengan mengatur matriks koordinat outpunya kita bisa mengatur kerapatannya.

Dan sekian artikel kali ini. Sisanya jika ingin lebih jelas silahkan lihat dokumentasi library matplotlib dan lain-lainnya yang dipakai pada script ini. Sekian.

referensi:
http://basemaptutorial.readthedocs.io/en/latest/utilities.html#interp, di akses pada 15 Maret 2017
https://matplotlib.org/basemap/users/mapcoords.html, di akses pada 15 Maret 2017

Komentar

Posting Komentar