Validasi MODIS dan WRF

---

Dokumentasi ini dibuat bagaimana saya mengolah data satelit MODIS untuk memvalidasi luaran simulasi WRF.
Produk MODIS yang saya gunakan adalah MODIS atmosphere dan MODIS Land Surface Temperature. Saya menggunakan parameter Retrieved Temperature Profile dan Retrieved Dew Point Temperature Profile pada MODIS atmosphere. Kedua parameter tersebut dibutuhkan untuk menghitung kelembapan relatif dengan persamaan Clausius-Clapeyron.

\[e_s(T) = 6.11*exp[\frac{2.453*10^6}{461}*(\frac{1}{273}-\frac{1}{T})] \] \[e_s(T_d) = 6.11*exp[\frac{2.453*10^6}{461}*(\frac{1}{273}-\frac{1}{T_d})]\] \[ RH = \frac{e_s(T_d)}{e_s(T)}*100{\%} \]

Pengolahan data MODIS

Setelah mencari informasi dari berbagai sumber seperti github, stackexchange, stackoverflow, saya mengolah data MODIS di R dengan bantuan package rgdal, sp, gdalUtils, gstat, dan raster. Kumpulan sub data yang berada di dalam file hdf diambil menggunakan gdalUtils dengan fungsi get_subdatasets(), gdal_translate(), dan gdal_warp()

1. Buka library

#Spasial dan Raster
library(raster); library(sp); library(rgdal); library(gdalUtils); library(lattice); library(latticeExtra)

## Loading required package: sp

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.4-0

library(ncdf4) ; library(PCICt)

#Manipulasi dan Visualisasi Data
library(dplyr); library(ggplot2); library(plotly)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:latticeExtra':
## 
##     layer

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:raster':
## 
##     select

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

2. Buka file MODIS

#Untuk ekstrak variabel suhu udara dan suhu titik embun
mod07.path    <- list.files(path = "hdf_MOD07", pattern = "*.hdf", full.names = TRUE)
myd07.path    <- list.files(path = "hdf_MYD07", pattern = "*.hdf", full.names = TRUE)

#Untuk ekstrak variabel suhu permukaan -> MOD11A1 dan MYD11A1
skt.mod.files <- list.files("hdf_MOD11A1", "*.hdf", full.names = T)
skt.myd.files <- list.files("hdf_MYD11A1", "*.hdf", full.names = T)

Beberapa parameter di dalam file hdf memiliki banyak variabel yang memiliki 3D (longitude, latitude, dan pressure) atau 2D (longitude dan latitude). Variasi ketinggian pada tekanan 5 hingga 1000 mb, tetapi yang akan saya gunakan hanya pada ketinggian 50 hingga 1000 mb untuk memvalidasi model WRF.

3. Mengatur penamaan file

Langkah ini digunakan untuk menentukan nama file hasil konversi dari hdf ke tif

#atur nama MOD07 dan MYD07
nama.mod07 <- strsplit(mod07.path,"\\.")
nama.myd07 <- strsplit(myd07.path,"\\.")

#atur nama MOD11A1 dan MYD11A1
skt.mod.nama  <- strsplit(skt.mod.files,"\\.") 
skt.myd.nama  <- strsplit(skt.myd.files,"\\.")

4. Buat file TIF dari HDF

Pada langkah ini, beberapa sub dataset di dalam HDF akan diekstrak menjadi GeoTIFF. Kemudian, resolusi MOD07 dari 5 km menjadi 1 km. Langkah ini menggunakan library gdalUtils dengan fungsi get_subdatasets() dan gdalwarp().

#---Konversi HDF ke RASTER untuk MOD07 dan MYD07---#
for (mod in 1:length(mod07.path)){ #Untuk MOD07
  sds.temp <- get_subdatasets(mod07.path[mod])[15]   #->  Air Temperature profile
  sds.mixr <- get_subdatasets(mod07.path[mod])[16]   #->  Dew Temperature profile
  
  gdalwarp(sds.temp, paste0("tif_MOD07_Temp_Profile/",paste("MOD",nama.mod07[[mod]][2],nama.mod07[[mod]][3],
                                           nama.mod07[[mod]][4],nama.mod07[[mod]][5],sep = "_"),".tif"),
           tps = T,t_srs = "+proj=longlat +datum=WGS84 +ellps=WGS84", output_Raster = T, overwrite = T,
           r = "near",    #resampling method
           tr= c(0.02,0.02)) #resolution (0.02 from 0.05 degree)
  
  gdalwarp(sds.mixr, paste0("tif_MOD07_Tdew_Profile/",paste("MOD",nama.mod07[[mod]][2],nama.mod07[[mod]][3],
                                           nama.mod07[[mod]][4],nama.mod07[[mod]][5],sep = "_"),".tif"),
           tps = T,t_srs = "+proj=longlat +datum=WGS84 +ellps=WGS84", output_Raster = T, overwrite = T,
           r = "near",    #resampling method
           tr= c(0.02,0.02)) #resolution (0.02 from 0.05 degree)
}
for (myd in 1:length(myd07.path)){ #Untuk MYD07
  sds.temp <- get_subdatasets(myd07.path[myd])[15]   #->  Air Temperature profile
  sds.mixr <- get_subdatasets(myd07.path[myd])[16]   #->  Dew Temperature profile
  
  gdalwarp(sds.temp, paste0("tif_MYD07_Temp_Profile/",paste("MYD",nama.myd07[[myd]][2],nama.myd07[[myd]][3],
                                           nama.myd07[[myd]][4],nama.myd07[[myd]][5],sep = "_"),".tif"),
           tps = T,t_srs = "+proj=longlat +datum=WGS84 +ellps=WGS84", output_Raster = T, overwrite = T,
           r = "near",    #resampling method
           tr= c(0.02,0.02)) #resolution (0.02 from 0.05 degree)
  
  gdalwarp(sds.mixr, paste0("tif_MYD07_Tdew_Profile/",paste("MYD",nama.myd07[[myd]][2],nama.myd07[[myd]][3],
                                           nama.myd07[[myd]][4],nama.myd07[[myd]][5],sep = "_"),".tif"),
           tps = T,t_srs = "+proj=longlat +datum=WGS84 +ellps=WGS84", output_Raster = T, overwrite = T,
           r = "near",    #resampling method
           tr= c(0.02,0.02)) #resolution (0.02 from 0.05 degree)
}

