Friday 16 June 2017

qgis - How to properly open an ASTER data?


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



Answer




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.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...