5

Yesterday i downloaded an Aster L1B image from USGS's website. If i load it in QGIS, it doesn't get displayed properly(location wise). Even after i've enabled "On the fly" option, it's not reading the WGS84 coordinate system (Figured out this type of glitch happened in ArcGIS also - How to import and project ASTER 1B data into Arcmap?). I tried to assign projection to the image but it didn't work (Got the message - hdf has no raster bands).

I tried the same in GRASS(r.in.gdal) and it got "Error: Selected band (1) doesn't exist"

But if i open the same image in Global Mapper, the image gets displayed perfectly.

How to solve this ? enter image description here

Improve this question
8
  • 1
    Can you add how you added the data? From the topic you linked I read that the file contains several subdatasets. It might be necessary to specify the subdataset you want, or extract it with gdalwarp in a first step. Running gdalinfo on the file might give further clues.
    – AndreJ
    Commented Sep 25, 2014 at 5:50
  • In QGIS, loaded just like a normal Raster layer -> Selected all the 15 subdatasets (4 VNIR, 6 SWIR, 5 TIR) In GRASS, as mentioned in my question using r.in.gdal. There is a command available in grass - r.in.aster. But from the description // Georeference, rectify and import Terra-ASTER imagery and relative DEM's using gdalwarp // i don't think(or know) it will serve my purpose.
    – joseph_k
    Commented Sep 25, 2014 at 6:22
  • The subdatasets might have different CRS or extent, so better try only one. From the image it looks like the raster is loaded without proper CRS information. gdalwarp can add CRS and extent if needed to make a raster file compatible with QGIS.
    – AndreJ
    Commented Sep 25, 2014 at 6:48
  • Tried both, gdalwarp & Assign Projection. As i stated in my question, if i run those two commands all i am getting is this one - hdf has no raster bands
    – joseph_k
    Commented Sep 25, 2014 at 7:02
  • But how come Global_Mapper can read the data properly ? That's what puzzling me.
    – joseph_k
    Commented Sep 25, 2014 at 7:02

2 Answers 2

Reset to default
7

The ASTER L1B files contain several subdatasets with different resolutions. That's why you can not easily add them to QGIS. You have to run gdalinfo and gdalwarp on it to get a tif file that QGIS can import:

gdalinfo AST_L1B.hdf >>info.txt

gives you a long list of metadata. Look out for the subdatasets:

Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData1
  SUBDATASET_1_DESC=[4200x4980] ImageData1 VNIR_Swath (8-bit unsigned integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData2
  SUBDATASET_2_DESC=[4200x4980] ImageData2 VNIR_Swath (8-bit unsigned integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData3N
  SUBDATASET_3_DESC=[4200x4980] ImageData3N VNIR_Swath (8-bit unsigned integer)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData3B

and take the full name of the first one:

gdalinfo HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData1 >>sds1.txt

This gives you detailed information of the CRS and extent of that subdataset:

LOWERLEFT=63.4513647453478, -151.617477277645
LOWERRIGHT=63.2162350711876, -150.220604276473
UPPERLEFT=63.9822445743476, -151.179093369548
UPPERRIGHT=63.7426354293633, -149.760000938454
SRS=PROJCS["UTM Zone 5, Northern Hemisphere",GEOGCS["Unknown datum based upon the GRS 1980 ellipsoid",DATUM["Not specified (based on GRS 1980 spheroid)",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

Next we transform it to a Geotiff:

gdalwarp HDF4_EOS:EOS_SWATH:"AST_L1B.hdf":VNIR_Swath:ImageData1 sds.tif
Creating output file that is 6063P x 5573L.
Processing input file HDF4_EOS:EOS_SWATH:AST_L1B.hdf:VNIR_Swath:ImageData1.
0...10...20...30...40...50...60...70...80...90...100 - done.

which can be loaded with corrcet CRS information and extent by QGIS. You have to do that with every subdataset you want.

Improve this answer
2
  • Thanks a lot Andre. Worked like a charm. One last question. We can run this in batch mode rite ? Through modeler..or through script
    – joseph_k
    Commented Sep 25, 2014 at 10:40
  • You have to get the exact name of the subdatasets (case-sensitive!). I don't know if you can extract these automatically. Adding the folder name is another problem if not running in a DOS shell within the same folder.
    – AndreJ
    Commented Sep 25, 2014 at 10:43
1

Another way to work with the images using the GUI is:

  1. First in a new project add the .hdf as a raster and select which band do you want to work with. QGIS will set a new user defined CRS with a number, for example in my case USER:100008.
  2. Go to Raster>Projections>Warp(Reproject). You will reproject each one of the bands given an different output name for each one. Select the USER:100008 as source SRS, and as target select the corresponding UTM projection. In my case EPSG:32719.
  3. After having a .tif file for each band you can do a layer stacking with them. I tried once and its possible to stack in a single process step the band 1,2 and 3 with a band 4 that is in a 60m resolution. It will resample the band 4 image to match the resolution of the other layers.

Conclusion: there is not a straight way to work with ASTER image in QGIS yet. A better solution might be to order ASTER image in .tif, I know it is possible through REVERB data center. Now that ASTER images and product are free available, it is the best option I think.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.