#---Konversi HDF ke RASTER untuk MOD11A1 dan MYD11A1---#
for (mod in 1:length(skt.mod.files)){ #Untuk MOD11A1
  skt.day    <- get_subdatasets(skt.mod.files[mod])[1]  #->  Suhu permukaan siang hari
  skt.night  <- get_subdatasets(skt.mod.files[mod])[5]  #->  Suhu permukaan malam hari
  time.day   <- get_subdatasets(skt.mod.files[mod])[3]  #->  Waktu siang hari
  time.night <- get_subdatasets(skt.mod.files[mod])[7]  #->  Waktu malam hari
  
  gdal_translate(skt.day, paste0("hdf_MOD11A1/LST_day/",paste("MOD",skt.mod.nama[[mod]][2],sep = "_"),".tif"), 
                 output_Raster = T)
  
  gdal_translate(skt.night, paste0("hdf_MOD11A1/LST_night/",paste("MOD",skt.mod.nama[[mod]][2],sep = "_"),".tif"), 
                 output_Raster = T)
  
  gdal_translate(time.day, paste0("hdf_MOD11A1/solar_time_day/",paste("MOD",skt.mod.nama[[mod]][2],sep = "_"),".tif"), 
                 output_Raster = T)
  
  gdal_translate(time.night, paste0("hdf_MOD11A1/solar_time_night/",paste("MOD",skt.mod.nama[[mod]][2],sep = "_"),".tif"),
                 output_Raster = T)  
}
for (myd in 1:length(skt.myd.files)){ #Untuk MYD11A1
  skt.day    <- get_subdatasets(skt.myd.files[myd])[1]  #->  Suhu permukaan siang hari
  skt.night  <- get_subdatasets(skt.myd.files[myd])[5]  #->  Suhu permukaan malam hari
  time.day   <- get_subdatasets(skt.myd.files[myd])[3]  #->  Waktu siang hari
  time.night <- get_subdatasets(skt.myd.files[myd])[7]  #->  Waktu malam hari
  
  gdal_translate(skt.day, paste0("hdf_MYD11A1/LST_day/",paste("MYD",skt.myd.nama[[myd]][2],sep = "_"),".tif"), 
           output_Raster = T)
  
  gdal_translate(skt.night, paste0("hdf_MYD11A1/LST_night/",paste("MYD",skt.myd.nama[[myd]][2],sep = "_"),".tif"), 
           output_Raster = T)
  
  gdal_translate(time.day, paste0("hdf_MYD11A1/solar_time_day/",paste("MYD",skt.myd.nama[[myd]][2],sep = "_"),".tif"), 
           output_Raster = T)
  
  gdal_translate(time.night, paste0("hdf_MYD11A1/solar_time_night/",paste("MYD",skt.myd.nama[[myd]][2],sep = "_"),".tif"), 
           output_Raster = T)  
}

5. Cek lokasi wilayah kajian

File perlu dicek lokasi kajian terlebih dahulu karena tidak semua digunakan untuk proses interpolasi. Gunung Kelud berada di antara tiga kabupaten, yaitu Kediri, Blitar, dan Malang. Data shapefile ketiga kabupaten tersebut digunakan untuk mengecek dan ekstraksi nilai raster MODIS yang sudah di resample

#Masukkan data shapefile Gunung Kelud dan Jawa Timur
kelud = shapefile("shp_kelud/daerah_gunung_kelud.shp")
kelud = kelud[kelud$WADMKK %in% c("Blitar","Malang","Kediri"), ]   # --> Memilih tiga kabupaten
jatim = shapefile("../Jatim/jatim.shp")

#Buka File raster MOD07 dan MYD07
##File suhu udara
mod07.temp.path = list.files(path = "tif_MOD07_Temp_Profile", pattern = "*.tif", full.names = TRUE)
myd07.temp.path = list.files(path = "tif_MYD07_Temp_Profile", pattern = "*.tif", full.names = TRUE)

##File suhu titik embun
mod07.tdew.path = list.files(path = "tif_MOD07_Tdew_Profile", pattern = "*.tif", full.names = TRUE)
myd07.tdew.path = list.files(path = "tif_MYD07_Tdew_Profile", pattern = "*.tif", full.names = TRUE)

#Buka File raster MOD11A1 dan MYD11A1
skt.mod.raster.day     <- list.files("hdf_MOD11A1/LST_day"         , "*.tif", full.names = F)
skt.mod.raster.night   <- list.files("hdf_MOD11A1/LST_night"       , "*.tif", full.names = F)
time.mod.raster.day    <- list.files("hdf_MOD11A1/solar_time_day"  , "*.tif", full.names = F)
time.mod.raster.night  <- list.files("hdf_MOD11A1/solar_time_night", "*.tif", full.names = F)

skt.myd.raster.day     <- list.files("hdf_MYD11A1/LST_day"         , "*.tif", full.names = F)
skt.myd.raster.night   <- list.files("hdf_MYD11A1/LST_night"       , "*.tif", full.names = F)
time.myd.raster.day    <- list.files("hdf_MYD11A1/solar_time_day"  , "*.tif", full.names = F)
time.myd.raster.night  <- list.files("hdf_MYD11A1/solar_time_night", "*.tif", full.names = F)

cek.koordinat.mod07.temp = c(); cek.koordinat.mod07.tdew = c()
cek.koordinat.myd07.temp = c(); cek.koordinat.myd07.tdew = c()

