Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal

The answer is surprisingly simple.

sgr_lst <- readGDAL(sds, as.is = TRUE)

solves the issue. The only thing left to do then is transform the resulting SpatialGridDataFrame to a proper Raster* object using raster(). For this purpose, it is necessary to retrieve the west, east, south, and north bounding coordinates (see ?extent: xmin, xmax, ymin, ymax) from the metadata e.g. via

meta <- attr(gdalinfo, "mdata")

## string patterns
crd_str <- paste0(c("WEST", "EAST", "SOUTH", "NORTH"), "BOUNDINGCOORDINATE")
## search for patterns in metadata
crd_id <- sapply(crd_str, function(i) grep(i, meta))
## extract information
crd <- meta[crd_id]
## create 'extent' object
crd <- as.numeric(sapply(strsplit(crd, "="), "[[", 2))
ext <- extent(crd)

The thus determined extent (together with the EPSG code) can then be passed on to the finally created RasterLayer which, after applying the scale factor and offset, definitely looks fine now.

## rasterize, apply offset and scale factor
rst_lst <- (raster(sgr_lst) - -1.5e+04) * 1.0e-02
## set extent and coordinate reference system
extent(rst_lst) <- crd
projection(rst_lst) <- "+init=epsg:4326"

I have not managed to do it in R either.. and I have also spent countless hours. What I do now is this:

1. Use the ModisReprojectionTool to extract the layers and the
   subsets I need as binary files
2. read the binary files and if necessary convert them to "raster" objects. I mainly use them as matrix and in the end transform them into raster-objects to write them as TIFF

the code to run MRT:

# Reproject from HDF to plain binary with MRT
setwd('/Volumes/Archive1/MOD11C3V5/')
DIR <- getwd()
# Run the Modis Reprojection Tool once with one HDF-file and save the desired parameterization in a .prm-file
# Use this .prm-file here:
ReprojectionParamter <- 'Pamir0.05_binray.prm'

# Getting the file list that you want to process
FileList <- list.files()
FileListHDF <- FileList[which(regexpr(pattern='hdf$',FileList)>0)] # only the hdf-files but not the hdf.xml files

### Setting environmental variables for MRT_DATA_DIR
Sys.setenv(MRT_DATA_DIR='/Volumes/DATA/ModisReprojectionTool/data')
Sys.setenv(MRTDATADIR='/Volumes/DATA/ModisReprojectionTool/data')

for (i in FileListHDF){
  system(command=paste('/Volumes/DATA/ModisReprojectionTool/bin/resample -p ',ReprojectionParamter,' -i ',DIR,'/',i,' -o ',DIR,'/binary_small/',i,'.hdr',sep=''),wait=T,)
}

Then you have to set the "what" argument of readBin to int, numeric, etc... I always have the same extent of my files (defined in the "prm"-file) I get the extent and resolution directly from there:

ReprojectionParamter <- scan('Pamir0.05_binary2.prm', nmax=90,what='character')
SpatExtent.minLon <- as.numeric(ReprojectionParamter[31])
SpatExtent.minLat <- as.numeric(ReprojectionParamter[36])
SpatExtent.maxLon <- as.numeric(ReprojectionParamter[37])
SpatExtent.maxLat <- as.numeric(ReprojectionParamter[30])
SpatExtent.RES <- as.numeric(ReprojectionParamter[72])
# readBinary function
UInt8_LST <- function(f, ...) readBin(f, what = "integer", signed = FALSE, endian = "little", size = 2, ...)
# Finally read the data; Nlat,Nlon etc you can calculate easily from the info in the prm-file; 0.02 was my scale factor
NC.LST_night <- matrix(UInt8_LST(f=LST_night.filename, n= Nlon*Nlat)),nrow=Nlat,ncol=Nlon,byrow=T)*0.02-273.15

Maybe it helps.