Pada artikel ini, saya akan menunjukkan beberapa tahapan dalam pengolahan data proyeksi perubahan iklim skenario RCP4.5 dan RCP8.5 model MIROC5. Data MIROC5 diunduh dari http://worldclim.org dengan resolusi 1 km.

1. Buka Library

library(raster)
## Loading required package: sp
library(foreach)
library(doParallel)
## Loading required package: iterators
## Loading required package: parallel

2. Masukkan data shapefile wilayah kajian

shp <- shapefile("batas_wilayah.shp")
shp <- spTransform(shp, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ")

3. Memotong data hasil unduhan WorldClim pada model MIROC5 sesuai dengan wilayah kajian

#Masukkan data
#1. Curah Hujan
pr.miroc <- list.files("data proyeksi/pr",".tif",full.names = T)
#2. Suhu Udara Minimum
tn.miroc <- list.files("data proyeksi/tn",".tif",full.names = T)
#3. Suhu Udara Maksimum
tx.miroc <- list.files("data proyeksi/tx",".tif",full.names = T)

#Potong data 
#Melakukan komputasi parallel dengan 4 prosesor. Hal ini dilakukan karena data global dengan resolusi 1 km agar waktu proses lebih cepat 
cl <- makeCluster(4); registerDoParallel(cl)

#Baca list per file raster 
#1. Curah Hujan
pr <- foreach(i=seq_along(pr.miroc)) %dopar% {
  library(raster)
  pr <- raster(pr.miroc[i]); pr <- crop(pr, shp); pr <- mask(pr, shp)
}
pr <- brick(pr)
#2. Suhu udara minimum
tn <- foreach(i=seq_along(tn.miroc)) %dopar% {
  library(raster)
  tn <- raster(tn.miroc[i]); tn <- crop(tn, shp); tn <- mask(tn, shp)
}
tn <- brick(tn)/10
#3. Suhu udara maksimum
tx <- foreach(i=seq_along(tx.miroc)) %dopar% {
  library(raster)
  tx <- raster(tx.miroc[i]); tx <- crop(tx, shp); tx <- mask(tx, shp)
}
tx <- brick(tx)/10

#Stop Cluster
stopCluster(cl)

4. Perhitungan suhu udara rata-rata, curah hujan, dan kelembapan relatif tahunan

Akumulasi curah hujan tahunan.

pr.rcp45.2050 <- calc(pr[[ 1:12]], sum)
pr.rcp45.2070 <- calc(pr[[13:24]], sum)
pr.rcp85.2050 <- calc(pr[[25:36]], sum)
pr.rcp85.2070 <- calc(pr[[37:48]], sum)
pr.mod <- brick(pr.rcp45.2050, pr.rcp45.2070, pr.rcp85.2050, pr.rcp85.2070)

Suhu udara minimum dan maksimum digunakan untuk menghitung suhu udara rata-rata

ta <- (tn + tx)/2
ta.rcp45.2050 <- calc(ta[[ 1:12]], mean)
ta.rcp45.2070 <- calc(ta[[13:24]], mean)
ta.rcp85.2050 <- calc(ta[[25:36]], mean)
ta.rcp85.2070 <- calc(ta[[37:48]], mean)
ta.mod <- brick(ta.rcp45.2050, ta.rcp45.2070, ta.rcp85.2050, ta.rcp85.2070)

Persamaan Clausius-Clapeyron digunakan untuk mengestimasi kelembapan relatif dengan asumsi suhu udara minimum sama dengan suhu titik embun.

#Fungsi RH
rh.fun <- function(tx, tn){
  ea <- 6.1078*exp(17.269*tn/(237.3 + tn))
  es <- 6.1078*exp(17.269*tx/(237.3 + tx))
  return(ea/es*100)
}
rh <- rh.fun(tx, tn)
rh.rcp45.2050 <- calc(rh[[ 1:12]], mean)
rh.rcp45.2070 <- calc(rh[[13:24]], mean)
rh.rcp85.2050 <- calc(rh[[25:36]], mean)
rh.rcp85.2070 <- calc(rh[[37:48]], mean)
rh.mod <- brick(rh.rcp45.2050, rh.rcp45.2070, rh.rcp85.2050, rh.rcp85.2070)

5. Proyeksi Kesesuaian Lahan Tanaman Kopi

Berikut ini adalah penilaian kesesuaian lahan tanaman kopi berdasarkan kriteria Balai Besar Litbang Sumber Daya Pertanian

Parameter Cuaca S1 S2 S3 N
Suhu udara 16-20 15-16 dan 20-22 14-15 dan 22-24 <14 dan >24
Curah hujan 1200-1800 1000-1200 dan 1800-2000 2000-3000 dan 800-1000 >3000 dan <800 
#1. Suhu udara
ta.class <- matrix(c(16, 20, 1,  #S1
                     15, 16, 2,  #S2
                     20, 22, 2,
                     14, 15, 3,  #S3
                     22, 24, 3,
                      0, 14, 4,  #N
                     24,Inf, 4), nrow = 7, ncol = 3, byrow = T)
ta.reclass <- reclassify(ta.mod, ta.class)
names(ta.reclass) <- c("RCP45.2050","RCP45.2070","RCP85.2050","RCP85.2070")
#2. Curah hujan
pr.class <- matrix(c(1200, 1800, 1,
                     1000, 1200, 2,
                     1800, 2000, 2,
                      800, 1000, 3,
                     2000, 3000, 3,
                        0,  800, 4,
                     3000,  Inf, 4), nrow = 7, ncol = 3, byrow = T)
pr.reclass <- reclassify(pr.mod, pr.class)
names(pr.reclass) <- c("RCP45.2050","RCP45.2070","RCP85.2050","RCP85.2070")
## class      : RasterBrick 
## dimensions : 95, 133, 12635, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 96.25833, 97.36667, 4.175, 4.966667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : memory
## names      : RCP45.2050, RCP45.2070, RCP85.2050, RCP85.2070 
## min values :          1,          1,          1,          1 
## max values :          4,          4,          4,          4
## class      : RasterBrick 
## dimensions : 95, 133, 12635, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 96.25833, 97.36667, 4.175, 4.966667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : memory
## names      : RCP45.2050, RCP45.2070, RCP85.2050, RCP85.2070 
## min values :          2,          3,          2,          3 
## max values :          4,          4,          4,          4

6. Plot

#Suhu Udara
plot(ta.reclass, legend = F, xlab = "Longitude") 
legend(x = 0.377, y = 0.9, legend = c("S1","S2","S3","N"), fill = rev(terrain.colors(4)))

#Curah hujan
plot(pr.reclass, legend = F) 
legend(x = 0.377, y = 0.9, legend = c("S1","S2","S3","N"), fill = rev(terrain.colors(4)))