# ------Cek lokasi data satelit yang merekam sekitar gunung kelud------- #
# 1. Suhu udara
for (nfile in 1:length(mod07.temp.path)){
  ekstrak <- raster(mod07.temp.path[nfile], band = 17)
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){cek.koordinat.mod07.temp[nfile] <- FALSE}
  else (cek.koordinat.mod07.temp[nfile] <- TRUE)
}

for (nfile in 1:length(myd07.temp.path)){
  ekstrak <- raster(myd07.temp.path[nfile], band = 17)
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){cek.koordinat.myd07.temp[nfile] <- FALSE}
  else (cek.koordinat.myd07.temp[nfile] <- TRUE)
}

# 2. Suhu titik embun
for (nfile in 1:length(mod07.tdew.path)){
  ekstrak <- raster(mod07.tdew.path[nfile], band = 17)
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){cek.koordinat.mod07.tdew[nfile] <- FALSE}
  else (cek.koordinat.mod07.tdew[nfile] <- TRUE)
}

for (nfile in 1:length(myd07.tdew.path)){
  ekstrak <- raster(myd07.tdew.path[nfile], band = 17)
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){cek.koordinat.myd07.tdew[nfile] <- FALSE}
  else (cek.koordinat.myd07.tdew[nfile] <- TRUE)
}

simpulan.mod07.temp = cbind(mod07.temp.path,cek.koordinat.mod07.temp)
simpulan.myd07.temp = cbind(myd07.temp.path,cek.koordinat.myd07.temp)

simpulan.mod07.tdew = cbind(mod07.tdew.path,cek.koordinat.mod07.tdew)
simpulan.myd07.tdew = cbind(myd07.tdew.path,cek.koordinat.myd07.tdew)


#FILE YANG AKAN DIGUNAKAN UNTUK INTERPOLASI
mod07.pakai.temp = na.omit(ifelse(simpulan.mod07.temp[,2] == TRUE, simpulan.mod07.temp[,1], NA))
myd07.pakai.temp = na.omit(ifelse(simpulan.myd07.temp[,2] == TRUE, simpulan.myd07.temp[,1], NA))

mod07.pakai.tdew = na.omit(ifelse(simpulan.mod07.tdew[,2] == TRUE, simpulan.mod07.tdew[,1], NA))
myd07.pakai.tdew = na.omit(ifelse(simpulan.myd07.tdew[,2] == TRUE, simpulan.myd07.tdew[,1], NA))

#Yang akan digunakan untuk interpolasi
mod.pakai.temp = c(mod07.pakai.temp, myd07.pakai.temp)               
mod.pakai.tdew = c(mod07.pakai.tdew, myd07.pakai.tdew)

#Mengatur nama file
mod.pisah.nama.temp = unlist(strsplit(mod.pakai.temp, "/"))
mod.pisah.nama.temp = mod.pisah.nama.temp[seq(2,length(mod.pisah.nama.temp),2)]
mod.pisah.nama.temp = unlist(strsplit(mod.pisah.nama.temp,"_"))

mod.pisah.nama.tdew = unlist(strsplit(mod.pakai.tdew, "/"))
mod.pisah.nama.tdew = mod.pisah.nama.tdew[seq(2,length(mod.pisah.nama.tdew),2)]
mod.pisah.nama.tdew = unlist(strsplit(mod.pisah.nama.tdew,"_"))

tanggal.temp     = mod.pisah.nama.temp[seq(2,length(mod.pisah.nama.temp),5)]
jam.temp         = mod.pisah.nama.temp[seq(3,length(mod.pisah.nama.temp),5)]
mod.or.myd.temp  = mod.pisah.nama.temp[seq(1,length(mod.pisah.nama.temp),5)]

tanggal.tdew     = mod.pisah.nama.tdew[seq(2,length(mod.pisah.nama.tdew),5)]
jam.tdew         = mod.pisah.nama.tdew[seq(3,length(mod.pisah.nama.tdew),5)]
mod.or.myd.tdew  = mod.pisah.nama.tdew[seq(1,length(mod.pisah.nama.tdew),5)]

# 3. Suhu Permukaan

6. Interpolasi

#Interpolasi untuk raster suhu udara ----
for (nfile in 1:length(mod.pakai.temp)){
  for (ketinggian in 3:18){
    ekstrak           <- raster(mod.pakai.temp[nfile], band = ketinggian)
    ekstrak[]         <- 0.009999999776482582*(ekstrak[] + 15000) - 273
    
    ekstrak.df        <- data.frame(rasterToPoints(ekstrak))
    names(ekstrak.df) <- c("lon","lat","suhu")
    
    write.csv(ekstrak.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(
      gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("tif_Temp_P_IDW/", 
                                                               paste(mod.or.myd.temp[nfile], tanggal.temp[nfile],
                                                                     jam.temp[nfile], ketinggian, sep = "_"),".tif"),
                a="linear:nodata=NA",
                of="GTiff",output_Raster=TRUE))
  }
}

#Interpolasi untuk raster suhu titik embun ----
for (nfile in 1:length(mod.pakai.tdew)){
  for (ketinggian in 5:18){
    ekstrak           <- raster(mod.pakai.tdew[nfile], band = ketinggian)
    ekstrak[]         <- 0.009999999776482582*(ekstrak[] + 15000) - 273
    
    ekstrak.df        <- data.frame(rasterToPoints(ekstrak))
    names(ekstrak.df) <- c("lon","lat","suhu")
    
    write.csv(ekstrak.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(
      gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("tif_Tdew_P_IDW/", 
                                                               paste(mod.or.myd.tdew[nfile], tanggal.tdew[nfile],
                                                                     jam.tdew[nfile], ketinggian, sep = "_"),".tif"),
                a="linear:nodata=NA",
                of="GTiff",output_Raster=TRUE))
  }
}
#Interpolasi untuk raster suhu permukaan ----
##Masukkan file raster
skt.mod.raster.day     <- list.files("hdf_MOD11A1/LST_day"         , "*.tif", full.names = T)
skt.mod.raster.night   <- list.files("hdf_MOD11A1/LST_night"       , "*.tif", full.names = T)
time.mod.raster.day    <- list.files("hdf_MOD11A1/solar_time_day"  , "*.tif", full.names = T)
time.mod.raster.night  <- list.files("hdf_MOD11A1/solar_time_night", "*.tif", full.names = T)

