Ekstrak Data Radar GEMATRONIK (Single Pol)

Halo sobat radar cuaca! Siapa yang tau radar GEMATRONIK?

Radar Gematronik ini adalah radar cuaca pabrikan Jerman, sebelum tahun 2019 nama pabriknya Selex ES GmbH, namun saat ini sudah berganti nama menjadi Leonardo.

Radar Gematronik ini cukup spesial, kalau radar EEC dan BARON hanya ada satu rak kabinet, radar Gematronik sampai 3 rak kabinet! Bayangin tuh betapa terlihat kokohnya sistem radar ini (pendapat pribadi ya). Struktur datanya juga cukup unik, jika radar lain menyatukan seluruh moment (Z, V, W) menjadi satu file, radar gematronik memisahkan dalam file yang berbeda. Ada plus minusnya sih, salah satu plusnya adalah ukuran data yang tidak besar, jadi memudahkan proses transfer data apabila data tersebut harus didistribusikan ke beberapa lokasi.

Minusnya, dalam pengolahan data menggunakan python jadi sedikit ribet, apalagi kalau si radar mengaplikasikan jenis filter tertentu, sistem radar akan menulis file baru di lokasi folder lain, namun nama filenya masih sama seperti rawdatanya. Cukup ribet, tapi masih bisa diakalin.

OK, sebelum kita mulai, seperti biasa ada beberapa kondisi yang sudah harus sobat penuhi untuk bisa mengolah data radar ini :

1. Telah menginstall python (baca juga : Tutorial Instalasi Python Menggunakan Anaconda)
2. Telah menginstall modul wradlib (baca juga : Tutorial Instalasi Modul Wradlib)
3. Memahami 3 langkah utama dalam melakukan ekstraksi data radar cuaca (baca bagian awal-awal pada postingan ini : Ekstrak Data Radar BARON)

Sebelum kalian membuka editor, apabila kalian menginstall modul wradlib pada environtment tertentu, jangan lupa diaktifkan dulu environmentnya ya!

Bersiaplah, ini tidak akan mudah, tapi saat kamu sudah berhasil, ini akan menjadi sangat menyenangkan!

Ali, 2020

Import modul-modul yang diperlukan

from datetime import datetime
from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt
import wradlib as wrl
import numpy as np
import warnings, math
warnings.filterwarnings("ignore")
warnings.filterwarnings("ignore", category=DeprecationWarning) 
warnings.filterwarnings("ignore", category=RuntimeWarning)

Tentukan lokasi data, data sampel dapat kalian download di link ini.

radarFile='D:/project_webprogramming/wxradarexplore/radarDataExtraction/data/BKS20180516191000dBZ.vol'
f = wrl.util.get_wradlib_data_file(radarFile)
raw = wrl.io.read_rainbow(f)

Dari potongan script diatas, fungsi utama untuk melakukan esktrak data radar gematronik adalah wrl.io.read_rainbow. Apa itu rainbow? Rainbow adalah nama aplikasi pengolahan data radar gematronik. Nah berbeda dengan radar BARON, fungsi wrl.io.read_rainbow hanya menghasilkan satu variabel output, admin namakan raw. Masih ingat variabel dictionary (dict) kan? Ya, admin menamakan variabel keluarannya sebagai rbdict karena memang outputnya berupa variabel tipe dictionary.

OK mari kita lihat isi raw (jangan lupa untuk running script biar variabelnya muncul ya!).

Masih ingat 5 parameter kunci untuk mengekstrak data radar? Yup salah satunya lokasi radar. Nah dalam data yang admin pakai, lokasi radar terdapat pada raw[‘volume’][‘sensorinfo’].

Tapi kalian harus memastikan lagi pada data yang kalian pakai, karena bisa saja berbeda bergantung pada generasinya. Jadi disini admin akan memberikan semua kemungkinan pada data rainbow (dengan try-except) OK mari kita ekstrak lokasi radarnya.

try:
    radarLon=float(raw['volume']['sensorinfo']['lon'])
    radarLat=float(raw['volume']['sensorinfo']['lat'])
    radarAlt=float(raw['volume']['sensorinfo']['alt'])
except:
    radarLon=float(raw['volume']['radarinfo']['@lon'])
    radarLat=float(raw['volume']['radarinfo']['@lat'])
    radarAlt=float(raw['volume']['radarinfo']['@alt'])
sitecoords=(radarLon,radarLat,radarAlt)

Langkah selanjutnya sama dengan postingan ekstrak data radar BARON, yaitu menyiapkan container data yang sudah di regriding.

res=250. # resolusi data yang diinginkan dalam meter
resCoords=res/111229. # resolusi data dalam derajat
rmax=250000./111229. # range maksimum
lonMax,lonMin=radarLon+(rmax),radarLon-(rmax) 
latMax,latMin=radarLat+(rmax),radarLat-(rmax)
nGrid=int(np.floor((lonMax-lonMin)/resCoords))+1 # jumlah grid
lonGrid=np.linspace(lonMin,lonMax,nGrid) # grid longitude
latGrid=np.linspace(latMin,latMax,nGrid) # grid latitude          
dataContainer = np.zeros((len(lonGrid),len(latGrid))) # penampung data

