3

I'm definitely not the most experienced GDAL user, but I already struggle for quite a while now getting my head around the hdf4 file format in GDAL. The thing is that I have a directory with many hdf4 files from Sentinel-2, all from the same day. In the end I want to build a mosaic of all the scenes. It's just that these scenes are from South Africa and when I create a VRT from all files and then use gdal_translate to make it a GeoTIFF it's all in the UTM zone 34N, what apparently doesn't make too much sense as it should be in 34S.

When I do gdalinfo on one of the .hdf files I get the following confusing output (sorry for the size of that):

gdalinfo HLS.S30.T34JDN.2018301.v1.4.hdf 
Driver: HDF4/Hierarchical Data Format Release 4
Files: HLS.S30.T34JDN.2018301.v1.4.hdf
Size is 512, 512
Coordinate System is `'
Metadata:
  ACCODE=LaSRCS2AV3.5.5
  AngleBand=0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
  arop_ave_xshift(meters)=0
  arop_ave_yshift(meters)=0
  arop_ncp=0
  arop_rmse(meters)=0
  arop_s2_refimg=NONE
  cloud_coverage=1
  DATASTRIP_ID=S2B_OPER_MSI_L1C_DS_MPS__20181028T120816_S20181028T084332_N02.06
  HLS_PROCESSING_TIME=2018-11-01T23:32:43Z
  HORIZONTAL_CS_CODE=EPSG:32734
  HORIZONTAL_CS_NAME=WGS84 / UTM zone 34S
  L1C_IMAGE_QUALITY=NONE
  L1_PROCESSING_TIME=2018-10-28T12:47:06.529844Z
  MEAN_SUN_AZIMUTH_ANGLE(B01)=57.9076150593934
  MEAN_SUN_ZENITH_ANGLE(B01)=27.018691721451
  MEAN_VIEW_AZIMUTH_ANGLE(B01)=177.57328049295
  MEAN_VIEW_ZENITH_ANGLE(B01)=3.06199611054485
  MSI band 01 bandpass adjustment slope and offset=0.995900, -0.000200
  MSI band 02 bandpass adjustment slope and offset=0.977800, -0.004000
  MSI band 03 bandpass adjustment slope and offset=1.007500, -0.000800
  MSI band 04 bandpass adjustment slope and offset=0.976100, 0.001000
  MSI band 11 bandpass adjustment slope and offset=1.000000, -0.000300
  MSI band 12 bandpass adjustment slope and offset=0.986700, 0.000400
  MSI band 8a bandpass adjustment slope and offset=0.996600, 0.000000
  NBAR_Solar_Zenith=43.7861821326806
  NCOLS=3660
  NROWS=3660
  PROCESSING_BASELINE=02.06
  PRODUCT_URI=S2B_MSIL1C_20181028T082039_N0206_R121_T34JDN_20181028T120816.SAFE
  SENSING_TIME=2018-10-28T08:43:32.629Z
  SPACECRAFT_NAME=Sentinel-2B
  spatial_coverage=100
  SPATIAL_RESOLUTION=30
  TILE_ID=S2B_OPER_MSI_L1C_TL_MPS__20181028T120816_A008584_T34JDN_N02.06
  ULX=399960
  ULY=-3199980
Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B01
  SUBDATASET_1_DESC=[3660x3660] B01 Grid (16-bit integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B02
  SUBDATASET_2_DESC=[3660x3660] B02 Grid (16-bit integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B03
  SUBDATASET_3_DESC=[3660x3660] B03 Grid (16-bit integer)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B04
  SUBDATASET_4_DESC=[3660x3660] B04 Grid (16-bit integer)
  SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B05
  SUBDATASET_5_DESC=[3660x3660] B05 Grid (16-bit integer)
  SUBDATASET_6_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B06
  SUBDATASET_6_DESC=[3660x3660] B06 Grid (16-bit integer)
  SUBDATASET_7_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B07
  SUBDATASET_7_DESC=[3660x3660] B07 Grid (16-bit integer)
  SUBDATASET_8_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B08
  SUBDATASET_8_DESC=[3660x3660] B08 Grid (16-bit integer)
  SUBDATASET_9_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B8A
  SUBDATASET_9_DESC=[3660x3660] B8A Grid (16-bit integer)
  SUBDATASET_10_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B09
  SUBDATASET_10_DESC=[3660x3660] B09 Grid (16-bit integer)
  SUBDATASET_11_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B10
  SUBDATASET_11_DESC=[3660x3660] B10 Grid (16-bit integer)
  SUBDATASET_12_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B11
  SUBDATASET_12_DESC=[3660x3660] B11 Grid (16-bit integer)
  SUBDATASET_13_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:B12
  SUBDATASET_13_DESC=[3660x3660] B12 Grid (16-bit integer)
  SUBDATASET_14_NAME=HDF4_EOS:EOS_GRID:"HLS.S30.T34JDN.2018301.v1.4.hdf":Grid:QA
  SUBDATASET_14_DESC=[3660x3660] QA Grid (8-bit unsigned integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Could someone who has a little more experience in working with this file format give a hint on what to do with that?

I used python, and my commands were: build_vrt = gdal.BuildVRT("final.vrt", "HLS.S30.T34JDN.2018301.v1.4.hdf ") final_tif = gdal.Translate("final.tiff", "final.vrt")

  • Please, include the commands used. – Gabriel De Luca Jan 29 at 12:59
  • The comment is: gdalinfo HLS.S30.T34JDN.2018301.v1.4.hdf – Robin Kohrs Jan 30 at 11:04
  • Hi Robin, please include in the post the commands used to create VRT and GeoTIFF files. – Gabriel De Luca Jan 30 at 11:15
  • Oh yes of course, sorry: I used python, and my commands were: build_vrt = gdal.BuildVRT("final.vrt", "HLS.S30.T34JDN.2018301.v1.4.hdf ") final_tif = gdal.Translate("final.tiff", "final.vrt") – Robin Kohrs Jan 30 at 11:39
  • Can you provide the output of running gdalinfo on your final tiff? – Ali Feb 2 at 16:06
1

Regarding the use of HDF format with GDAL the key concept is that the file contains subdatasets as opposed to raster bands.

With a GeoTiff, for example, you can open it, and pick a band to read

f = ''
ds = gdal.Open(f)

if ds is None:
    raise Exception('Failed to open dataset')
b = ds.GetRasterBand(1)

# etc...

When using a HDF file, the file itself does not contain raster bands - it's only really a pointer to the subdatasets within. To interrogate a subdataset, it first has to be opened, similar to above.

f = ''
ds = gdal.Open(f) # ds is now the HDF file

# Error catching here..

print(ds.GetMetdata()) # This will print the metadata you see in your gdalinfo output

# The GetSubDatasets method returns a list of all the subdatasets in the HDF file
# Each element is a tuple object of the form (subdataset_path,subdataset_desc)
sds_list = ds.GetSubDatasets()

# Let's see what the subdatasets are...
for sds in sds_list:
    print(sds[1]) # i.e. the desc

# Say the first subdataset corresponds to reflectance in the blue band and we now
# want to access the data. We open the subdataset as we would any other dataset
blue_sds = sds_list[0] # this is a tuple (path,desc)
blue = gdal.Open(blue_sds[0]) # open using the path

prj = blue.GetProjection()
bnd = blue.GetRasterBand(1)

# etc...
| 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.