[gdal-dev] Bulding VRT files from HDF sds

[Jose Gomez-Dans]

Hi,
I want to put together a time series of MODIS datasets in VRT, the originals
being in HDF EOS format. I do that by defining the SDS for each date and
sticking it into a file:
HDF4_EOS:EOS_GRID:"../blablah.hdf":blah:blah
HDF4_EOS:EOS_GRID:"../blablah.hdf":blah:blah
HDF4_EOS:EOS_GRID:"../blablah.hdf":blah:blah
HDF4_EOS:EOS_GRID:"../blablah.hdf":blah:blah
[...]

I then use
    gdalbuildvrt -separate -input_files_list myfile.txt output_vrt.vrt

which creates the VRT file in output_vrt.vrt. The XML appears sensible, and
all things are as they should.

<VRTRasterBand dataType="Int16" band="32">
    <Metadata>
      <MDI key="STATISTICS_MINIMUM">-310</MDI>
      <MDI key="STATISTICS_MAXIMUM">32766</MDI>
      <MDI key="STATISTICS_MEAN">321.45436944444</MDI>
      <MDI key="STATISTICS_STDDEV">2990.0525001527</MDI>
    </Metadata>
    <NoDataValue>3.27670000000000E+04</NoDataValue>
    <ComplexSource>
      <SourceFilename
relativeToVRT="1">HDF4_EOS:EOS_GRID:&quot;../XXX&quot;:XXXX:XXX</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1200" RasterYSize="1200"
DataType="Int16" BlockXSize="1200" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="1200" ySize="1200"/>
      <DstRect xOff="0" yOff="0" xSize="1200" ySize="1200"/>
      <NODATA>32767</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="33">
    <NoDataValue>3.27670000000000E+04</NoDataValue>
    <ComplexSource>
      <SourceFilename
relativeToVRT="1">HDF4_EOS:EOS_GRID:&quot;../XXX&quot;:XXXX:XXX</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1200" RasterYSize="1200"
DataType="Int16" BlockXSize="1200" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="1200" ySize="1200"/>
      <DstRect xOff="0" yOff="0" xSize="1200" ySize="1200"/>
      <NODATA>32767</NODATA>
    </ComplexSource>
  </VRTRasterBand>


If I try to read this into python, using the GDAL bindings, it can deal with
individual bands, but using

g = gdal.Open ( VRT_file )
data = g.ReadAsArray()

it fails claiming that it can't find the original HDF datafile. this error
pops up with the band that corresponds to band number 33.

Running gdalinfo -stats on the new filename proceeds reporting the expected
information, band by band... slowly but surely, until it reaches band 33
(it's band 33 everytime, even with different datasets):

Band 33 Block=128x128 Type=Int16, ColorInterp=Undefined
ERROR 4: `HDF4_EOS:EOS_GRID:"../BlahBlah.hdf":blah:blah' does not exist in
the file system,
and is not recognised as a supported dataset
name.


ERROR 1: IReadBlock failed at X offset 0, Y offset 0

I can read individual bands OK (although I haven't tried looping through
them), including band 33. I have stacks of memory, and data size should not
be an issue.

Any clues?

Thanks!
Jose
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100416/bd3d8a19/attachment.html