4

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.

2

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)

dem with spplot

0

Unfortunately, your test example is not online anymore, so I took a gobal Land-Sea-Temperature sample from ftp://landsaf.ipma.pt/pub/Example_of_Product/

As mentioned in my answer to Transforming geostationary satellite image to lon/lat, the size of global geos reprojected files is 3712x3712 pixels (3-km-resolution), but the actual data array is smaller. GDAL knows the GEOS projection, but has a different resolution using meters, depending on the satellite height above the ellipsoid. But you can add the proper extent expected by GDAL, and then reproject to WGS84 with:

gdal_translate -a_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs" -a_ullr -5568000 5568000 5568000 -5568000 HDF5:"HDF5_LSASAF.hdf"://LST temp.tif
gdalwarp -t_srs EPSG:4326 -wo SOURCE_EXTRA=100 temp.tif LSASAF.tif

to get this picture:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.