7

The ICARE Data and Services Center distributes the SEV_AERUS-AEROSOL-D3 product in HDF5 format (http://www.icare.univ-lille1.fr/archive/?dir=GEO/MSG+0000/SEV_AERUS-AEROSOL-D3.v1.03/).

The main component of the HDF file is a 3712 x 3712 matrix, which represents measurements on a geostationary view of the earth, such as in the following image:

In R, I am able to read the latter into a matrix object named m using the h5 package:

library(h5)
setwd("/home/michael/Dropbox/BGU/Alex/HDF5")
filename = "SEV_AERUS-AEROSOL-D3_2005-01-01_V1-03.h5"
h = h5file(filename)
m = h["ANGS_06_16"][]

Here are the matrix dimensions:

dim(m)
## [1] 3712 3712

And, for example, these are the values in the first 1-5 rows/columns:

m[1:5, 1:5]
##        [,1]   [,2]   [,3]   [,4]   [,5]
## [1,] -32768 -32768 -32768 -32768 -32768
## [2,] -32768 -32768 -32768 -32768 -32768
## [3,] -32768 -32768 -32768 -32768 -32768
## [4,] -32768 -32768 -32768 -32768 -32768
## [5,] -32768 -32768 -32768 -32768 -32768

The HDF5 file also provides metadata with the lon/lat of two corners for the bounding box of m (below) as well as the projection type (specified simply as "GEOS"). Here is a reproduction of the bounding box coordinates exactly as given in the metadata:

library(raster)
ext = new("Extent"
          , xmin = -75.1001529854046
          , xmax = 75.1001529854046
          , ymin = -77.9875583358503
          , ymax = 77.9875583358503
)
ext
## class       : Extent 
## xmin        : -75.10015 
## xmax        : 75.10015 
## ymin        : -77.98756 
## ymax        : 77.98756

The problem is that I need to find out the spatial location of each measurement in the matrix. In order to do that, the matrix needs to be converted into a point layer with lon/lat coordinates for each pixel.

My idea was:

  1. Convert the latter bounding box coordinates of m to "GEOS" projection coordinates.

  2. Convert m to a RasterLayer object, using the "GEOS" projection bounding box and a "GEOS" CRS definition.

  3. Transform the resulting RasterLayer to a point layer.

  4. Reproject the point layer to lon/lat.

Unfortunately, I am stuck in step (1). For example, the following code takes the Southernmonst/Westernmost corner in order to calculate its coordinates in the "GEOS" system:

ext = data.frame(x = ext[1], y = ext[3])
ext
##           x         y
## 1 -75.10015 -77.98756
coordinates(ext) = ~ x + y
proj4string(ext) = '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0'
spTransform(ext, "+proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")

The last expression gives the following error message:

Error in spTransform(x, CRS(CRSobj), ...) : 
error in pj_transform: tolerance condition error

In other words, the (-75.10015, -77.98756) point can't be converted from lon/lat to "GEOS". I also tried different variations of the "GEOS" PROJ.4 string, none worked. The "tolerance condition error" is apparently a general PROJ.4 issue, since I also get the same error when trying to reproject the (-75.10015, -77.98756) point from lon/lat to "GEOS" using pyproj in Python, as suggested here.

When using coordinates closer to the equator, such as (10,20), the code works fine (no "tolerance condition error"). But this doesn't help, since I can't subset the m matrix without knowing the new bounding box coordinates.

Another option could have been translating the matrix to lon/lat using Longitute/Latitude static files, as suggested here. However I couldn't find any such files matching the 3712 x 3712 dimensions.

Does anyone have experience with going from a 3712 x 3712 matrix representing a geostationary satellite image to a point layer in lon/lat? Are there any static files providing predetermined lon/lat values of each cell in such a grid?

|improve this question
  • 1
    I suggst you to use python. Piproj. First you have to transform lines and columns in the image intermediate coordinates as explained. Or maybe your file contain it : search Gdal array.... and gdal projection string inside the file........ – J L Lliso Valverde Nov 9 '18 at 19:35
13

Your attempt is designed to fail. If you look at the image, you see the data arranged as a circle, with black triangles in the corners of the square, where the satellite view goes right into orbit. In your test data, you see only NODATA -32768 for those parts of the image.

The extent is between +/-75 and +/- 78, but these values are only reached in the middle of the egdes. So you can not reproject those black triangles to Earth surface coordinates.

UPDATE

The Metadata of the HDF file reveals some mysteries:

Altitude=42164 
Ancillary_Files=MSG+0000.3km.lat

So the satellite height is the same as mentioned in http://geotiff.maptools.org/proj_list/geos.html, and I assume they took the same ellipsoid (not exactly WGS84).

With the help of http://www.cgms-info.org/documents/pdf_cgms_03.pdf and http://publications.jrc.ec.europa.eu/repository/bitstream/JRC52438/combal_noel_msg_final.pdf, I found that the size of 3712px is not the real extent covered by the data. The size provides a scanning angle of the satellite of about +/-8.915 degree, but the angle that was used is smaller.

Proj.4 calculates the extent by multiplying the satellite's scanning angle by the height above ground (see http://proj4.org/projections/geos.html). So with a bit of try, an extent of +/- 5568000m (3712*3000m/2 or 8.915*pi()/180*35785831m) fits to the 3712px used in the 3-km-resolution hdf.

So the correct translation commands are:

gdal_translate -a_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs" -a_ullr -5568000 5568000 5568000 -5568000 HDF5:"SEV_AERUS-AEROSOL-D3_2006-01-01_V1-03.h5"://ANGS_06_16 temp.tif
gdalwarp -t_srs EPSG:4326 -wo SOURCE_EXTRA=100 temp.tif output.tif

And the result looks good:

enter image description here

As an alternative, you can take the lat and lon subdatasets from http://www.icare.univ-lille1.fr/archive/?dir=GEO/STATIC/ in file MSG+0000.3km.hdf

|improve this answer
  • Thank you very much for the quick answer. The explanation makes sense: if the bounding box corners are outside of the "GEOS" grid then there's no way I can find their coordinates in that CRS. The gdalwarp command also works, but doesn't seem to find the right transformation either. The resulting image looks exacly like the origin, only with a different bounding box, in lon/lat: 0.0000000000000000,0.0000000010817455 0.0000089843362493,0.0000090436947705 which is clearly incorrect... – Michael Dorman May 8 '16 at 7:51
  • See my extended answer. – AndreJ May 8 '16 at 9:46
  • I don't yet understand the required geometric calculation to get exact geo-referencing with this method, so not sure I am going to use it over the subsetting method (unless I must get data for an area that was not predefined in the static files). But in principle this is the result I wanted, so I'm marking it as the right answer. Thank you very much! – Michael Dorman May 8 '16 at 10:13
  • Great! Thank you for the update, this is what I was looking for. – Michael Dorman May 14 '16 at 5:45
  • I'd like to ask the rationale behind the +/- 5568000m value? I'm sure it's obvious but I can't find any reference to it? – tda Dec 19 '18 at 15:48
2

UPDATE:

The above-mentioned approach by Oscar Perpiñán found here seems to work, if I only take a subset of the matrix the way Oscar did.

For that to work I needed to know the subset boundaries and the static files, both available online for a few regions (such as North Africa). 

m = m[700:1850, 1240:3450, drop = FALSE]

lon = raster("NAFR_LON.img")
lat = raster("NAFR_LAT.img")

lat <- lat[]
lon <- lon[]

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

writeOGR(ll, ".", "output_eu", driver = "ESRI Shapefile", overwrite_layer = TRUE)

This is what the resulting layer looks like in QGIS:

Nevertheless, if there is a way to transform the entire matrix, as asked in the original question, I'll be happy to learn about it...

UPDATE: User AndreJ provides the complete solution, above.

|improve this answer

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.