Analisis Deforestasi dengan R Menggunakan Hansen Dataset dan gfcanalysis Package

Data Global Forest Watch dari Hansen et al (2013) ¹ merupakan data yang sangat populer dan sangat sering dipakai dalam berbagai analisis, terutama dalam bidang kehutanan dan ekologi.

Dataset Hansen (biasanya disebut demikian) terdiri atas beberapa data yaitu:

• Tree canopy cover for year 2000 (treecover2000)
• Global forest cover gain 2000–2016 (gain)
• Year of gross forest cover loss event (lossyear)
• Data mask (datamask)
• Circa year 2000 Landsat 7 cloud-free image composite (first)
• Circa year 2017 Landsat cloud-free image composite (last)

Kita dapat mengakses data hansen website Global Forest Change.

Untuk mengakses data Hansen melalui R, kita dapat memanfaatkan package gfcanalysis.

Sebelum melakukan analisis, kita perlu memahami definisi yang digunakan.

Pertama, deforestasi dalam analisis ini adalah hilangnya tutupan hutan. Sedangkan hutan, adalah piksel dengan tutupan pohon lebih dari nilai threshold yang kita gunakan.

Dalam latihan kali ini, kita akan gunakan threshold 90, jadi nilai piksel sama dengan atau lebih dari 90 (berarti memiliki 90% tutupan pohon) kita definisikan sebagai hutan.

Apa yang akan kita lakukan

Kita akan melakukan analisis kehilangan tutupan pohon dari tahun 2000-2017 untuk daerah Kota dan Kabupaten Jambi.

Analisis dilakukan menggunakan data Hansen, sehingga untuk lebih jelas memahami analisis yang akan kita lakukan, sebaiknya baca dulu data yang terdapat di dataset Hansen.

Secara teknis, kita akan melakukan langkah-langkah berikut

1. Instal dan load package yang dibutuhkan
2. Download data batas administrasi Provinsi Jambi (atau bisa juga gunakan area yang lain)
3. Download data Hansen
4. Penentuan threshold
5. Perhitungan luas hilangnya tutupan pohon
6. Membuat animasi sederhana

Apa yang kita butuhkan

Dalam tutorial ini, kita membutuhkan:

1. Laptop terinstal R dan RStudio
2. gfcanalysis package
3. Koneksi internet

Langkah 1: Persiapan

Instal dan load package yang dibutuhkan

Kita membutuhkan beberapa package, yaitu gfcanalysis, raster, sp dan rgdal.

Jalankan skrip berikut untuk mengecek apakah package sudah terinstal.

Jika package sudah terinstal, maka akan dilanjutkan untuk load library.

Jika belum, maka RStudio akan instal package terlebih dahulu, baru kemudian me-load package-nya.

if (!require(gfcanalysis)) install.packages('gfcanalysis')
if (!require(rgdal)) install.packages('rgdal')
if (!require(raster)) install.packages('raster')
if (!require(sp)) install.packages('sp')

Tentukan folder output

Untuk menggunakan folder yang sama dengan file R yang digunakan, tulis skrip ini:

output_folder <- "."

Jika menginginkan hasil analisis disimpan di folder terpisah, gunakan skrip di bawah ini. Ganti direktori folder sesuai keinginan.

#pastikan gunakan slash miring kanan.
output_folder <- "C:/R for blog" 

Langkah 2: siapkan area of interest (aoi)

Kita akan menggunakan fungsi pada package raster. Cara untuk melakukan download pernah dibahas pada tulisan ini.

Pertama kita download batas administrasi tingkat kota/ kabupaten.

# download GADM data
idn1 <- getData('GADM', country='IDN', level=1,
                 download = TRUE) #level provinsi
#plot data beserta judulnya
plot(idn1, main = "Provinsi di Indonesia") 

Lalu gunakan skrip ini untuk melihat data atributnya. Perhatikan bahwa nama kabupaten/ kota ada di kolom NAME_1. Cari dan catat nama kabupaten/ kota yang dikehendaki. Pastikan penulisannya sama.

# lihat data atribut
idn1@data

Untuk latihan ini, kita gunakan Provinsi Jambi.

# memilih polygon Kota Jambi dan menyimpanya 
# sebagai objek aoi
aoi <- idn1[idn1$NAME_1 == "Jambi",] 
                                          
