4

I am developing a python plugin for QGIS. The purpose of the plugin is to properly open an HDF5 dataset as a raster. Opening and displaying the raster can be done without the need for a plugin, however, QGIS does not recognize the origin (top left corner coordinates) of the dataset, or specific data file that I am importing. Top left corner coordinate information is embedded in the HDF5 metadata and I would like the plugin to recognize and feed it to QGIS.

When I run the code below through the plugin I get the raster loaded but not displayed, and the top left corner coordinate is (0, 0). How can I change that to whatever (x, y) values I want using PyQGIS syntax?

This is the code I am working with:

from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication, QFileInfo
from PyQt4.QtGui import QAction, QIcon, QFileDialog
import resources
from second_try_dialog import SecondTryDialog
import os.path
from osgeo import gdal, ogr, osr
from qgis.core import QgsRasterLayer, QgsMapLayerRegistry, QgsCoordinateReferenceSystem
from qgis.core import * #this is redundant with the line above
from qgis.gui import *

#### to open windows browser and search for HDF5 files only:
myfileNames = QFileDialog.getOpenFileNames(self.dlg,
                                            self.tr("HDF5 File Selector"), 
                                            "", 
                                            self.tr("HDF5 (*.hdf5 *.h5)"))
#### to loop through all selected HDF5 files, just handling one for now
for i in myfileNames:
    #### bind the HDF5 file
    myfile = QgsRasterLayer(i)

    #### to deal with multiple datasets within file
            if (len(myfile.subLayers()) > 1):
                #### to open dataset with desired data and name it Reflectance
                mydset = QgsRasterLayer(myfile.subLayers()[0], 'Reflectance')

                #### to set proper Coordinate Reference System (crs) info
                crs = mydset.crs()
                crs.createFromId(32611)
                mydset.setCrs(crs)

                #### to set raster bands to load as RGB
                mydset.renderer().setGreenBand(34)
                mydset.renderer().setRedBand(52)
                mydset.renderer().setBlueBand(18)

                #### trying to set top left coordinates
                #### This part isn't doing what I expected
                rect = QgsRectangle(326380.0,4103390.0-2853,326380.0+1332,4103390.0)
                mydset.setExtent(rect)

                #### to add selected dataset/raster to map canvas
                QgsMapLayerRegistry.instance().addMapLayer(mydset)

The result of calling this procedure is not what I expected. I want the raster be displayed in the map canvas and the top left corner to be located at (326380, 4103390). This is not what happens. The problem is that the raster doesn't load in the map canvas, even though the layer appears in the Layers Panel. The coordinates in the blank canvas seemed to have changed though. If I remove the two lines of code from the part that is not working, the raster does load, albeit with the top left corner at (0, 0).

Also, I keep getting an annoying message box that says: "CRS was undefined: defaulting to CRS EPSG...". I would like to resolve this problem as well.

How do I solve this?


Originally I thought the raster wasn't being displayed because the map canvas was blank but upon close inspection with the identify tool it became clear that the raster was there, but all the values for every band were overwritten with 'no data', making it appear blank. Also, since querying my raster's extent returns bottom left and top right corner coordinates, I have updated the values inside QgsRectangle so they correspond to those positions.

  • about CRS, what is the CRS shown with $> gdalinfo HDF5:<file>:<sublayer> ? – Luigi Pirelli Jul 1 '16 at 8:41
  • @LuigiPirelli Using gdalinfo I get a line that just says: Coordinate System is `' Also, I get 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). These coordinates seem strange because they do not correspond to my raster coordinates at all. – Adrian Jul 5 '16 at 20:42
  • are you talking about the HDF5 file or one of it's sublayers? – Luigi Pirelli Jul 7 '16 at 8:49
  • Open this file with gdal: f = gdal.Open('path/to/file.hdf5') and check this value: f.GetGeoTransform() it should return a tuple with 1st and 4th numbers different than 0.0 if origin is other than (0, 0) – dmh126 Jul 11 '16 at 16:20
1

I'm not 100% sure, but I think you may have to set the coordinate system of the canvas as well. There is a post here that indicates this: How to set CRS in QGIS using pyQGIS

For example, the answer there shows code like the following:

mycrs = QgsCoordinateReferenceSystem(32611)
iface.mapCanvas().mapRenderer().setDestinationCrs(mycrs)

I would suggest doing that first, before loading the sublayer.

  • Thanks for your response! I looked into it but QGIS seems to already set the CRS of the Map Canvas to the CRS of the layer. I double checked by investigating the project settings after importing the HDF5 dataset (albeit in the wrong location) through the plugin. – Adrian Jul 14 '16 at 20:28

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.