Selanjutnya something challenging! Yaitu mengekstrak azimuth, range, data radar pada setiap elevasi. Hal ini agak sedikit lebih membutuhkan pemahaman, anyway nanti admin jelaskan. Ini codenya :

# menentukan waktu (end of observation)
nSlices=len(raw['volume']['scan']['slice'])
date=(raw['volume']['scan']['slice'][nSlices-1]['slicedata']['@date'])
time=(raw['volume']['scan']['slice'][nSlices-1]['slicedata']['@time'])
try:timeEnd=datetime.strptime('{}{}'.format(date,time),"%Y-%m-%d%H:%M:%S")
except:timeEnd=datetime.strptime('{}{}'.format(date,time),"%Y-%m-%d%H:%M:%S.%f")

nElevation=len(raw['volume']['scan']['slice']) # jumlah seluruh elevasi
for i in range(nElevation):
    try:elevation = float(raw['volume']['scan']['slice'][i]['posangle'])
    except:elevation = float(raw['volume']['scan']['slice'][0]['posangle'])

    print('Extracting radar data : SWEEP-{0} at Elevation Angle {1:.1f} deg ...'.format(i+1,elevation))
    
    # ekstrak azimuth data
    try:
        azi = raw['volume']['scan']['slice'][i]['slicedata']['rayinfo']['data']
        azidepth = float(raw['volume']['scan']['slice'][i]['slicedata']['rayinfo']['@depth'])
        azirange = float(raw['volume']['scan']['slice'][i]['slicedata']['rayinfo']['@rays'])  
    except:
        azi0 = raw['volume']['scan']['slice'][i]['slicedata']['rayinfo'][0]['data']
        azi1 = raw['volume']['scan']['slice'][i]['slicedata']['rayinfo'][1]['data']
        azi = (azi0/2) + (azi1/2)
        del azi0, azi1
        azidepth = float(raw['volume']['scan']['slice'][i]['slicedata']['rayinfo'][0]['@depth'])
        azirange = float(raw['volume']['scan']['slice'][i]['slicedata']['rayinfo'][0]['@rays'])            
    try:
        azires = float(raw['volume']['scan']['slice'][i]['anglestep'])
    except:
        azires = float(raw['volume']['scan']['slice'][0]['anglestep'])
    azi = (azi * azirange / 2**azidepth) * azires
    
    flag=0
    if np.size(azi) >= 999:
        flag=2
        azi = azi/3
        for ii in range(int(np.floor(np.size(azi)/3))):
            azi[ii] = azi[3*ii]+azi[3*ii+1]+azi[3*ii+2]
        azi = azi[range(int(np.floor(np.size(azi)/3)))]
    elif np.size(azi) >= 500:
        flag=1
        azi = azi/2
        for ii in range(int(np.floor(np.size(azi)/2))):
            azi[ii] = azi[2*ii]+azi[2*ii+1]
        azi = azi[range(int(np.floor(np.size(azi)/2)))]
    
    # esktrak range data
    try:
        stoprange = float(raw['volume']['scan']['slice'][i]['stoprange'])
        rangestep = float(raw['volume']['scan']['slice'][i]['rangestep'])
    except:
        stoprange = float(raw['volume']['scan']['slice'][0]['stoprange'])
        rangestep = float(raw['volume']['scan']['slice'][0]['rangestep'])
    r = np.arange(0, stoprange, rangestep)*1000
    
    data_ = raw['volume']['scan']['slice'][i]['slicedata']['rawdata']['data']
    datadepth = float(raw['volume']['scan']['slice'][i]['slicedata']['rawdata']['@depth'])
    datamin = float(raw['volume']['scan']['slice'][i]['slicedata']['rawdata']['@min'])
    datamax = float(raw['volume']['scan']['slice'][i]['slicedata']['rawdata']['@max'])
    data_ = datamin + data_ * (datamax - datamin) / 2 ** datadepth

    if flag==2:
        data_ = data_/3
        for jj in range(int(np.floor(np.size(data_[:,1])/3))):
            data_[jj,:] = data_[3*jj,:] + data_[3*jj+1,:] + data_[3*jj+2,:]
        data_ = data_[range(int(np.floor(np.size(data_[:,1])/3))),:]
    elif flag==1:
        data_ = data_/2
        for jj in range(int(np.floor(np.size(data_[:,1])/2))):
            data_[jj,:] = data_[2*jj,:] + data_[2*jj+1,:]
        data_ = data_[range(int(np.floor(np.size(data_[:,1])/2))),:]

    # If len(azi) == 447 will generate error in wrl.ipol.interpolate_polar
    # "ValueError: operands could not be broadcast together with shapes"
    if len(azi) == 175:
        azi = azi[:-1]
        data_ = data_[:-1,:]

    delta = len(r) - len(np.transpose(data_))
    if delta > 0:
        r = r[:-delta]
    data=data_

    # transformasi dari koordinat bola ke koordinat kartesian
    rangeMesh, azimuthMesh =np.meshgrid(r,azi) # meshgrid azimuth dan range
    lonlatalt = wrl.georef.polar.spherical_to_proj(
        rangeMesh, azimuthMesh, elevation, sitecoords
    ) 
    x, y = lonlatalt[:, :, 0], lonlatalt[:, :, 1]
        

    # proses regriding ke data container yang sudah dibuat sebelumnya
    lonMesh, latMesh=np.meshgrid(lonGrid,latGrid)
    gridLatLon = np.vstack((lonMesh.ravel(), latMesh.ravel())).transpose()
    xy=np.concatenate([x.ravel()[:,None],y.ravel()[:,None]], axis=1)
    radius=r[np.size(r)-1]
    center=[x.mean(),y.mean()]
    gridded = wrl.comp.togrid(
        xy, gridLatLon,
        radius, center, data.ravel(),
        wrl.ipol.Linear
    )
    griddedData = np.ma.masked_invalid(gridded).reshape((len(lonGrid), len(latGrid)))
    dataContainer=np.dstack((dataContainer,griddedData))