plot(aoi, main = "Provinsi Jambi") #plot

Langkah 3: Download data Hansen

Untuk mengetahui scene atau tile mana yang akan kita download, kita dapat menghitung dulu berapa jumlah tile yang dibutuhkan. Kita dapat plot datanya untuk melihat lokasi relatif aoi kita terhadap tile GFC.

# menghitung jumlah tile yang dibutuhkan
tiles <- calc_gfc_tiles(aoi)
print(length(tiles))

# plot aoi dan tile
plot(tiles)
plot(aoi, add=TRUE, lty=2, col="#00ff0050")

Daerah kita cukup sempit sehingga hanya dibutuhkan 1 tile saja.

Kita lanjutkan ke proses download data Hansen.

Sebelum, download, kembali kita harus mengingat bahwa dataset ini terdiri atas beberapa data dan kita tidak membutuhkan semuanya. Secara default, proses download tidak akan download data first dan last. Jika ingin download datanya, silakan tuliskan first_and_last=TRUE

download_tiles(tiles, output_folder)

Hasilnya:

## 1 tiles to download/check.
## 0 file(s) succeeded, 4 file(s) skipped, 0 file(s) failed.

Proses download akan dimulai. Waktu yang dibutuhkan proses ini bervariasi tergantung kecepatan internet.

Jika proses download selesai, lihat hasil download di folder sesuai yang ditunjuk oleh output_folder (lihat kembali skrip di bagian atas jika bingung).

Kurang lebih hasil downloadnya seperti ini.

Langkah 4: Penentuan threshold

Logikanya seperti ini.

Data Hansen membutuhkan threshold kerapatan kanopi pohon (dalam persen) untuk menentukan mana kenampakan hutan mana yang bukan.

Misal kita isikan thresholdnya 90, maka piksel dengan kerapatan kanopi pohon di atas 90 akan dimasukkan sebagai kelas hutan, dan nilai dibawah 90 akan masuk ke kelas bukan hutan.

Dengan demikian, kita baru bisa menghitung luas hutan (dalam hal ini tree cover) yang hilang (loss) dan yang tumbuh (gain).

Penentuan threshold ini akan bervariasi tergantung wilayah kajian. Untuk keperluan latihan ini, kita gunakan saja nilai 90.

threshold <- 90

Pertama kita extract data dari folder, lalu menyimpannya ke format .tif

gfc_extract <- extract_gfc(aoi, output_folder, 
                           filename="Jambi_GFC_extract.tif", 
                           overwrite = TRUE)

Selanjutnya kita thresholding (nilai 90) dan menyimpannya ke format .tif.

gfc_thresholded <- threshold_gfc(gfc_extract, forest_threshold = threshold, 
                                 filename="Jambi_GFC_extract_thresholded.tif",
                                 overwrite = TRUE)
gfc_thresholded

Hasilnya:

## class      : RasterBrick 
## dimensions : 8152, 13476, 109856352, 5  (nrow, ncol, ncell, nlayers)
## resolution : 0.00025, 0.00025  (x, y)
## extent     : 101.1285, 104.4975, -2.77125, -0.73325  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/R for blog/Jambi_GFC_extract_thresholded.tif 
## names      : forest2000, lossyear, gain, lossgain, datamask 
## min values :          0,        0,    0,        0,        1 
## max values :          1,       17,    1,        1,        2

Membaca hasil thresholding

Data hasil thresholding memiliki 5 layer (band), dengan rincian sebagai berikut.

plot(gfc_thresholded)

BAND 1 (forest2000)

Band 1 memuat informasi kelas Hutan (1) dan Bukan Hutan (0). Penentuan kelas Hutan dan Bukan hutan berdasarkan threshold yang sudah kita masukkan di proses sebelumnya (=> 90 –> Hutan, <90 Bukan hutan).

Sekali lagi ingat, tergantung nilai thresholdnya.

BAND 2 (lossyear)

Nilai piksel pada band ini menunjukkan tahun ke berapa hutan dalam suatu piksel hilang. Nilai 0 berarti tidak ada hutan yang hilang.

