LSA-SAF data reprojection with GDAL and R

The Land Surface Analysis Satellite Applications Facility (LSA SAF) distributes several products derived from the EUMETSAT satellites data. I am working with R to read and process the data (the file of this example is available from the Static Auxiliary Data page of the LSA-SAF service.)

library(rgdal)
library(raster)


r <- raster("HDF5_LSASAF_USGS_DEM_Euro.hdf")

Although the file is read correctly, the projection info is missing.

> r
class       : RasterLayer 
dimensions  : 651, 1701, 1107351  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 1701, 0, 651  (xmin, xmax, ymin, ymax)
coord. ref. : NA 

data source : /home/oscar/temp/HDF5_LSASAF_USGS_DEM_Euro.hdf 
names       : HDF5_LSASAF_USGS_DEM_Euro 
values      : -32768, 32767  (min, max)

If I use gdalsrsinfo to parse the file and obtain the projection, I get ERROR 1: ERROR - failed to load SRS definition from ... which is more or less the same if I use GDALspatialRef in R:

> GDALSpatialRef('HDF5_LSASAF_USGS_DEM_Euro.hdf')
[1] ""

As far as I know these files use the GEOS projection. However, if I set this PROJ.4 string and the file is transformed to a latitude-longitude projection, the results are erroneous (look at the extent of the projectExtent output):

projection(r) <- CRS("+proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")


> projectExtent(r, CRS(projLL))
class       : RasterLayer 
dimensions  : 651, 1701, 1107351  (nrow, ncol, ncell)
resolution  : 8.983153e-06, 9.043695e-06  (x, y)
extent      : 0, 0.01528034, 0, 0.005887445  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

It seems that the LSA-SAF data are not yet geo-corrected and the images are expressed in the "raw" Column-Row system, so this could be the root of the problem. There are Latitude and Longitude static files which could be used for remapping the rasters. However, I would like to use a more direct approach with a correct PROJ.4 string and a projection transform with gdalwarp or rgdal::spTransform in R.

Answer

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)