skt.myd.raster.day     <- list.files("hdf_MYD11A1/LST_day"         , "*.tif", full.names = T)
skt.myd.raster.night   <- list.files("hdf_MYD11A1/LST_night"       , "*.tif", full.names = T)
time.myd.raster.day    <- list.files("hdf_MYD11A1/solar_time_day"  , "*.tif", full.names = T)
time.myd.raster.night  <- list.files("hdf_MYD11A1/solar_time_night", "*.tif", full.names = T)

#### MOD11A1
for (nfile in 1:length(skt.mod.raster.day)){
  ekstrak   <- raster(readGDAL(skt.mod.raster.day[nfile], as.is = T))
  ekstrak[] <- 0.02*ekstrak[] - 273
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MOD11A1/ekstrak_kelud_skt_day/", "MOD_",
                                                               skt.mod.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
              output_Raster=TRUE))
    
  }
}

for (nfile in 1:length(skt.mod.raster.night)){
  ekstrak   <- raster(readGDAL(skt.mod.raster.night[nfile], as.is = T))
  ekstrak[] <- 0.02*ekstrak[] - 273
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MOD11A1/ekstrak_kelud_skt_night/", "MOD_",
                                                               skt.mod.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
              output_Raster=TRUE))
  }
}

for (nfile in 1:length(time.mod.raster.day)){
  ekstrak   <- raster(readGDAL(time.mod.raster.day[nfile], as.is = T))
  ekstrak[] <- 0.1*ekstrak[]
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MOD11A1/ekstrak_kelud_solar_time_day/", "MOD_",
                                                               skt.mod.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
              output_Raster=TRUE))
  }
}

for (nfile in 1:length(time.mod.raster.night)){
  ekstrak   <- raster(readGDAL(time.mod.raster.night[nfile], as.is = T))
  ekstrak[] <- 0.1*ekstrak[]
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MOD11A1/ekstrak_kelud_solar_time_night/", "MOD_",
                                                               skt.mod.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
              output_Raster=TRUE))
    }
}

#### MYD11A1
for (nfile in 1:length(skt.myd.raster.day)){
  ekstrak   <- raster(readGDAL(skt.myd.raster.day[nfile], as.is = T))
  ekstrak[] <- 0.02*ekstrak[] - 273
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MYD11A1/ekstrak_kelud_skt_day/", "MYD_",
                                                                         skt.myd.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
                        output_Raster=TRUE))
    
  }
}

for (nfile in 1:length(skt.myd.raster.night)){
  ekstrak   <- raster(readGDAL(skt.myd.raster.night[nfile], as.is = T))
  ekstrak[] <- 0.02*ekstrak[] - 273
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MYD11A1/ekstrak_kelud_skt_night/", "MYD_",
                                                                         skt.myd.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
                        output_Raster=TRUE))
    
  }
}

for (nfile in 1:length(time.myd.raster.day)){
  ekstrak   <- raster(readGDAL(time.myd.raster.day[nfile], as.is = T))
  ekstrak[] <- 0.1*ekstrak[]
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MYD11A1/ekstrak_kelud_solar_time_day/", "MYD_",
                                                                         skt.myd.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
                        output_Raster=TRUE))
  }
}

for (nfile in 1:length(time.myd.raster.night)){
  ekstrak   <- raster(readGDAL(time.myd.raster.night[nfile], as.is = T))
  ekstrak[] <- 0.1*ekstrak[]
  
  #Sesuaikan sistem koordinat shapefile dengan raster
  ekstrak <- projectRaster(ekstrak, crs = crs(jatim))
  
  #Cek error -> jika error berarti tidak termasuk di wilayah kelud -> FALSE
  ekstrak          <- crop(ekstrak,jatim)
  ekstra.df        <- data.frame(rasterToPoints(ekstrak))
  names(ekstra.df) <- c("lon","lat","suhu")
  
  if (dim(ekstra.df)[1] == 0){
    print(nfile)
  } else {
    write.csv(ekstra.df, "interp.csv", row.names = FALSE)
    
    vrt_header <- c(
      '<OGRVRTDataSource>',
      '\t<OGRVRTLayer name="interp">',
      '\t<SrcDataSource>interp.csv</SrcDataSource>',
      '\t<GeometryType>wkbPoint</GeometryType>', 
      '\t<GeometryField encoding="PointFromColumns" x="lon" y="lat" z="suhu"/>',
      '\t</OGRVRTLayer>',
      '\t</OGRVRTDataSource>')
    
    vrt_filecon <- file("interp.vrt","w")
    writeLines(vrt_header,con=vrt_filecon)
    close(vrt_filecon)
    
    #Buat file interpolasi menggunakan metode linear
    setMinMax(gdal_grid(src_datasource="interp.vrt", dst_filename=paste0("hdf_MYD11A1/ekstrak_kelud_solar_time_night/", "MYD_",
                                                                         skt.myd.nama[[nfile]][2],".tif"), a="nearest:nodata=NA",
                        output_Raster=TRUE))
  }
}

7. Perhitungan Kelembapan Relatif

#Masukkan file raster yang telah diinterpolasi
tdew.files <- list.files("tif_Tdew_P_IDW/","*.tif",full.names = T)
tair.files <- list.files("tif_Temp_P_IDW/","*.tif",full.names = T)

