1

I have a set of HDF5 files, and I would like to cut and convert them to tiff files, but I am having troubles with the coordinate reference system. An example file is: "HDF5_LSASAF_MSG_LST_MSG-Disk_201901010600" if I give

gdalinfo "HDF5_LSASAF_MSG_LST_MSG-Disk_201901010600"

I get
enter image description here

enter image description here

So, if I try to convert it with gdal_translate, returning

gdal_translate -a_srs "WGS84" -of GTiff NETCDF:"HDF5_LSASAF_MSG_LST_MSG-Disk_201901010600"://LST LST_201901011200.tif

Note that I am forced to write NETCDF: instead of HDF5:. Only with the first option the system recognizes the file as existing. To write NETCDF: come to me looking at the gdalinfo output, which returns

SUBDATASET_1_NAME=NETCDF:"HDF5_LSASAF_MSG_LST_MSG-Disk_201901010600":LST

So, I get the tif file, but it has a strange coordinate system. For example, if I try to open it in QGIS, I see this
enter image description here

enter image description here

enter image description here

enter image description here

it seems to be a small square and does not have synchronicity with the OSM standard basemap, although they have the same Coordinate Reference System (CRS), EPSG:3857 - WGS 84 / Pseudo-Mercator.

Does anybody know why this happens?

You can download an example file here.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Andrea
  • 69
  • 3
  • Are you aware of this? https://gis.stackexchange.com/a/392388/88814 and https://gis.stackexchange.com/a/383437/88814 – Babel Dec 19 '22 at 16:10
  • That looks like an azimuthal projection showing about half the earth (like a view from a high orbit satellite). Much like this: https://stackoverflow.com/questions/53643062/lsa-saf-satellite-hdf5-plot-in-python-cartopy - for some reason the file you link isn't recognised as HDF5 by my gdal or command-line hdfview... – Spacedman Dec 19 '22 at 16:48
  • @Spacedman in fact, the command gdalinfo recognizes it as NETCDF, returning SUBDATASET_1_NAME=NETCDF:"HDF5_LSASAF_MSG_LST_MSG-Disk_201901010600":LST I don't know why – Andrea Dec 20 '22 at 09:53
  • HDF and NetCDF are closely related. I've a nasty feeling the download is getting corrupted. How long in bytes is the file you are using? My download is 82629665 exactly. Looking at it in binary the first few bytes are subtly different to an HDF5 file I've got that does work... – Spacedman Dec 20 '22 at 12:37
  • @Spacedman the file is correct. I think the problem is also that my gdal installation does not have the HDF5 library: in fact, gdalinfo --formats does not return "HDF5" within the list, and maybe for this reason it recognizes it as NETCDF, making a mistake and not correctly processing it with gdal_translate. So, I think I have to install this but I don't know how on Windows... – Andrea Dec 20 '22 at 13:54
  • I'm pretty sure what I'm downloading isn't correct. Looking at the HDF5 spec https://docs.hdfgroup.org/hdf5/develop/_f_m_t11.html#FileMetaData I don't see the \r \n sequence, I only see a \n, which is a sign that the file has been treated as text at some point and had line-endings converted. An HDF5 I have that does work for me does have the header format from the spec in the link. Have you got an original source with a link to that data? – Spacedman Dec 20 '22 at 14:10
  • I can get a file with the same name here: https://datalsasaf.lsasvcs.ipma.pt/PRODUCTS/MSG/MLST/HDF5/2019/01/01/ (login required) which looks like yours but is much smaller (6Mb vs 82Mb). It reads into QGIS the right way round but has no coordinate reference so its not in the right place. You probably no need a GDAL with the HDF5 driver, which might be a pain to do in Windows unless someone has compiled it all into a version for Windows... – Spacedman Dec 20 '22 at 14:27
  • With the one I've downloaded from the site I linked to, after two gdal commands (a gdal_translate and a gdalwarp) I get a perfectly registered file in EPSG:4326 lat-long. Even the metadata listing you give is the same as mine though, but the file you link to download is broken for me! – Spacedman Dec 20 '22 at 15:25

1 Answers1

1

I solved the problem by using the following python script:

import os

filename = "HDF5_LSASAF_MSG_LST_MSG-Disk_201901010000"
if os.path.exists(filename):
    print("File exists, ok.\n")
else:
    print("ERROR: file does NOT exist!\n")
    raise SystemExit


f_out = filename + ".tif"
f_rep = filename + "_rep.tif"
print("%s\n%s\n" % (f_out, f_rep))


os.system('gdal_translate -of GTiff -a_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs"

 -a_ullr -5568000 5568000 5568000 -5568000  HDF5:"'+ filename + '"://LST '+ f_out)


Mapping the new values and filling the new _rep.tif file


os.system('gdalwarp -t_srs EPSG:4326 -wo SOURCE_EXTRA=100 ' + f_out + ' ' + f_rep)
Andrea
  • 69
  • 3