Trabajar con archivos HDF4 en GDAL

Question

Definitivamente no soy el usuario más experimentado de GDAL, pero ya hace bastante tiempo que me cuesta hacerme a la idea del hdf4 en GDAL. El caso es que tengo un directorio con muchos hdf4 archivos de Sentinel-2, todos del mismo día. Al final quiero construir un mosaico de todas las escenas. Lo que pasa es que estas escenas son de Sudáfrica y cuando creo un VRT a partir de todos los archivos y luego uso gdal_translate para hacerlo un GeoTIFF está todo en la zona UTM 34N, lo que aparentemente no tiene mucho sentido ya que debería estar en 34S.

Cuando lo haga gdalinfo en uno de los .hdf obtengo el siguiente resultado confuso (perdón por el tamaño):

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)

¿Podría alguien que tenga un poco más de experiencia en el trabajo con este formato de archivo dar una pista sobre qué hacer con eso?

Usé python, y mis comandos fueron: build_vrt = gdal.BuildVRT("final.vrt", "HLS.S30.T34JDN.2018301.v1.4.hdf ") final_tif = gdal.Translate("final.tiff", "final.vrt")

Answer 1

En cuanto al uso del formato HDF con GDAL, el concepto clave es que el archivo contiene subconjuntos de datos en lugar de bandas ráster.

Con un GeoTiff, por ejemplo, puedes abrirlo y elegir una banda para leerla

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

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

# etc...

Cuando se utiliza un archivo HDF, el archivo en sí no contiene bandas de trama, sino que sólo es un puntero a los subconjuntos de datos que contiene. Para consultar un subconjunto de datos, primero hay que abrirlo, de forma similar a la descrita anteriormente.

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