code_01_processing_MODIS_LST_from_HDF_to_table.R

library(magrittr)
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.3-4, (SVN revision 766)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.3-1

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:magrittr':
## 
##     extract

library(gdalUtils)
library(fst)

setwd("/home/michael/Dropbox/BGU/Itai_Kloog/p_55_processing_MODIS_LST_from_HDF_to_table/")

# Start
start = Sys.time()

# HDF file name
filename = "MOD11A1.A2016183.h18v04.006.2016241045804.hdf"

# Get subdataset paths & variable names
s = get_subdatasets(filename)
vars = gsub(".*:", "", s)

# Read from HDF to multi-band raster
r = stack()
for(i in s) {
  tmp = i %>% readGDAL %>% raster
  r = stack(r, tmp)
}

## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:QC_Day has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Day_view_time has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Day_view_angl has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:LST_Night_1km has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:QC_Night has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Night_view_time has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Night_view_angl has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Emis_31 has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Emis_32 has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Clear_day_cov has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
## HDF4_EOS:EOS_GRID:MOD11A1.A2016183.h18v04.006.2016241045804.hdf:MODIS_Grid_Daily_1km_LST:Clear_night_cov has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns

names(r) = vars

# Examine raster
plot(r)

# To points
pnt = rasterToPoints(r, spatial = TRUE)
pnt = spTransform(pnt, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

# To table
dat = as.data.frame(pnt)

# Rename 'x-y' to 'lon-lat'
dat$lon = dat$x; dat$lat = dat$y
dat$x = NULL; dat$y = NULL

# Add metadata
dat$date = 
  filename %>% 
  strsplit("\\.") %>% 
  sapply("[", 2) %>% 
  gsub("A", "", .) %>% 
  as.Date(format = "%Y%j")

# Subset variables to keep
to_keep = c(
  "date",
  "lon",
  "lat",
  "LST_Day_1km", 
  "QC_Day", 
  "Day_view_time", 
  # "Day_view_angl", 
  "LST_Night_1km", 
  "QC_Night", 
  "Night_view_time", 
  # "Night_view_angl", 
  "Emis_31", 
  "Emis_32"#, 
  # "Clear_day_cov", 
  # "Clear_night_cov",
)
dat = dat[, to_keep]

# Examine result
head(dat)

##         date         lon      lat LST_Day_1km QC_Day Day_view_time
## 1 2016-07-01 0.006481621 49.99583          NA      3            NA
## 2 2016-07-01 0.019444863 49.99583          NA      3            NA
## 3 2016-07-01 0.032408104 49.99583          NA      3            NA
## 4 2016-07-01 0.045371346 49.99583          NA      3            NA
## 5 2016-07-01 0.058334588 49.99583          NA      3            NA
## 6 2016-07-01 0.071297830 49.99583          NA      3            NA
##   LST_Night_1km QC_Night Night_view_time Emis_31 Emis_32
## 1            NA        3              NA      NA      NA
## 2            NA        3              NA      NA      NA
## 3            NA        3              NA      NA      NA
## 4            NA        3              NA      NA      NA
## 5            NA        3              NA      NA      NA
## 6            NA        3              NA      NA      NA

# Write table
write.fst(dat, gsub(".hdf", ".fst", filename, fixed = TRUE))

# End
end = Sys.time()

# Overall timing
end - start

## Time difference of 7.342018 secs