LSA-SAF data reprojection with GDAL and R

Finally I found a solution using the Latitude and Longitude static files together with the DEM raster:

library(raster)
library(rgdal)

dem <- raster("HDF5_LSASAF_USGS_DEM_Euro.hdf")
lat <- raster("HDF5_LSASAF_MSG_LAT_Euro_4bytesPrecision.hdf")
lon <- raster("HDF5_LSASAF_MSG_LON_Euro_4bytesPrecision.hdf")

## Only the numeric content is needed
dem <- dem[]
## Scale lat and long values
lat <- lat[]/10000
lon <- lon[]/10000

## It is important to exclude 90º angles to avoid problems
## with reprojection. Besides, I reduce the data to a region
valid <- (lon < 10  & lon > -10) & (lat < 50 & lat > 30)
lat <- lat[valid]
lon <- lon[valid]
dem <- dem[valid]
dem[dem<0] <- NA

I define an SpatialPointsDataFrame object (irregularly spaced points) using the latitude and longitude values as the coordinates, the DEM file as the data, and with the long/lat projection.

projLL <- CRS('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0')
ll <- SpatialPointsDataFrame(cbind(lon, lat), data=data.frame(dem),
                             proj4string=projLL)

Next, I reproject this SpatialPointsDataFrame to the GEOS projection with spTransform.

projGeos <- CRS("+proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")
xy <- spTransform(ll, projGeos)

Finally, I define an SpatialPixelsDataFrame (a regular grid) with the new coordinates, which can be easily converted to a RasterLayer object.

demGrid <- SpatialPixelsDataFrame(xy, data.frame(dem), tolerance=0.311898)
dem <- raster(demGrid)

> dem
class       : RasterLayer 
dimensions  : 360, 592, 213120  (nrow, ncol, ncell)
resolution  : 2986.185, 2992.444  (x, y)
extent      : -883596.1, 884225.5, 3470127, 4547407  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
data source : in memory
names       : dem 
values      : 0, 3866  (min, max)