#atur penamaan (salah satu parameter saja karena namanya sama)
nama <- unlist(strsplit(tdew.files, "/"))
counter = 3

#Buat file raster RH
for (nfile in 1:length(tdew.files)) {
  ekstrak.1 <- raster(tdew.files[nfile])
  ekstrak.1[] <- 6.11*exp((2.453*10^6)/461*((1/273)-(1/(ekstrak.1[]+273))))

  ekstrak.2 <- raster(tair.files[nfile])
  ekstrak.2[] <- 6.11*exp((2.453*10^6)/461*((1/273)-(1/(ekstrak.2[]+273))))
  
  extent(ekstrak.1) <- extent(ekstrak.2)

  rh <- ekstrak.1/ekstrak.2
  
  writeRaster(rh, filename = paste0("tif_RH_P/",nama[counter]), overwrite=T)
  
  counter = counter + 3
}

8. Cek ulang raster dengan shapefile Kelud

#Cek data NA yang masuk daerah gunung kelud ----
#Masukkan file raster RH
rhum.files <- list.files("tif_RH_P/","*.tif",full.names = T)

# 1. Suhu udara
counter = 3
for (nfile in 1:length(tair.files)){
  ekstrak <- raster(tair.files[nfile])
  
  #Cek error 
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){print(nfile)}
  else {
    ekstrak <- crop(ekstrak, kelud)
    writeRaster(ekstrak, filename = paste0("cek_Tair_kelud/",nama[counter]), overwrite=T)
  }
  counter = counter + 3
}

# 2. Suhu titik embun
counter = 3
for (nfile in 1:length(tdew.files)){
  ekstrak <- raster(tdew.files[nfile])
  
  #Cek error
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){print(nfile)}
  else {
    ekstrak <- crop(ekstrak, kelud)
    writeRaster(ekstrak, filename = paste0("cek_Tdew_kelud/",nama[counter]), overwrite=T)
  }
  counter = counter + 3
}

# 3. Kelembapan relatif
counter = 3
for (nfile in 1:length(rhum.files)){
  ekstrak <- raster(rhum.files[nfile])
  
  #Cek error
  ek.tr   <- try(crop(ekstrak,kelud), silent = T)
  
  if (inherits(ek.tr, 'try-error')){print(nfile)}
  else {
    ekstrak <- crop(ekstrak, kelud)
    writeRaster(ekstrak, filename = paste0("cek_rh_kelud/",nama[counter]), overwrite=T)
  }
  counter = counter + 3
}

Validasi data MODIS dan WRF

Data hasil simulasi WRF terlebih dahulu diolah dengan NCL untuk mendapatkan variabel suhu udara dan kelembapan relatif. Waktu WRF disesuaikan dengan data MODIS. Parameter suhu udara dan kelembapan relatif perlu diinterpolasikan terlebih dahulu untuk menyesuaikan ketinggian antara MODIS dan WRF. Rata-rata ketinggian di G. Kelud sekitar 850 mb. Oleh karena itu, ketinggian yang diambil dari MODIS maupun WRF antara 850 - 70 mb.

1. Buka file WRF hasil interpolasi ketinggian dan MODIS

#Buka file WRF
ta.wrf   <- ncvar_get(nc_open("temp_interp_for_modis.nc"))
rh.wrf   <- ncvar_get(nc_open("rh_interp_for_modis.nc"))
skt.wrf     <- ncvar_get(nc_open("../WRF_chem/wrfout_d03_2014-01-10_16.nc"), "TSK")[14,22,] - 273

#Buka file MODIS
ta.raster <- list.files("cek_Tair_kelud/","*.tif",full.names = T)
rh.raster <- list.files("cek_rh_kelud/"  ,"*.tif",full.names = T)

#Mengatur urutan nama file MODIS 
ta.rearrange.raster = list()
rh.rearrange.raster = list()
counter_1        = 1     ; counter_2     = 12

for (nfile in 1:485){
  ta.rearrange.raster[[nfile]] <- ta.raster[counter_1:counter_2]
  rh.rearrange.raster[[nfile]] <- rh.raster[counter_1:counter_2]
  
  t.x.raster                  <- ta.rearrange.raster[[nfile]][-c(1:8)]
  t.y.raster                  <- c(t.x.raster, ta.rearrange.raster[[nfile]][1:8])
  
  r.x.raster                  <- rh.rearrange.raster[[nfile]][-c(1:8)]
  r.y.raster                  <- c(r.x.raster, rh.rearrange.raster[[nfile]][1:8])
  
  ta.rearrange.raster[[nfile]] <- rev(t.y.raster)
  rh.rearrange.raster[[nfile]] <- rev(r.y.raster)
  
  counter_1 = counter_1 + 12
  counter_2 = counter_2 + 12
}

ta.raster = unlist(ta.rearrange.raster)   #Setelah diurutkan berdasarkan tekanan
rh.raster = unlist(rh.rearrange.raster)

2. Ambil Variabel Waktu dan Ketinggian dari WRF

waktu <- ncvar_get(nc_open("temp_interp_for_modis.nc"), "time")
waktu <- as.PCICt(waktu*60*60, "proleptic_gregorian", origin = "2014-01-10") + 16*3600
waktu <- waktu + 7*3600  #Sudah waktu lokal -> WIB

lev   <- ncvar_get(nc_open("temp_interp_for_modis.nc"), "lev")

3. Pengolahan Data

# Masukkan file WRF
ta.wrf.df  <- data.frame(date.time =  as.POSIXct(rep(as.character(waktu), each = 16)), `suhu udara` = c(ta.wrf),
                           ketinggian = rep(lev, times = 1000), stringsAsFactors = F, check.names = F)
rh.wrf.df  <- data.frame(date.time = as.POSIXct(rep(as.character(waktu), each = 16)), `kelembapan relatif` = c(rh.wrf), 
                           ketinggian = rep(lev, times = 1000), stringsAsFactors = F, check.names = F)
