1. Apa itu data NetCDF?

NetCDF (ekstensi .nc) adalah seperangkat pustaka perangkat lunak dan self-describing. NetCDF dikembangkan dan dikelola oleh Unidata. Self-describing artinya ada header yang mendeskripsikan tata letak file, khususnya larik data, serta metadata file arbitrer dalam bentuk atribut nama/nilai. NetCDF mendukung pembuatan, akses, dan berbagi data ilmiah berorientasi larik (array). Larik di sini adalah data yang memiliki \(\geq\) 1 dimensi. Saat ini, netCDF mendukung dalam membaca data HDF5. Bahkan, data netCDF dapat langsung diolah secara online menggunakan Open-source Project for a Network Data Access Protocol (OPeNDAP). Beberapa perangkat lunak dan bahasa pemrograman yang mendukung dalam membaca data netCDF seperti:

  1. NCAR Command Language (NCL)
  2. GrADS
  3. Matrix Laboratory (MATLAB)
  4. ArcMAP
  5. ncBrowse
  6. Panoply
  7. QGIS
  8. R
  9. Python
  10. Julia
Gambar 1. Struktur data netCDF (Sumber: https://simulatingcomplexity.files.wordpress.com/2014/11/netcdf-file-structure.png)

Gambar 1. Struktur data netCDF (Sumber: https://simulatingcomplexity.files.wordpress.com/2014/11/netcdf-file-structure.png)

Data netCDF sering digunakan untuk penyimpanan data-data dalam ilmu kebumian, khususnya dalam bidang meteorologi dan klimatologi. Sebagai contoh data-data luaran model seperti ERA5 dan ERA-Interim, bahkan data observasi seperti CHIRPS menggunakan format netCDF.

2. Package R untuk Pengolahan Data NetCDF

Berdasarkan hasil penelusuran, saya menemukan beberapa package yang digunakan untuk membaca file netCDF di R. Beberapa di antaranya adalah:

  1. ncdf4 (Pierce 2019)
  2. RNetCDF (Michna and Woods 2020)

Sementara itu, beberapa package tambahan pendukung yang direkomendasikan oleh penulis untuk kedua package di atas adalah:

  1. raster (Hijmans 2020)
  2. sp (Bivand, Pebesma, and Gomez-Rubio 2013)
  3. rgdal (Bivand, Keitt, and Rowlingson 2020)
  4. tidync (Sumner 2020)
  5. PCICt (Pacific Climate Impacts Consortium; portions based on code written by the R-Core team and Drepper. 2018)
  6. rasterVis (Perpiñán and Hijmans 2020)

Sebenarnya, tidak ada perbedaan yang terlalu mencolok antara package ncdf4 dan RNetCDF. Hanya saja, saya terbiasa menggunakan package ncdf4 untuk pengolahan data-data nc dan ditambah dengan package raster. Di dalam artikel ini, saya menggunakan package ncdf4 untuk membaca data nc. Selanjutnya, package raster dan rasterVis untuk menampilkan grafik dan memodifikasi data nc. Package nomor 1, 2, 3, 5, dan 6 sering saya gunakan, tetapi untuk nomor 4 belum pernah. Berdasarkan pembuat package, tidync berguna untuk mengetahui isi dari data nc dengan informasi yang lebih mudah dipahami jika dibandingkan hanya menggunakan ncdf4. Anda bisa mengikuti tutorial dari package tidync di laman https://ropensci.org/blog/2019/11/05/tidync/#overview. Saya menggunakan R versi 4.0 pada sistem operasi Ubuntu 20.04.

> R.version.string
[1] "R version 4.0.2 (2020-06-22)"
> osVersion
[1] "Ubuntu 20.04.1 LTS"

3. Tahapan Instalasi NetCDF

Anda hanya perlu mengetik perintah berikut.

> install.packages(c("ncdf4","raster","rasterVis","PCICt"))

atau pada aplikasi RStudio dapat dilakukan dengan memilih menu Tools > Install Packages. Kemudian, muncul kotak dialog seperti di bawah ini.

Gambar 2. Instalasi package pada RStudio

Gambar 2. Instalasi package pada RStudio

Setelah selesai melakukan instalasi, buka package tersebut dengan mengetik perintah berikut.

> library(ncdf4); library(raster); library(rasterVis); library(PCICt)

4. Membuka Data NetCDF

Saya menggunakan data curah hujan permukaan bulanan (1981-2019) dari CHIRPS dengan resolusi \(1^o \times 1^o\) yang diunduh dari halaman web IRI Data Library (IRIDL). Anda bisa langsung mengimpor data tersebut dengan langsung mengetik/copas web hyperlink CHIRPS karena IRIDL mendukung OPeNDAP.

a. Package ncdf4

> chirps.ncdf4 <- nc_open("chirps_monthly.nc")

Anda dapat melihat isi dari data nc tersebut seperti di bawah ini. Cara ini hampir sama dengan perintah ncdump di terminal Bash

> chirps.ncdf4
File chirps_monthly.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float precipitation[X,Y,T]   
            pointwidth: 1
            standard_name: lwe_precipitation_rate
            long_name: CHIRPS precipitation
            expires: 1596240000
            ncolor: 254
            CE: 1500
            CS: 0
            maxncolor: 254
            colormap: [13882323 16777215 16777215 [16777184 5] [15453831 3] [16748574 16] [13434880 22] [4026644 22] [7451452 21] [9419919 43] [9234160 43] [6333684 42] [2237106 43] 2237106]
            colorscalename: prcp_monthlyrate_max1500_smooth
            file_missing_value: -9999
            history: [(Boxes with less than 25.0% dropped) (Boxes with less than 25.0% dropped)]
            missing_value: NaN
            units: mm/month

     3 dimensions:
        X  Size:360
            pointwidth: 1
            gridtype: 1
            units: degree_east
        Y  Size:99
            pointwidth: 1
            gridtype: 0
            units: degree_north
        T  Size:468
            standard_name: time
            pointwidth: 1
            long_name: Time
            expires: 1596240000
            calendar: 360
            gridtype: 0
            units: months since 1960-01-01
> ncdump -h chirps_monthly.nc
netcdf chirps_monthly {
dimensions:
    X = 360 ;
    Y = 99 ;
    T = 468 ;
variables:
    float X(X) ;
        X:pointwidth = 1.f ;
        X:gridtype = 1 ;
        X:units = "degree_east" ;
    float Y(Y) ;
        Y:pointwidth = 1.f ;
        Y:gridtype = 0 ;
        Y:units = "degree_north" ;
    float T(T) ;
        T:standard_name = "time" ;
        T:pointwidth = 1.f ;
        T:long_name = "Time" ;
        T:expires = 1596240000 ;
        T:calendar = "360" ;
        T:gridtype = 0 ;
        T:units = "months since 1960-01-01" ;
    float precipitation(T, Y, X) ;
        precipitation:pointwidth = 1.f ;
        precipitation:standard_name = "lwe_precipitation_rate" ;
        precipitation:long_name = "CHIRPS precipitation" ;
        precipitation:expires = 1596240000 ;
        precipitation:ncolor = 254 ;
        precipitation:CE = 1500.f ;
        precipitation:CS = 0.f ;
        precipitation:maxncolor = 254 ;
        precipitation:colormap = "[13882323 16777215 16777215 [16777184 5] [15453831 3] [16748574 16] [13434880 22] [4026644 22] [7451452 21] [9419919 43] [9234160 43] [6333684 42] [2237106 43] 2237106]" ;
        precipitation:colorscalename = "prcp_monthlyrate_max1500_smooth" ;
        precipitation:file_missing_value = -9999.f ;
        precipitation:history = "[(Boxes with less than 25.0% dropped) (Boxes with less than 25.0% dropped)]" ;
        precipitation:missing_value = NaNf ;
        precipitation:units = "mm/month" ;
}

b. Package raster

Selain menggunakan fungsi nc_open(), fungsi brick() dari package raster juga dapat digunakan untuk membaca data nc. Akan tetapi, package ncdf4 juga diperlukan sebagai prasyarat dalam menjalankan fungsi tersebut.

> chirps.raster <- brick("chirps_monthly.nc", varname="precipitation")
> chirps.raster
class      : RasterBrick 
dimensions : 99, 360, 35640, 468  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : -0.5, 359.5, -49.5, 49.5  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /home/agung/Documents/CDAV/chirps_monthly.nc 
names      : X252.5, X253.5, X254.5, X255.5, X256.5, X257.5, X258.5, X259.5, X260.5, X261.5, X262.5, X263.5, X264.5, X265.5, X266.5, ... 
T (months since 1960-01-01): 252.5, 719.5 (min, max)
varname    : precipitation

c. Package tidync

Saya jarang menggunakan library(tidync) untuk melihat informasi data nc. Mari kita coba! ^_^

> chirps.tidync <- tidync("chirps_monthly.nc")
> chirps.tidync

Data Source (1): chirps_monthly.nc ...

Grids (4) <dimension family> : <associated variables> 

[1]   D0,D1,D2 : precipitation    **ACTIVE GRID** ( 16679520  values per variable)
[2]   D0       : X
[3]   D1       : Y
[4]   D2       : T

Dimensions 3 (all active): 
  
  dim   name  length   min   max start count  dmin  dmax unlim coord_dim 
  <chr> <chr>  <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl> <lgl> <lgl>     
1 D0    X        360    0   359      1   360    0   359  FALSE TRUE      
2 D1    Y         99  -49    49      1    99  -49    49  FALSE TRUE      
3 D2    T        468  252.  720.     1   468  252.  720. FALSE TRUE

Data CHIRPS di atas memiliki nilai curah hujan pada variabel precipitation yang memiliki 3 dimensi, yaitu dimensi X, Y, dan T. Dimensi X, Y, dan T secara berturut-turut adalah longitude, latitude, dan waktu. Artinya, data CHIRPS memiliki nilai masing-masing terhadap longitude (360), latitude (99), dan waktu (468). Panjang waktu data ini sebanyak 468 bulan (39 tahun). Untuk memahami isi dari informasi nc, Anda perlu melihat ilustrasi dari Gambar 1.

Menurut saya, tidync berguna dalam meringkas metadata menjadi data frame sehingga dapat memudahkan pembacaan informasi tersebut.

5. Memilih variabel

Umumnya, data netCDF dapat berisi lebih dari 1 variabel sehingga ukurannya dapat mencapai beberapa GB. Kita dapat memilih variabel yang ingin dianalisis atau divisualisasikan dengan fungsi ncvar_get() pada package ncdf4. Lebih mudahnya, kita dapat menggunakan fungsi brick() seperti pada poin 4b. Perbedaan pada kedua fungsi tersebut adalah tipe dan ukuran. Fungsi ncvar_get() menghasilkan tipe larik sehingga ukurannya lebih besar jika dibandingkan dengan fungsi brick() dengan tipe RasterBrick.

> #Package ncdf4
> var.precip.ncdf4 <- ncvar_get(chirps.ncdf4, varid = "precipitation") 
> #Package raster
> chirps.raster <- brick("chirps_monthly.nc", varname="precipitation")
> crs(chirps.raster) <- crs(raster()) #Menambahkan sistem koordinat referensi
> print(object.size(var.precip.ncdf4), units = "MB")
127.3 Mb
> print(object.size(chirps.raster), units = "MB")
0.1 Mb

Saya lebih suka menggunakan fungsi brick() karena tidak terlalu memakan banyak RAM daripada ncvar_get(). Selain itu, variabel lokasi (longitude dan latitude) dan waktu tidak perlu dipanggil dengan objek baru ketika kita ingin membuat grafik. Lain halnya dengan ncvar_get(), kita perlu mendeklarasikan/memanggil variabel lokasi dan waktu.

> latitude  <- ncvar_get(chirps.ncdf4, varid = "Y")
> longitude <- ncvar_get(chirps.ncdf4, varid = "X")
> waktu     <- ncvar_get(chirps.ncdf4, varid = "T")
> head(latitude) #Hanya menampilkan 6 elemen
[1] -49 -48 -47 -46 -45 -44
> head(longitude)
[1] 0 1 2 3 4 5
> head(waktu)
[1] 252.5 253.5 254.5 255.5 256.5 257.5

Hasil cetakan objek latitude dan longitude di atas menunjukkan data tersebut memiliki resolusi spasial sebesar \(1^o \times 1^o\). Untuk objek waktu, hasil cetakan berupa numerik dan perlu dikonversi menjadi format tanggal. Package PCICt berperan di sini untuk mengubah waktu dalam format numerik menjadi date. Anda perlu melihat informasi data netCDF pada variabel T. Tipe kalender yang digunakan pada data ini adalah 360 yang berarti semua bulan dalam setahun diasumsikan memiliki jumlah tanggal sebanyak 30. Selain tipe kalender 360, ada tipe yang lain seperti gregorian, proleptic_gregorian, noleap, dan 365 (dapat dibaca pada dokumentasi package PCICt atau ketik ?PCICt::as.PCICT()).

Fungsi as.PCICt() memiliki argumen origin=... yang perlu dicantumkan di dalam fungsi ketika kelas objek waktu adalah numerik. Pada unit atribut waktu (lihat hasil cetakan chirps.ncdf4), tercantum months since 1960-01-01 yang berarti data ini adalah data bulanan yang terhitung sejak tanggal 1 Januari 1960. Tanggal tersebut sebagai acuan dari argumen origin=.... Nilai dari objek waktu dikonversi dari bulanan ke detik (1 bulan = 2592000 detik).

> konv.waktu <- as.PCICt(waktu*2592000, cal = "360", origin = "1960-01-01")
> head(konv.waktu)
[1] "1981-01-16" "1981-02-16" "1981-03-16" "1981-04-16" "1981-05-16"
[6] "1981-06-16"

Tips: Jika Anda sudah mengetahui waktu data netCDF, Anda hanya perlu mendeklarasikan sendiri waktu tersebut dengan fungsi seq(). Contoh: seq(as.Date("1981-01-01"), as.Date("2019-12-01"), by = "month"). Dengan fungsi tersebut, data bulanan akan terbentuk secara berurutan terhadap bulan (by = “month”). Anda juga dapat membuat urutan waktu harian dengan hanya mengganti by = "month" menjadi by = "day".

Untuk tutorial ini, saya fokus pada objek yang berbentuk raster.

6. Menambahkan Waktu ke Objek Raster

Objek chirps.raster sebetulnya sudah memiliki waktu, tetapi masih berbentuk numerik. Kita bisa mengubah waktu dengan format numerik menjadi format date dari objek konv.waktu. Anda bisa melihat label time memiliki format tanggal. Hal ini berguna untuk memanipulasi data jika kita ingin mengetahui klimatologi bulanan atau musiman.

> chirps.raster <- setZ(chirps.raster, konv.waktu)
> chirps.raster
class      : RasterBrick 
dimensions : 99, 360, 35640, 468  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : -0.5, 359.5, -49.5, 49.5  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /home/agung/Documents/CDAV/chirps_monthly.nc 
names      : X252.5, X253.5, X254.5, X255.5, X256.5, X257.5, X258.5, X259.5, X260.5, X261.5, X262.5, X263.5, X264.5, X265.5, X266.5, ... 
time       : 1981-01-16, 2019-12-16 (min, max)
varname    : precipitation

7. Memilih Lokasi Spesifik

Anda dapat memilih wilayah kajian tertentu dengan fungsi crop(). Jika Anda memiliki data shapefile batas administratif, Anda bisa langsung menerapkannya pada fungsi tersebut. Sebagai contoh, saya menggunakan batas garis pantai Indonesia yang diunduh dari https://gadm.org/ level 0 dengan format .rds.

> indo.rds <- readRDS("gadm36_IDN_0_sp.rds")
> chirps.raster.indo <- crop(chirps.raster, indo.rds)

8. Memanipulasi Data

Saya memberikan contoh di bawah ini bagaimana cara menjumlahkan data bulanan untuk setiap tahun dari 1981-2019 dan menghitung bulanan klimatologi (1981-2019). Data klimatologi di sini maksudnya adalah merata-ratakan setiap bulan dalam periode tertentu, misalnya Jan 1981-2019, Feb 1981-2019, …, Des 1981-2019 sehingga diperoleh luaran datanya selama 12 bulan. Saya menggunakan objek chirps.raster.indo.

a. Data Bulanan Menjadi Data Tahunan

> chirps.raster.indo.yearly <- stackApply(chirps.raster.indo, rep(seq(1, nlayers(chirps.raster.indo)/12), each = 12), sum)
> names(chirps.raster.indo.yearly) <- seq(1981, 2019) #Memberikan nama tahun

b. Data Bulanan Menjadi Data Klimatologi Bulanan

> chirps.raster.indo.monthly <- stackApply(chirps.raster.indo, rep(seq(1, 12), each = nlayers(chirps.raster.indo)/12), mean)
> names(chirps.raster.indo.monthly) <- month.name #Memberikan nama bulan

9. Membuat Grafik

Fungsi standar yang digunakan untuk membuat grafik dari data netCDF adalah plot(). Secara default, fungsi plot() menghasilkan gambar maksimal sejumlah 16 grafik.

> plot(chirps.raster)

Pada package rasterVis, tampilan grafik spasial menjadi lebih teratur dengan fungsi levelplot() dan dapat ditambahkan garis pantai dengan fungsi layer().

> levelplot(chirps.raster.indo.monthly) + layer(sp.polygons(indo.rds))

Referensi

Bivand, Roger, Tim Keitt, and Barry Rowlingson. 2020. Rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. https://CRAN.R-project.org/package=rgdal.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. https://asdar-book.org/.

Hijmans, Robert J. 2020. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Michna, Pavel, and Milton Woods. 2020. RNetCDF: Interface to ’Netcdf’ Datasets. https://CRAN.R-project.org/package=RNetCDF.

Pacific Climate Impacts Consortium; portions based on code written by the R-Core team, David Bronaugh for the, and Ulrich Drepper. 2018. PCICt: Implementation of Posixct Work-Alike for 365 and 360 Day Calendars. https://CRAN.R-project.org/package=PCICt.

Perpiñán, Oscar, and Robert Hijmans. 2020. rasterVis. http://oscarperpinan.github.io/rastervis/.

Pierce, David. 2019. Ncdf4: Interface to Unidata netCDF (Version 4 or Earlier) Format Data Files. https://CRAN.R-project.org/package=ncdf4.

Sumner, Michael. 2020. Tidync: A Tidy Approach to ’Netcdf’ Data Exploration and Extraction. https://CRAN.R-project.org/package=tidync.