EKSTRAKSI POLA KEKERINGAN PERTANIAN PULAU JAWA MENGGUNAKAN DATA SATELIT NOAA-18 AVHRR
Weling Suseno , Rokhmatuloh1, dan Agus Wibowo
Abstrak
Kekeringan pertanian merupakan salah satu masalah serius yang membutuhkan perhatian khusus karena dampaknya dapat mempengaruhi kondisi ketahanan pangan nasional, terutama kekeringan pertanian di pulau Jawa sebagai lumbung padi nasional dimana lebih dari 50% produksi padi nasional berasal dari pulau ini. Studi tentang kekeringan telah banyak dilakukan diantaranya untuk identifikasi sebaran kekeringan, analisis faktor-faktor penyebab, hingga prediksi kekeringan, yang tercakup dalam mitigasi bencana kekeringan. Salah satu hambatan besar dari proses tersebut adalah pada tahap pemetaan sebaran kekeringan atau penyediaan informasi kekeringan secara spasial yang up-to-date atau real time. Keterbatasan tersebut kini dapat diatasi dengan menggunakan aplikasi penginderaan jauh. Penelitian ini membahas mengenai esktraksi pola pertanian kekeringan Pulau Jawa menggunakan citra satelit NOAA-18 AVHRR dengan metode triangle. Dari hasil ekstraksi pola kekeringan pulau Jawa menggunakan data satelit NOAA-18 AVHRR, pada bulan Mei 2008 seluruh wilayah pertanian di pulau Jawa berada pada kondisi normal. Pola kekeringan pertanian pulau Jawa bulan September 2008 sebagian besar tersebar di wilayah Jawa Timur dan beberapa Kabupaten di Jawa Tengah.
Kata kunci : Kekeringan pertanian, penginderaan jauh, NOAA-18 AVHRR, metode triangle,
I. Pendahuluan
Kekeringan merupakan sebuah fenomena alam yang biasa terjadi akibat dari pengaruh iklim (White, 1990). Letak Indonesia yang berada di dekat khatulistiwa menjadikan Indonesia memiliki iklim yang panas sehingga rentan terhadap bencana kekeringan. Bencana ini dapat terjadi di berbagai wilayah termasuk wilayah pertanian atau disebut dengan kekeringan pertanian.
Pulau Jawa merupakan salah satu pulau besar di Indonesia dimana sektor pertanian masih menjadi sektor andalan. Hal ini digambarkan dari penggunaan tanah di pulau Jawa yang didominasi oleh pertanian (gambar 1). Pertanian di pulau Jawa memegang peranan yang sangat penting dalam mendukung ketahanan pangan di pulau tersebut dan ketahanan pangan secara nasional karena lebih 50% produksi padi nasional dihasilkan di pulau ini. Sebagai wilayah yang memiliki iklim panas, pertanian di pulau Jawa pun rentan terhadap kekeringan. Dari data yang dicatat oleh Dinas Pertanian Jawa Barat dan Jawa Tengah (2008), hingga Juli 2008 terdapat 20 Kabupaten di Jawa Barat dan 25 Kabupaten di Jawa Tengah yang mengalami kekeringan pertanian diikuti dengan 48.720 ha sawah di Jawa Barat dan 6.870 ha sawah di Jawa Tengah yang mengalami gagal panen. Peran pulau Jawa yang sangat sentral dalam mendukung ketahanan pangan nasional sehingga kejadian kekeringan pada wilayah ini perlu ditanggapi secara serius dan dibutuhkan studi-studi tentang kekeringan dalam mendukung mitigasi bencana kekeringan.
Studi tentang kekeringan telah banyak dilakukan diantaranya untuk identifikasi sebaran kekeringan, analisis faktor-faktor penyebab, hingga prediksi kekeringan, yang tercakup dalam mitigasi bencana kekeringan. Salah satu hambatan besar dari proses tersebut adalah pada tahap pemetaan sebaran kekeringan atau penyediaan informasi kekeringan secara spasial yang up-to-date atau real time. Keterbatasan tersebut kini dapat diatasi dengan menggunakan aplikasi penginderaan jauh.
Salah satu data penginderaan jauh yang banyak digunakan dalam penelitian kekeringan yaitu citra satelit NOAA (National Oceanic and Atmospheric Administration). NOAA merupakan satelit orbit polar yang diciptakan untuk kepentingan meteorologi. Satelit ini membawa sensor AVHRR (Advanced Very High Resolution Radiometer) yang memiliki 5 (lima) saluran, terdiri dari saluran panjang gelombang tampak dan saluran panjang gelombang infra merah termal. Kemampuan sensor AVHRR dalam mendeteksi termal sangat baik digunakan dalam menentukan suhu permukaan tanah yang dapat dijadikan sebagai salah satu indikator yang menggambarkan kekeringan tanah. Disamping itu, satelit NOAA memiliki resolusi temporal hingga dua kali sehari yang berguna dalam mendukung observasi meteorologis yang memiliki tingkat kedinamisan temporal yang tinggi. Hal ini sangat sesuai dengan kebutuhan dalam pemetaan sebaran kekeringan yang memang sangat dipengaruhi oleh keberadaan iklim.
Penelitian ini membahas mengenai esktraksi pola pertanian kekeringan Pulau Jawa menggunakan citra satelit NOAA-18 AVHRR. Identifikasi kekeringan yang dilakukan yaitu pada bulan Mei yang mewakili awal musim kemarau dan bulan September 2008 yang mewakili akhir musim kemarau.
Sumber : Badan Pertanahan Nasional (2000)
Gambar 1. Peta Wilayah Pertanian Pulau Jawa
II. Tujuan Penelitian
Tujuan penelitian ini adalah mengekstrak pola sebaran kekeringan pertanian Pulau Jawa menggunakan citra NOAA-18 AVHRR. Ekstraksi pola kekeringan pertanian yang dilakukan menggunakan metode triangle.
III. Metode Triangle
Metode triangle atau segitiga merupakan metode yang dikembangkan oleh Carlson, et al (1995) dalam mengidentifikasi kekeringan. Metode ini menjelaskan hubungan antara suhu permukaan (LST) dengan indeks vegetasi (NDVI). Segitiga yang dimaksud pada metode ini merupakan hasil bentukan scatterplot dari kedua parameter segitiga tersebut (gambar 2).
Gambar 2. Scatterplot suhu permukaan tanah (LST) dan NDVI (Claps (2004))
Tepi-tepi pada scatterplot menggambarkan jangkauan (range) variasi suhu permukaan tanah. Pada batas atas (tepi panas) menunjukkan kondisi tanah tanpa kelembaban atau kering, sedangkan pada bagian bawah (tepi basah) menunjukkan tanah sangat basah. Dari kedua parameter segitiga tersebut dapat ditentukan indeks kelembaban tanah dari ratio antara jarak sebuah titik pada scatterplot tersebut ke tepi basah (A, dalam gambar 3) dengan jarak antar tepi-tepi (B) yang disebut dengan TVDI (Temperature Vegetation Dryness Index) (Sandholt, 2002). TVDI merupakan indeks kekeringan, dimana nilai maksimum TVDI=1 menunjukkan piksel kering dan piksel basah ditunjukkan dengan nilai TVDI mendekati 0.
(1)
Gambar 3. TVDI dalam segitiga LST dan NDVI (Claps (2004))
Beberapa aspek penting pada metode triangle yaitu scatter plot pada gambar 2, menurun ke kanan (menuju suhu yang lebih rendah) pada saat mengalami peningkatan fraksi vegetasi. Hal ini menjelaskan bahwa vegetasi yang disinari matahari umumnya lebih dingin dibandingkan dengan tanah kosong atau gundul. Segitiga ini menunjukkan sangat kecil variasi dari suhu permukaan pada vegetasi kerapatan tinggi. Tepi panas dan dingin, menunjukkan batas suhu permukaan tanah (suhu tertinggi dan terendah) terhadap fraksi tutupan vegetasi (NDVI), dapat dinyatakan bahwa suhu vegetasi tidak terlalu berubah dan variasi suhu dalam gambaran segitiga hanya pada permukaan tanah. Jadi tepi dingin dan panas menggambarkan piksel terbasah dan terkering.
IV. Pengolahan Citra NOAA-18 AVHRR
• Pra Pengolahan
Pra pengolahan citra dilakukan untuk mendapatkan tampilan citra berada dalam kondisi idealnya, tereferensi geometris terhadap sistem proyeksi koordinat bumi, dan pemilihan wilayah penelitian yang diinginkan pada citra sehingga dapat dilakukan proses pengolahan untuk mendapatkan data yang diinginkan. Adapun tahap-tahap yang dilakukan antara lain :
a. Koreksi Geometrik
Data hasil perekaman citra satelit perlu diregistrasikan kedalam suatu sistem proyeksi atau koordinat bumi. Hal ini perlu dilakukan supaya informasi data citra sesuai keberadaannya di bumi. Pada proses koreksi geometrik ini terdapat dua tahapan. Tahap pertama adalah relokasi piksel ke posisi seharusnya dan proses resampling nilai piksel (Sobirin dkk, 2007).
Dalam melakukan koreksi geometrik hal pertama yang harus dilakukan adalah menentukan titik kontrol. Titik kontrol ini berupa objek yang terlihat pada citra sekaligus terlihat pada peta rujukan yang akan digunakan dalam koreksi geometrik. Penentuan letak titik kontrol berupa lokasi yang identik seperti pada persilangan antara sungai dan jalan ataupun persimpangan jalan, garis pantai, dan beberapa objek lain yang tampak dengan jelas di citra maupun peta rujukan. Peta rujukan yang digunakan dalam melakukan koreksi geometrik yaitu Peta Rupa Bumi Indonesia skala 1:250.000.
b. Subset
Subset merupakan proses pemotongan citra untuk mendapatkan wilayah penelitian yang dinginkan pada citra. Untuk melakukan pemotongan pada citra diperlukan region of interest sebagai pemotong citra, dalam penelitian ini adalah Pulau Jawa.
c. Cloud Mask
Salah satu kekurangan dari data penginderaan jauh multispektral adalah ketidakmampuannya dalam menembus awan, sehingga objek yang ada dipermukaan bumi tertutup awan, tidak dapat terlihat atau mengalami anomali pada respon albedo. Cloud mask pada citra NOAA-AVHRR dilakukan untuk mendeteksi wilayah yang tertutup awan, dengan menggunakan rumus (Saunders dan Kribel, 1988) :
(2)
Jika nilai piksel hasil rasio (r) band 2 dan band 1 lebih besar dari 0,6, maka piksel tersebut menunjukkan awan.
d. Klasifikasi Wilayah Pertanian
Untuk mengklasifikasi wilayah pertanian di Pulau Jawa pada citra NOAA-AVHRR, dilakukan dengan cara memotong citra menggunakan region of interest lahan pertanian di Pulau Jawa yang diperoleh dari peta penggunaan tanah Pulau Jawa tahun 1999, 2000, dan 2001 Badan Pertanahan Nasional.
e. Kalibrasi saluran 1 dan 2 citra NOAA-AVHRR
Kalibrasi data NOAA-AVHRR saluran 1 dan 2 dilakukan berdasarkan referensi dari NOAA user’s guide, dengan menggunakan rumus :
A = S x C10 + I (3)
dimana
A = Albedo (persen)
S = Slope (persen albedo/count)
I = Intercept (persen albedo)
C10 = nilai piksel data AVHRR (10 bit)
f. Kalibrasi saluran 4 dan 5 Citra NOAA-AVHRR
Kalibrasi saluran inframerah termal (saluran 4 dan 5) dilakukan berdasarkan referensi dari NOAA user’s guide.
Pada tiap baris scan in-orbit, sensor AVHRR mengamati tiga tipe target yang berbeda. Pertama, memberikan keluaran 10 count ketika mengamati angkasa (cold space), kemudian count tunggal untuk setiap 2048 piksel target permukaan bumi, dan 10 count ketika mengamati target blackbody internal (sebenarnya hanya cermin scan AVHRR yang berotasi). Target angkasa (cold space) dan blackbody internal digunakan untuk kalibrasi AVHRR, karena nilai radiansi dapat secara independen ditentukan untuk tiap target.
Adapun langkah-langkah dalam mengkalibrasi saluran 4 dan 5 citra NOAA-AVHRR, adalah sebagai berikut :
- Menghitung suhu blackbody internal (TBB)
Suhu target blackbody internal diukur dengan 4 PRT. Pada tiap baris scan, data word 18, 19, dan 20 dalam format frame minor HRPT mempunyai 3 nilai dari 4 PRT. PRT yang berbeda disampling pada tiap baris scan; setiap baris scan ke-lima, semua ketiga nilai PRT adalah 0 yang menunjukkan bahwa satu set 4 data PRT telah disampling. Nilai count CPRT dari tiap PRT dihitung menjadi suhu dengan formula :
TPRT = d0 + d1 CPRT + d2 C2PRT + d3 C3PRT + d4 C4PRT (4)
Nilai koefisien d0, d1, d2, d3, dan d4 untuk tiap PRT ditampilkan pada tabel 1 untuk satelit NOAA-18.
Untuk menghitung suhu blackbody internal TBB, NESDIS menggunakan perata-rataan :
(5)
Tabel 1. Koefisien Konversi NOAA-18 AVHRR/3
PRT d0 d1 d2 d3 d4
1 276,601 0,05090 1,657 E-06 0 0
2 276,683 0,05101 1,482 E-06 0 0
3 276,565 0,05117 1,313 E-06 0 0
4 276,615 0,05103 1,484 E-06 0 0
Sumber : National Climatic Data Center U.S Department of Commerce, 2005
- Menghitung radiansi blackbody internal (NBB)
Radiansi NBB pada tiap saluran termal dari blackbody internal pada suhu TBB adalah rataan terbobot fungsi Planck pada respon spektral saluran tersebut. Fungsi respon spektral untuk tiap saluran diukur pada sekitar 200 internal panjang gelombang dan disediakan bagi NESDIS oleh pembuat instrument. Secara praktis, suatu look-up tabel yang menghubungkan radiansi dengan suhu dibuat untuk tiap saluran. Tiap tabel menunjukkan radiansi pada tiap 1/10 derajat Kelvin antara 180 dan 340 K. Tabel ini disebut “Tabel Energi”. Didapatkan bahwa persamaan dua-langkah berikut secara akurat menghasilkan Tabel Energi setara dengan suhu blackbody dengan ketelitian ± 0,01 K pada range 180 sampai 340K. Tiap saluran termal mempunyai satu persamaan, yang menggunkan bilangan gelombang pusat (centroid wavenumber), vC, dan suhu blackbody efektif, TB. Persamaan dua-langkah tersebut adalah:
T*BB = A + B TBB (6)
(7)
Dimana konstanta radiasi c1 dan c2 adalah :
c1 = 1,1910427 x 10-5 mW/(m2-sr-cm-4)
c2 = 1,4387752 cm-K
Nilai vC dan koefisien A dan B untuk saluran 3B, 4, dan 5 NOAA-18 ditampilkan pada Tabel 2. Bilangan gelombang pusat tunggal untuk tiap saluran menggantikan metode sebelumnya, yang menggunakan bilangan gelombang pusat yang berbeda untuk tiap empat range suhu.
Tabel 2. Koefisien Radiansi NOAA-18 AVHRR/3
vC A B
Saluran 3B 2659,7952 1,698704 0,996960
Saluran 4 928,1460 0,436645 0,998607
Saluran 5 833,2532 0,253179 0,999057
Sumber : National Climatic Data Center U.S Department of Commerce, 2005
- Menghitung radiansi permukaan bumi (NE) menggunakan korelasi non-linear
Keluaran dari dua target kalibrasi in-orbit digunakan untuk menghitung perkiraan linear dari radiansi permukaan bumi NE. Tiap baris scan, AVHRR mengukur target blackbody internal dan mengeluarkan 10 nilai count untuk tiap tiga dtektor saluran termal; yang terletak pada words 23 sampai 52 dalam susunan data HRPT. Ketika AVHRR mengarah ke angkasa (cold space), 10 count dari tiap lima saluran dikeluarkan dan disimpan pada word 53 sampai 102. Nilai count tiap saluran dirata-ratakan untuk menghaluskan noise acak; seringkali counts dari lima baris scan yang berurutan dirata-ratakan karena diperlukan lima baris untuk memperoleh satu set pengukuran seluruh 4 PRT. Count blackbody rata-rata, CBB, dan count angkasa (space) rata-rata, CS, bersama dengan radiansi blackbody NBB dan radiance angkasa, NS, digunakan untuk menghitung perkiraan radiansi linear, NLIN,
(8)
Dimana CE adalah keluaran count AVHRR pada target permukaan bumi (2048 count tiap baris scan).
Detektor saluran termal 3B mempunyai respon linear terhadap radians yang dating sehingga radiansi linear yang dihitung dengan persamaan (8) merupakan nilai sebenarnya untuk saluran 3B. Untuk saluran ini, nilai radiansi angkasa NS adalah 0 sehingga tidak diperlukan koreksi non-linear.
Detektor Mercury-Cadmium-Telluride yang digunakan untuk saluran 4 dan 5 mempunyai respon non-linear terhadap radiansi yang dating. Pengukuran laboratorium pada pre-launch menunjukkan bahwa :
a) Radiansi scene adalah fungsi non-linear (kuadratik) dari count keluaran AVHRR.
b) Ketidaklinearan tersebut tergantung pada suhu operasi AVHRR.
Diasumsikan bahwa respon non-linear akan tetap ada pada saat mengorbit. Untuk seri satelit NOAA KLM (NOAA-15, 16, 17), NESDIS menggunakan metode koreksi non-linear berdasarkan radiansi. Pada metode ini, perkiraan radiansi linear mula-mula dihitung menggunakan radiansi angkasa non-zero, NS pada persamaan (8). Kemudian, nilai radiansi linear dimasukkan ke dalam persamaan kuadrat untuk menghasilkan koreksi radiansi non-linear, NCOR:
NCOR = b0 + b1 NLIN + b2 N2LIN (9)
Akhirnya, radiansi permukaan bumi diperoleh dengan menambahkan NCOR dan NLIN,
NE = NLIN + Ncor (10)
Menetapkan nilai radiansi angkasa non-zero merupakan cara matematis yang mempunyai dua keuntungan utama. Pertama, hanya diperlukan satu persamaan koreksi kuadrati per saluran; koefisien kuadratik adalah tidak bergantung pada suhu operasi AVHRR. Kedua, metode ini menghasilkan pengukuran pre-launch dengan sangat baik; perbedaan RMS antara data fitted dan hasil pengukuran adalah sekitar 0,1 K untuk kedua saluran 4 dan 5. Nilai NS dan koefisien kuadratik b0, b1, dan b2 ditampilkan pada tabel 3 ntuk NOAA-18.
Tabel 3. Koefisien Radiansi Non-linear NOAA-18 AVHRR/3
NS b0 b1 b2
Saluran 4 -5,53 5,82 -0,11069 0,00052337
Saluran 5 -2,22 2,67 -0,04360 0,00017715
Sumber : National Climatic Data Center U.S Department of Commerce, 2005
- Konversi radiansi permukaan bumi (NE) menjadi suhu blackbody (TE)
Suhu TE didefinisikan dengan membuat invers langkah-langkah yang digunakan untuk menghitung radiansi NE yang diukur oleh saluran AVHRR dari blackbody pada suhu TE. Proses dua-langkah tersebut adalah :
(11)
(12)
Nilai vC dan koefisien A dan B ditampilkan pada tabel 2.3.
• Esktraksi Pola Kekeringan Pertanian
Pada pengolahan citra dalam penelitian ini digunakan untuk memperoleh data kejadian kekeringan di Pulau Jawa. Adapun tahap-tahap yang dilakukan antara lain :
a) Menghitung Suhu Permukaan Tanah
Suhu permukaan tanah dapat diukur dengan metode penginderaan jauh menggunakan algoritma split-window (Singh, 1984), dengan rumus :
TS = 1,699 (T4) - 0,699 (T5) – 0,24 (13)
dimana :
TS = Suhu permukaan tanah (°K)
T4 = Saluran 4 citra NOAA-AVHRR yang telah dikalibrasi
T5 = Saluran 5 citra NOAA-AVHRR yang telah dikalibrasi
b) Menghitung Normalized Difference Vegetation Index (NDVI)
Sebelum menentukan nilai Temperature Vegetation Dryness Index (TVDI), terlebih dahulu dilakukan penghitungan nilai Normalized Difference Vegetation Index (NDVI), dengan rumus :
(14)
c) Menghitung Temperature Vegetation Dryness Index (TVDI)
TDVI dapat dihitung menggunakan NDVI dan Ts yang diperoleh dari satelit dengan rumus : (Sandholt, 2002)
(15)
dimana
Ts = Suhu permukaan tanah (°K)
Tsmin = Suhu permukaan tanah minimum (°K)
NDVI = Normalized Difference Vegetation Index
a dan b = Konstanta regresi Ts dan NDVI pada tepi panas segitiga
V. Hasil dan Pembahasan
Dari hasil ekstraksi pola kekeringan pulau Jawa menggunakan data satelit NOAA-18 AVHRR (gambar 4), pada tanggal 3 Mei 2008 seluruh wilayah pertanian di pulau Jawa berada pada kondisi normal. Hal ini diperkirakan bahwa pada masa ini yaitu awal musim kemarau kelembaban tanah pada wilayah pertanian masih tinggi pengaruh dari curah hujan pada bulan sebelumnya.
Pola kekeringan pertanian pulau Jawa tanggal 21 September 2008 sebagian besar tersebar di wilayah Jawa Timur, yaitu Kabupaten Gresik, Mojokerto, Sidoarjo, Pasuruan, Probolinggo, Lumajang, Bondowoso, Situbondo, Jember, Banyuwangi, Ponorogo, Madiun, Magetan, Nganjuk, Kediri dan Bojonegoro. Wilayah Jawa Tengah, yaitu Kabupaten Rembang,Demak, Semarang, Batang, Temanggung, dan Kendal.
Gambar 4. (A) Sebaran Kekeringan Pertanian Pulau Jawa tanggal 3 Mei 2008 hasil ekstraksi citra NOAA-18 AVHRR; (B) Sebaran Kekeringan Pertanian Pulau Jawa tanggal 21 September 2008 hasil ekstraksi citra NOAA-18 AVHRR
VI. Saran
Terdapat beberapa hal yang harus diperhatikan dalam penelitian ini, antara lain :
1. Perlunya observasi lapangan untuk mengetahui tingkat kevalidan hasil ekstraksi pola kekeringan pertanian di pulau Jawa dari data satelit NOAA-18 AVHRR.
2. Perlunya analisa lebih lanjut (analisa spasial) untuk menjelaskan persebaran kekeringan pertanian yang dihasilkan.
Penelitian ini merupakan hasil awal dari penelitian yang dilakukan oleh penulis mengenai perbandingan metode triangle dengan metode termal inersia dalam identifikasi kekeringan pertanian di pulau Jawa menggunakan data satelit NOAA-18 AVHRR.
Daftar Pustaka
Carlson, T.N., Gillies, R.R., Schumugge, T.J. (1995). An Interpretation of Methodologies for Indirect Measurement of Soil Water Content. Agricultural and Forest Meteorology, 77, 191-205.
Claps, P., Laguardia, G. (2004). Assessing Spatial Variability of Soil Water Content Through Thermal Inertia and NDVI. The International Society for Optical Engineering. (28 Agustus 2008). http://www.idrologia.polito.it/~claps/Papers/spie-barcellona2.pdf.
National Climatic Data Center. (2005). NOAA KLM User’s Guide. U.S Department of Commerce. (28 Maret 2006). http://www2.ncdc.noaa.gov/docs.klm. Purwadhi, S. F. (2001). Interpretasi Citra Digital. PT Grasindo, Jakarta.
Saunders, R. W., dan K. T. Kriebel. (1988). An Improved Method for Detecting Clear Sky and Cloudy Radiences from AVHRR Data. Int. Journal of Remote Sensing, vol. 9 no. 1, 123-150.
Sandholt, I., Rasmussen, K., Andersen, J. (2002). A Simple Interpretation of the Surface Temperature-Vegetation Index Space for Assessment of Surface Moisture Status. Remote Sensing of Environment, 79, 213-224.
Singh, S. M. (1984). Removal of Atmospheric effects on a Pixel by Pixel Basis from the Thermal Infrared Data from Instruments on Satellites. The Advanced Very High Resolution Radiometer (AVHRR). Int. Journal of Remote Sensing, 5, 161-183; erratum, ibid., 5, 618.
Sobirin, R. Hermina, D. I. Sari, & Suprayogi. (2007). Modul Praktikum Citra Digital, Menggunakan ER. Mapper 6.4. Departemen Geografi Fakultas Matematika dan Ilmu Pengetahuan Alam Universitas Indonesia, Depok.
White, F.H. (1990). A Study of the Feasibility of Using Simulation Models and Mathematical Programs as Aids to Drought Monitoring and Management. Bureau of Rural Resources, Canbera.
Tidak ada komentar:
Posting Komentar