skt.wrf.df <- data.frame(date.time = as.POSIXct(as.character(waktu)), `suhu permukaan` = skt.wrf)

#Hilangkan data NA khusus untuk suhu udara dan kelembapan relatif -> NA berarti yang ketinggiannya 50, 920, 950, dan 1000 mb
ta.wrf.df <- na.omit(ta.wrf.df)
rh.wrf.df <- na.omit(rh.wrf.df)

#Masukkan informasi tanggal dan ketinggian dari file MODIS
ta.raster.df <- data.frame(ta.files = c(ta.raster), ketinggian = ta.wrf.df$ketinggian[1:12], stringsAsFactors = F)
rh.raster.df <- data.frame(rh.files = c(rh.raster), ketinggian = rh.wrf.df$ketinggian[1:12], stringsAsFactors = F)

for (nfile in 1:length(ta.raster)) {
  ta.raster.df$jam[nfile]      <- strsplit(ta.raster[nfile],"[[:punct:]]|A")[[1]][8]
  ta.raster.df$yyyyddd[nfile]  <- strsplit(ta.raster[nfile],"[[:punct:]]|A")[[1]][7]
  
  ta.raster.df$yyyy[nfile]     <- as.numeric(stringr::str_c(strsplit(ta.raster.df$yyyyddd[nfile],"")[[1]][1:4], collapse = ""))
  ta.raster.df$ddd[nfile]      <- as.numeric(stringr::str_c(strsplit(ta.raster.df$yyyyddd[nfile],"")[[1]][5:7], collapse = ""))
}

ta.raster.df$jam  <- round(as.numeric(ta.raster.df$jam)/100)

ta.raster.df$tanggal <- format(as.Date(ta.raster.df$ddd, origin = "2012-12-31"), "%d")
ta.raster.df$bulan   <- format(as.Date(ta.raster.df$ddd, origin = "2012-12-31"), "%m")

ta.raster.df$date <- paste(ta.raster.df$yyyy, ta.raster.df$bulan, ta.raster.df$tanggal, sep = "-")
ta.raster.df$date <- paste(ta.raster.df$date, ta.raster.df$jam)
ta.raster.df$date.time <- as.POSIXct(ta.raster.df$date, format = "%Y-%m-%d %H")
ta.raster.df      <- ta.raster.df[,-c(3:9)]  #remove column

# Lakukan pada parameter RH
rh.raster.df$date.time <- ta.raster.df$date.time

# Ekstrak data titik dari raster (ta.raster.df dan rh.raster.df)
for (nfile in 1:length(ta.raster)){
  ekstrak <- raster(ta.raster.df$ta.files[nfile])
  
  koordinat <- sp::SpatialPoints(expand.grid(lon = xFromCol(ekstrak, seq(ekstrak@ncols)), 
                                             lat = yFromRow(ekstrak, seq(ekstrak@nrows))))
  crs(koordinat) <- crs(ekstrak)
  kel.koor <- as.data.frame(koordinat@coords)
  kel.koor <- filter(kel.koor ,lon >= 112 & lon <= 112.5, lat <= -7.8 & lat >= -8.2)
  
  #Cek error
  if (inherits(extract(ekstrak, kel.koor), 'try-error') | dim(kel.koor)[1] == 0){
    ta.raster.df$`suhu udara`[nfile] <- NA
  } else {
    ta.raster.df$`suhu udara`[nfile] <- median(extract(ekstrak, kel.koor),na.rm = T)
  }
}
for (nfile in 1:length(rh.raster)){
  ekstrak <- raster(rh.raster.df$rh.files[nfile])
  
  koordinat <- sp::SpatialPoints(expand.grid(lon = xFromCol(ekstrak, seq(ekstrak@ncols)), 
                                             lat = yFromRow(ekstrak, seq(ekstrak@nrows))))
  crs(koordinat) <- crs(ekstrak)
  kel.koor <- as.data.frame(koordinat@coords)
  kel.koor <- filter(kel.koor ,lon >= 112 & lon <= 112.5, lat <= -7.8 & lat >= -8.2)
  
  if (inherits(median(extract(ekstrak, kel.koor),na.rm = T), 'try-error') | dim(kel.koor)[1] == 0){
    rh.raster.df$`kelembapan relatif`[nfile] <- NA
  } else {
    rh.raster.df$`kelembapan relatif`[nfile] <- median(extract(ekstrak, kel.koor),na.rm = T)
  }
}
rh.raster.df$`kelembapan relatif` <- rh.raster.df$`kelembapan relatif`*100 #konversi ke persentase 

# Ekstrak data titik dari suhu permukaan (ta.raster.df dan rh.raster.df)
skt.mod.crop.day       <- list.files("hdf_MOD11A1/ekstrak_kelud_skt_day"         , "*.tif")
skt.mod.crop.night     <- list.files("hdf_MOD11A1/ekstrak_kelud_skt_night"       , "*.tif")
time.mod.crop.day      <- list.files("hdf_MOD11A1/ekstrak_kelud_solar_time_day"  , "*.tif")
time.mod.crop.night    <- list.files("hdf_MOD11A1/ekstrak_kelud_solar_time_night", "*.tif")

skt.myd.crop.day       <- list.files("hdf_MYD11A1/ekstrak_kelud_skt_day"         , "*.tif")
skt.myd.crop.night     <- list.files("hdf_MYD11A1/ekstrak_kelud_skt_night"       , "*.tif")
time.myd.crop.day      <- list.files("hdf_MYD11A1/ekstrak_kelud_solar_time_day"  , "*.tif")
time.myd.crop.night    <- list.files("hdf_MYD11A1/ekstrak_kelud_solar_time_night", "*.tif")