Selesai! seluruh proses ekstraksi data sudah selesai. Apakah kepala anda mulai berputar hebat? Hahaha tidak mengapa, dulu pada awalnya admin juga begitu. Memang pabrikan radar memiliki cara tersendiri untuk menyimpan datanya, oleh karena itu cara mengekstraknya pun juga berbeda-beda tingkat kesulitannya. Anda tidak perlu paham detail setiap line pada script, karena admin pun ga langsung paham, itu juga cara-caranya melihat dari sumber yang admin cantumkan dibawah.

Yang penting kita sudah tahu bahwa seluruh data sudah tersimpan dalam variabel dataContainer. Mengapa banyak sekali try-except pada potongan script diatas? Karena memang menyesuaikan dengan versi radar software, jadi admin jaga-jaga aja, so script ini bisa digunakan untuk ekstrak data radar Gematronik yang lama maupun yang baru. Selanjutnya visualisasi, misalkan produk CMAX (Column Maximum).

# plotting CMAX data menggunakan Basemap
cmaxData=np.nanmax(dataContainer[:,:,:],axis=2)
cmaxData[cmaxData<0]=np.nan;cmaxData[cmaxData>100]=np.nan

m=Basemap(llcrnrlat=latMin,urcrnrlat=latMax,\
            llcrnrlon=lonMin,urcrnrlon=lonMax,\
            resolution='i')
x0,y0=m(radarLon,radarLat)
x1,y1=m(lonMesh,latMesh)
clevsZ = [5,10,15,20,25,30,35,40,45,50,55,60,65,70]
colors=['#07FEF6','#0096FF','#0002FE','#01FE03','#00C703','#009902','#FFFE00','#FFC801','#FF7707','#FB0103','#C90002','#980001','#FF00FF','#9800FE']
plt.figure(figsize=(10,10))
m.plot(x0, y0, 'ko', markersize=3)
m.contourf(x1, y1, cmaxData,clevsZ,colors=colors)
m.colorbar(ticks=clevsZ,location='right',size='4%',label='Reflectivity (dBZ)')
m.drawparallels(np.arange(math.floor(latMin),math.ceil(latMax),1.5),labels=[1,0,0,0],linewidth=0.2)
m.drawmeridians(np.arange(math.floor(lonMin),math.ceil(lonMax),1.5),labels=[0,0,0,1],linewidth=0.2)
m.drawcoastlines() 
title='CMAX '+timeEnd.strftime("%Y%m%d-%H%M")+' UTC'
plt.title(title,weight='bold',fontsize=15)
plt.savefig('GEMATest.png',bbox_inches='tight',dpi=200,pad_inches=0.1)
plt.close()

Tada! Inilah hasil running dari full script tersebut :

Seluruh data radar GEMATRONIK dapat diekstrak dan di tampilkan dengan cara tersebut. Full script terdapat pada link ini.

Apakah saat mengesekusi script kalian mendapatkan error seperti ini?

Berarti kalian melewatkan satu langkah penting dalam instalasi modul wradlib. Kalian harus membaca postingan ini Tutorial Instalasi Modul Wradlib di bagian akhir. Selamat mencoba!

Baca juga :

1. Ekstrak Data Radar BARON (Single Pol)
2. Ekstrak Data Radar VAISLA (Single Pol)
3. Ekstrak Data Radar EEC CFRadial NetCDF (Single Pol)

Sumber : https://docs.wradlib.org/en/stable/notebooks/fileio/wradlib_radar_formats.html#Gematronik-Rainbow