Nilai 1 menunjukkan hutan hilang di tahun pertama (2001), nilai 2 menunjukkan hutan hilang di tahun kedua (2002), begitu seterusnya sampai nilai 17 (hutan hilang tahun 2017).

BAND 3 (gain)

Band ini menunjukkan piksel-piksel Bukan hutan pada tahun 2000 berubah menjadi hutan di tahun-tahun berikutnya (tidak dirinci pada tahun ke berapa). Nilai 0 menunjukkan tidak ada gain, nilai 1 menunjukkan gain.

BAND 4 (loss and gain)

Band ini menunjukkan piksel-piksel yang mengalami loss dan gain. Tidak dirinci apakah gain dulu baru loss atau sebaliknya. Susah diinterpretasi. Nilai 0 menunjukkan tidak ada loss dan gain, nilai 1 menunjukkan ada loss dan gain.

BAND 5

Band 5 merupakan band untuk ekstrak area daratan dan air. Nilai 0 = No data, Nilai 1 = daratan, nilai 2 = air.

Langkah 5: Perhitungan luas hilangnya tutupan pohon

Perhitungan dilakukan dengan menggunakan fungsi gfc_stats.

gfc_stats <- gfc_stats(aoi, gfc_thresholded)
gfc_stats

## $loss_table
##    year   aoi   cover     loss
## 1  2000 AOI 1 2009452       NA
## 2  2001 AOI 1 1987050 22401.69
## 3  2002 AOI 1 1966125 20925.54
## 4  2003 AOI 1 1954948 11176.39
## 5  2004 AOI 1 1908312 46636.81
## 6  2005 AOI 1 1871005 37306.39
## 7  2006 AOI 1 1811989 59015.92
## 8  2007 AOI 1 1753731 58258.43
## 9  2008 AOI 1 1687250 66481.14
## 10 2009 AOI 1 1623250 63999.30
## 11 2010 AOI 1 1602997 20253.83
## 12 2011 AOI 1 1576768 26228.59
## 13 2012 AOI 1 1495366 81402.40
## 14 2013 AOI 1 1462664 32701.75
## 15 2014 AOI 1 1428413 34250.92
## 16 2015 AOI 1 1402012 26401.06
## 17 2016 AOI 1 1325502 76509.41
## 
## $gain_table
##      period   aoi     gain lossgain
## 1 2000-2016 AOI 1 327949.7 275413.1

Hasilnya terdiri atas dua tabel.

Tabel yang pertama menunjukkan jumlah piksel kehilangan (loss) per tahun.

Tabel yang kedua menunjukkan jumlah piksel gain dan loss and gain.

Kita simpan hasil perhitungan ke format csv. Hasil penyimpanan dapat dilihat di output_folder

write.csv(gfc_stats$loss_table, 
          file='Jambi_GFC_extract_thresholded_losstable.csv', 
          row.names=FALSE
          )
write.csv(gfc_stats$gain_table, 
          file='Jambi_GFC_extract_thresholded_gaintable.csv', 
          row.names=FALSE
          )

Kita juga bisa membuat grafik sederhana untuk melihat trend perubahan luas hutan di area kajian.

barplot(height=gfc_stats$loss_table$cover, 
        names=gfc_stats$loss_table$year,
        xlab="Tahun", 
        ylab="Luas hutan dalam piksel (30m)", 
        main="Luas hutan di Provinsi Jambi")

Langkah 6: Membuat animasi sederhana

Pertama kita buat dulu layer per tahun. Kita simpan file-nya karena file ini akan bermanfaat untuk digunakan kembali di kemudian hari.

gfc_annual_stack <- annual_stack(gfc_thresholded)
writeRaster(gfc_annual_stack, filename="Jambi_GFC_extract_thresholded_annual.tif")
plot(gfc_annual_stack)

Keterangan untuk file tersebut:

Layer-layer tersebut dapat dikombinasikan untuk membuat animasi sederhana.

aoi$label <- "Jambi" # Label the polygon on the plot
Jambi_animate_annual(aoi, gfc_annual_stack, out_dir='.', site_name='Jambi')

Hasil animasinya akan tersimpan dalam file .html.

Animasi hasil analisis dapat dilihat di tautan ini!

1. Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 342, (15 November): 850–853. Data available on-line from: http://earthenginepartners.appspot.com/science-2013-global-forest↩