### list files yg telah dicrop
skt.mod.crop.day   <- data.frame(`MOD11A1` = c(skt.mod.crop.day  , time.mod.crop.day), 
                                 `Keterangan` = c(rep("raster",97), rep("time",96)), stringsAsFactors = F)

skt.myd.crop.day   <- data.frame(`MYD11A1` = c(skt.myd.crop.day  , time.myd.crop.day), 
                                 `Keterangan` = c(rep("raster",91), rep("time",91)), stringsAsFactors = F)

skt.mod.crop.night <- data.frame(`MOD11A1` = c(skt.mod.crop.night, time.mod.crop.night), 
                                 `Keterangan` = c(rep("raster",38), rep("time",38)), stringsAsFactors = F)

skt.myd.crop.night <- data.frame(`MYD11A1` = c(skt.myd.crop.night, time.myd.crop.night), 
                                 `Keterangan` = c(rep("raster",50), rep("time",50)), stringsAsFactors = F)
###

### list files sebelum di crop
skt.mod.day   <- data.frame(`MOD11A1` = c(skt.mod.raster.day, time.mod.raster.day), 
                            `Keterangan` = c(rep("raster",114), rep("time",114)), stringsAsFactors = F)
skt.myd.day   <- data.frame(`MYD11A1` = c(skt.myd.raster.day, time.myd.raster.day), 
                            `Keterangan` = c(rep("raster",113), rep("time",113)), stringsAsFactors = F)

skt.mod.night <- data.frame(`MOD11A1` = c(skt.mod.raster.night, time.mod.raster.night), 
                            `Keterangan` = c(rep("raster",114), rep("time",114)), stringsAsFactors = F)
skt.myd.night <- data.frame(`MYD11A1` = c(skt.myd.raster.night, time.myd.raster.night), 
                            `Keterangan` = c(rep("raster",113), rep("time",113)), stringsAsFactors = F)

## Merge
skt.mod.pakai.day <- intersect(skt.mod.crop.day, skt.mod.day, by = "MOD11A1")
skt.myd.pakai.day <- intersect(skt.myd.crop.day, skt.myd.day, by = "MYD11A1")

skt.mod.pakai.night <- intersect(skt.mod.crop.night, skt.mod.night, by = "MOD11A1")
skt.myd.pakai.night <- intersect(skt.myd.crop.night, skt.myd.night, by = "MYD11A1")

## Ekstrak nilai
skt.mod.crop.day       <- list.files("hdf_MOD11A1/ekstrak_kelud_skt_day"         , "*.tif", full.names = T)
skt.mod.crop.night     <- list.files("hdf_MOD11A1/ekstrak_kelud_skt_night"       , "*.tif", full.names = T)
time.mod.crop.day      <- list.files("hdf_MOD11A1/ekstrak_kelud_solar_time_day"  , "*.tif", full.names = T)
time.mod.crop.night    <- list.files("hdf_MOD11A1/ekstrak_kelud_solar_time_night", "*.tif", full.names = T)

skt.myd.crop.day       <- list.files("hdf_MYD11A1/ekstrak_kelud_skt_day"         , "*.tif", full.names = T)
skt.myd.crop.night     <- list.files("hdf_MYD11A1/ekstrak_kelud_skt_night"       , "*.tif", full.names = T)
time.myd.crop.day      <- list.files("hdf_MYD11A1/ekstrak_kelud_solar_time_day"  , "*.tif", full.names = T)
time.myd.crop.night    <- list.files("hdf_MYD11A1/ekstrak_kelud_solar_time_night", "*.tif", full.names = T)

for (nfile in 1:dim(skt.mod.pakai.day)[1]) {
  ekstrak <- raster(c(skt.mod.crop.day, time.mod.crop.day)[nfile])
  if (inherits(raster(c(skt.mod.crop.day, time.mod.crop.day)[nfile]), 'try-error')){
    skt.mod.pakai.day$value[nfile] <- NA
  } else {
    
    if (inherits(extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8))), 'try-error')){
      skt.mod.pakai.day$value[nfile] <- NA
      } else {
        ekstrak <- extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8)))
        skt.mod.pakai.day$value[nfile] <- quantile(ekstrak, probs = 0.95, na.rm = T)
        }
  }
}
for (nfile in 1:dim(skt.mod.pakai.night)[1]) {
  ekstrak <- raster(c(skt.mod.crop.night, time.mod.crop.night)[nfile])
  if (inherits(raster(c(skt.mod.crop.night, time.mod.crop.night)[nfile]), 'try-error')){
    skt.mod.pakai.night$value[nfile] <- NA
  } else {
    
    if (inherits(extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8))), 'try-error')){
      skt.mod.pakai.night$value[nfile] <- NA
      } else {
        ekstrak <- extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8)))
        skt.mod.pakai.night$value[nfile] <- quantile(ekstrak, probs = 0.95, na.rm = T)
        }
  }
}
for (nfile in 1:dim(skt.myd.pakai.day)[1]) {
  ekstrak <- raster(c(skt.myd.crop.day, time.myd.crop.day)[nfile])
  if (inherits(raster(c(skt.myd.crop.day, time.myd.crop.day)[nfile]), 'try-error')){
    skt.myd.pakai.day$value[nfile] <- NA
  } else {
    
    if (inherits(extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8))), 'try-error')){
      skt.myd.pakai.day$value[nfile] <- NA
      } else {
        ekstrak <- extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8)))
        skt.myd.pakai.day$value[nfile] <- quantile(ekstrak, probs = 0.95, na.rm = T)
        }
  }
}
for (nfile in 1:dim(skt.myd.pakai.night)[1]) {
  ekstrak <- raster(c(skt.myd.crop.night, time.myd.crop.night)[nfile])
  if (inherits(raster(c(skt.myd.crop.night, time.myd.crop.night)[nfile]), 'try-error')){
    skt.myd.pakai.night$value[nfile] <- NA
  } else {
    
    if (inherits(extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8))), 'try-error')){
      skt.myd.pakai.night$value[nfile] <- NA
      } else {
        ekstrak <- extract(ekstrak, cbind(c(112.28,112.28,112.341,112.341),c(-7.913,-8,-7.913,-8)))
        skt.myd.pakai.night$value[nfile] <- quantile(ekstrak, probs = 0.95, na.rm = T)
        }
  }
}

#Penyesuaian data MODIS dengan WRF rh.wrf.df <=> rh.raster.df dan ta.wrf.df <=> ta.raster.df
ta.wrf.mod <- left_join(ta.raster.df, ta.wrf.df, by = c("date.time","ketinggian"), suffix = c(".MODIS", ".WRF"))
ta.wrf.mod <- na.omit(ta.wrf.mod)
rh.wrf.mod <- left_join(rh.raster.df, rh.wrf.df, by = c("date.time","ketinggian"), suffix = c(".MODIS", ".WRF"))
rh.wrf.mod <- na.omit(rh.wrf.mod)

# Import data csv suhu permukaan ----
skt.modis <- read.csv("skt_modis.csv", stringsAsFactors = F)
skt.modis$date.time <- as.POSIXct(skt.modis$date.time, format = "%m/%d/%Y %H:%M")
skt.modis$date.time <- as.POSIXct(format.POSIXct(skt.modis$date.time, "%Y-%m-%d %H:%M:%S"))

# Merge MODIS dari WRF
skt.merge <- left_join(skt.modis, skt.wrf.df, by = "date.time", suffix = c(".MODIS",".WRF"))
skt.merge <- na.omit(skt.merge)

4. Plotting

load("variabelData.RData")

#Suhu Udara 850mb
ggplot(data = filter(ta.wrf.mod, ketinggian == 850), mapping = aes(x = date.time)) + 
  geom_line(aes(y = `suhu udara.MODIS`, col = "MODIS")) + 
  geom_line(aes(y = `suhu udara.WRF`, col = "WRF")) + 
  xlab("Waktu") + ylab(expression("Suhu Udara 850 mb" ~ "("^0~"C)")) + theme(legend.title = element_blank())

#Kelembapan relatif 850mb
ggplot(data = filter(rh.wrf.mod, ketinggian == 850), mapping = aes(x = date.time)) + 
  geom_line(aes(y = `kelembapan relatif.MODIS`, col = "MODIS")) + 
  geom_line(aes(y = `kelembapan relatif.WRF`, col = "WRF")) + 
  xlab("Waktu") + ylab("Kelembapan Relatif 850 mb (%)") + theme(legend.title = element_blank())

#Kontur Profil Suhu Udara
fig.ta.modis <- plot_ly(data = ta.wrf.mod, x = ~date.time, y = ~ketinggian, z = ~`suhu udara.MODIS`, 
                      type = "contour", contours = list(start = 30, end = -80, size = 10, showlabels = TRUE), 
                      colorscale = list(colorRampPalette(c("red","yellow","blue")))) %>%
  layout(xaxis = list(title = "Waktu (UTC)"), yaxis = list(title = "Ketinggian (mb)", autorange="reversed")) %>%
  colorbar(title = "Suhu Udara 850 mb \n (derajat C)" )

fig.ta.wrf <- plot_ly(data = ta.wrf.mod, x = ~date.time, y = ~ketinggian, z = ~`suhu udara.WRF`,
                      type = "contour", contours = list(start = 30, end = -80, size = 10, showlabels = TRUE), 
                      colorscale = list(colorRampPalette(c("red","yellow","blue"))), showlegend = FALSE) %>%
  layout(xaxis = list(title = "Waktu (UTC)"), yaxis = list(title = "Ketinggian (mb)", autorange="reversed")) %>% 
  hide_colorbar()

subplot(fig.ta.modis, fig.ta.wrf, shareY = T, titleX = T)

#Kontur Profil kelembapan relatif
fig.rh.modis <- plot_ly(data = rh.wrf.mod, x = ~date.time, y = ~ketinggian, z = ~`kelembapan relatif.MODIS`, 
                      type = "contour", contours = list(start = 0, end = 100, size = 25, showlabels = TRUE), 
                      colorscale = list(colorRampPalette(c("red","blue")))) %>%
  layout(xaxis = list(title = "Waktu (UTC)"), yaxis = list(title = "Ketinggian (mb)", autorange="reversed")) %>% 
  colorbar(title = "Kelembapan Relatif (%)")

fig.rh.wrf <- plot_ly(data = rh.wrf.mod, x = ~date.time, y = ~ketinggian, z = ~`kelembapan relatif.WRF`, 
                      type = "contour", contours = list(start = 0, end = 100, size = 25, showlabels = TRUE), 
                      colorscale = list(colorRampPalette(c("red","blue")))) %>%
  layout(xaxis = list(title = "Waktu (UTC)"), yaxis = list(title = "Ketinggian (mb)", autorange="reversed")) %>% 
  hide_colorbar()

subplot(fig.rh.modis, fig.rh.wrf, shareY = T, titleX = T)

# Plot Suhu Permukaan
ggplot(skt.merge, aes(date.time)) + 
  geom_line(aes(y = suhu.permukaan.MODIS, col = "MODIS")) +
  geom_line(aes(y = suhu.permukaan.WRF, col = "WRF")) +
  xlab("Waktu") + ylab("Suhu Permukaan (derajat C)") + theme(legend.title = element_blank())

plot_ly(skt.merge, x = ~date.time) %>%
  add_trace(y = ~suhu.permukaan.MODIS, name = "MODIS", mode = "lines") %>%
  add_trace(y = ~suhu.permukaan.WRF  , name = "WRF"  , mode = "lines") %>%
  layout(xaxis = list(title = "Waktu (WIB)"), 
         yaxis = list(title = "Suhu Permukaan (derajat C)"))

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter