Skip to content
Snippets Groups Projects
utils.py 53.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • # -*- coding: utf-8 -*-
    """
    /***************************************************************************
    
                                  -------------------
            begin                : 2015-08-20
            git sha              : $Format:%H$
            copyright            : (C) 2017 by HU-Berlin
            email                : benjamin.jakimow@geo.hu-berlin.de
     ***************************************************************************/
    
    /***************************************************************************
     *                                                                         *
     *   This program is free software; you can redistribute it and/or modify  *
     *   it under the terms of the GNU General Public License as published by  *
     *   the Free Software Foundation; either version 2 of the License, or     *
     *   (at your option) any later version.                                   *
     *                                                                         *
     ***************************************************************************/
    """
    # noinspection PyPep8Naming
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    import os, sys, math, re, io, fnmatch, uuid
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    from qgis.core import *
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    from qgis.gui import *
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    import qgis.utils
    
    from PyQt5.QtCore import *
    
    from PyQt5.QtWidgets import *
    
    from PyQt5.QtGui import *
    from PyQt5.QtXml import QDomDocument
    from PyQt5 import uic
    
    from timeseriesviewer import DIR_UI, DIR_REPO
    
    from timeseriesviewer import messageLog
    
    import timeseriesviewer
    MAP_LAYER_STORES = [QgsProject.instance()]
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    def qgisInstance():
        """
        If existent, returns the QGIS Instance.
        :return: QgisInterface | None
        """
    
        from timeseriesviewer.main import TimeSeriesViewer
        if isinstance(qgis.utils.iface, QgisInterface) and \
            not isinstance(qgis.utils.iface, TimeSeriesViewer):
            return qgis.utils.iface
        else:
            return None
    
    
    def findMapLayer(layer)->QgsMapLayer:
        """
        Returns the first QgsMapLayer out of all layers stored in MAP_LAYER_STORES that matches layer
        :param layer: str layer id or layer name or QgsMapLayer
        :return: QgsMapLayer
        """
        if isinstance(layer, QgsMapLayer):
            return layer
        elif isinstance(layer, str):
            #check for IDs
            for store in MAP_LAYER_STORES:
                l = store.mapLayer(layer)
                if isinstance(l, QgsMapLayer):
                    return l
            #check for name
            for store in MAP_LAYER_STORES:
                l = store.mapLayersByName(layer)
                if len(l) > 0:
                    return l[0]
        return None
    
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    def file_search(rootdir, pattern, recursive=False, ignoreCase=False, directories=False, fullpath=False):
    
        """
        Searches for files
        :param rootdir: root directory to search for files.
        :param pattern: wildcard ("my*files.*") or regular expression to describe the file name.
        :param recursive: set True to search recursively.
        :param ignoreCase: set True to ignore character case.
        :param directories: set True to search for directories instead of files.
        :return: [list-of-paths]
        """
    
        assert os.path.isdir(rootdir), "Path is not a directory:{}".format(rootdir)
        regType = type(re.compile('.*'))
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        for entry in os.scandir(rootdir):
    
            if directories == False:
                if entry.is_file():
                    if fullpath:
                        name = entry.path
                    else:
                        name =  os.path.basename(entry.path)
                    if isinstance(pattern, regType):
                        if pattern.search(name):
                            yield entry.path.replace('\\','/')
    
                    elif (ignoreCase and fnmatch.fnmatch(name, pattern.lower())) \
                            or fnmatch.fnmatch(name, pattern):
                        yield entry.path.replace('\\','/')
                elif entry.is_dir() and recursive == True:
                    for r in file_search(entry.path, pattern, recursive=recursive, directories=directories):
                        yield r
            else:
                if entry.is_dir():
                    if recursive == True:
                        for d in file_search(entry.path, pattern, recursive=recursive, directories=directories):
                            yield d
                            
                    if fullpath:
                        name = entry.path
                    else:
                        name = os.path.basename(entry.path)
                    if isinstance(pattern, regType):
                        if pattern.search(name):
                            yield entry.path.replace('\\','/')
    
                    elif (ignoreCase and fnmatch.fnmatch(name, pattern.lower())) \
                            or fnmatch.fnmatch(name, pattern):
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
                        yield entry.path.replace('\\','/')
    
    NEXT_COLOR_HUE_DELTA_CON = 10
    NEXT_COLOR_HUE_DELTA_CAT = 100
    def nextColor(color, mode='cat'):
        """
        Returns another color
        :param color: the previous color
        :param mode: 'cat' - for categorical color jump (next color looks pretty different to previous)
                     'con' - for continuous color jump (next color looks similar to previous)
        :return:
        """
        assert mode in ['cat','con']
        assert isinstance(color, QColor)
        hue, sat, value, alpha = color.getHsl()
        if mode == 'cat':
            hue += NEXT_COLOR_HUE_DELTA_CAT
        elif mode == 'con':
            hue += NEXT_COLOR_HUE_DELTA_CON
        if sat == 0:
            sat = 255
            value = 128
            alpha = 255
            s = ""
        while hue > 360:
            hue -= 360
    
        return QColor.fromHsl(hue, sat, value, alpha)
    
    
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    def createQgsField(name : str, exampleValue, comment:str=None):
    
        """
        Create a QgsField using a Python-datatype exampleValue
        :param name: field name
        :param exampleValue: value, can be any type
        :param comment: (optional) field comment.
        :return: QgsField
        """
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        t = type(exampleValue)
        if t in [str]:
            return QgsField(name, QVariant.String, 'varchar', comment=comment)
        elif t in [bool]:
            return QgsField(name, QVariant.Bool, 'int', len=1, comment=comment)
        elif t in [int, np.int32, np.int64]:
            return QgsField(name, QVariant.Int, 'int', comment=comment)
    
        elif t in [float, np.double, np.float16, np.float32, np.float64]:
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            return QgsField(name, QVariant.Double, 'double', comment=comment)
        elif isinstance(exampleValue, np.ndarray):
            return QgsField(name, QVariant.String, 'varchar', comment=comment)
        elif isinstance(exampleValue, list):
    
            assert len(exampleValue) > 0, 'need at least one value in provided list'
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            v = exampleValue[0]
            prototype = createQgsField(name, v)
            subType = prototype.type()
            typeName = prototype.typeName()
            return QgsField(name, QVariant.List, typeName, comment=comment, subType=subType)
        else:
            raise NotImplemented()
    
    
    
    def setQgsFieldValue(feature:QgsFeature, field, value):
        """
        Wrties the Python value v into a QgsFeature field, taking care of required conversions
        :param feature: QgsFeature
        :param field: QgsField | field name (str) | field index (int)
        :param value: any python value
        """
    
        if isinstance(field, int):
            field = feature.fields().at(field)
        elif isinstance(field, str):
            field = feature.fields().at(feature.fieldNameIndex(field))
        assert isinstance(field, QgsField)
    
        if value is None:
            value = QVariant.NULL
        if field.type() == QVariant.String:
            value = str(value)
        elif field.type() in [QVariant.Int, QVariant.Bool]:
            value = int(value)
        elif field.type() in [QVariant.Double]:
            value = float(value)
        else:
            raise NotImplementedError()
    
       # i = feature.fieldNameIndex(field.name())
        feature.setAttribute(field.name(), value)
    
    
    
    def appendItemsToMenu(menu, itemsToAdd):
        """
        Appends items to QMenu "menu"
        :param menu: the QMenu to be extended
        :param itemsToAdd: QMenu or [list-of-QActions-or-QMenus]
        :return: menu
        """
        assert isinstance(menu, QMenu)
        if isinstance(itemsToAdd, QMenu):
            itemsToAdd = itemsToAdd.children()
        if not isinstance(itemsToAdd, list):
            itemsToAdd = [itemsToAdd]
        for item in itemsToAdd:
            if isinstance(item, QAction):
                #item.setParent(menu)
    
    
                a = menu.addAction(item.text(), item.triggered, item.shortcut())
                a.setEnabled(item.isEnabled())
                a.setIcon(item.icon())
                menu.addAction(a)
                s = ""
            elif isinstance(item, QMenu):
                item.setParent(menu)
                menu.addMenu(menu)
            else:
                s = ""
    
        return menu
    
    
    
    def toVectorLayer(src) -> QgsVectorLayer:
        """
        Returns a QgsRasterLayer if it can be extracted from src
        :param src: any type of input
        :return: QgsRasterLayer or None
        """
        lyr = None
        try:
            if isinstance(src, str):
                lyr = QgsVectorLayer(src)
            elif isinstance(src, QgsMimeDataUtils.Uri):
                lyr, b = src.vectorLayer()
                if not b:
                    lyr = None
            if isinstance(src, ogr.DataSource):
                path = src.GetDescription()
                bn = os.path.basename(path)
                lyr = QgsVectorLayer(path, bn, 'ogr')
            elif isinstance(src, QgsVectorLayer):
                lyr = src
    
        except Exception as ex:
            print(ex)
    
        return lyr
    
    def toDataset(src, readonly=True)->gdal.Dataset:
        """
        Returns a gdal.Dataset if it can be extracted from src
        :param src: input source
        :param readonly: bool, true by default, set False to upen the gdal.Dataset in update mode
        :return: gdal.Dataset
        """
        ga = gdal.GA_ReadOnly if readonly else gdal.GA_Update
        if isinstance(src, str):
            return gdal.Open(src, ga)
        elif isinstance(src, QgsRasterLayer) and src.dataProvider().name() == 'gdal':
            return toDataset(src.source(), readonly=readonly)
        elif isinstance(src, gdal.Dataset):
            return src
        else:
            return None
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    def toQgsMimeDataUtilsUri(mapLayer:QgsMapLayer)->QgsMimeDataUtils.Uri:
        """
        Creates a QgsMimeDataUtils.Uri that describes the QgsMapLayer `mapLayer`
        :param mapLayer: QgsMapLayer
        :return: QgsMimeDataUtils.Uri
        """
        assert isinstance(mapLayer, QgsMapLayer)
        uri = QgsMimeDataUtils.Uri()
        uri.uri = mapLayer.source()
        uri.name = mapLayer.name()
        if uri.name == '':
            uri.name = os.path.basename(uri.uri)
        uri.providerKey = mapLayer.dataProvider().name()
    
        if isinstance(mapLayer, QgsRasterLayer):
            uri.layerType = 'raster'
        elif isinstance(mapLayer, QgsVectorLayer):
            uri.layerType = 'vector'
        elif isinstance(mapLayer, QgsPluginLayer):
            uri.layerType = 'plugin'
        else:
            raise NotImplementedError()
        return uri
    
    
    
    def toMapLayer(src)->QgsMapLayer:
        """
        Return a QgsMapLayer if it can be extracted from src
        :param src: any type of input
        :return: QgsMapLayer
        """
        lyr = toRasterLayer(src)
        if isinstance(lyr, QgsMapLayer):
            return lyr
        lyr = toVectorLayer(src)
        if isinstance(lyr, QgsMapLayer):
            return lyr
        return lyr
    
    def toRasterLayer(src) -> QgsRasterLayer:
        """
        Returns a QgsRasterLayer if it can be extracted from src
        :param src: any type of input
        :return: QgsRasterLayer or None
        """
        lyr = None
        try:
            if isinstance(src, str):
                lyr = QgsRasterLayer(src)
            elif isinstance(src, QgsMimeDataUtils.Uri):
                lyr, b = src.rasterLayer('')
                if not b:
                    lyr = None
            elif isinstance(src, gdal.Dataset):
                lyr = QgsRasterLayer(src.GetFileList()[0], '', 'gdal')
            elif isinstance(src, QgsMapLayer) :
                lyr = src
            elif isinstance(src, gdal.Band):
                return toRasterLayer(src.GetDataset())
    
        except Exception as ex:
            print(ex)
    
        if isinstance(lyr, QgsRasterLayer) and lyr.isValid():
            return lyr
        else:
            return None
    
    
    
    def allSubclasses(cls):
        """
        Returns all subclasses of class 'cls'
        Thx to: http://stackoverflow.com/questions/3862310/how-can-i-find-all-subclasses-of-a-class-given-its-name
        :param cls:
        :return:
        """
        return cls.__subclasses__() + [g for s in cls.__subclasses__()
                                       for g in allSubclasses(s)]
    
    def scaledUnitString(num, infix=' ', suffix='B', div=1000):
    
        """
        Returns a human-readable file size string.
        thanks to Fred Cirera
        http://stackoverflow.com/questions/1094841/reusable-library-to-get-human-readable-version-of-file-size
        :param num: number in bytes
        :param suffix: 'B' for bytes by default.
        :param div: divisor of num, 1000 by default.
        :return: the file size string
        """
        for unit in ['','K','M','G','T','P','E','Z']:
            if abs(num) < div:
    
                return "{:3.1f}{}{}{}".format(num, infix, unit, suffix)
    
            num /= div
    
        return "{:.1f}{}{}{}".format(num, infix, unit, suffix)
    
    class SpatialPoint(QgsPointXY):
    
        """
        Object to keep QgsPoint and QgsCoordinateReferenceSystem together
        """
    
        @staticmethod
        def fromMapCanvasCenter(mapCanvas):
            assert isinstance(mapCanvas, QgsMapCanvas)
            crs = mapCanvas.mapSettings().destinationCrs()
            return SpatialPoint(crs, mapCanvas.center())
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        @staticmethod
        def fromMapLayerCenter(mapLayer:QgsMapLayer):
            assert isinstance(mapLayer, QgsMapLayer) and mapLayer.isValid()
            crs = mapLayer.crs()
            return SpatialPoint(crs, mapLayer.extent().center())
    
    
        @staticmethod
        def fromSpatialExtent(spatialExtent):
            assert isinstance(spatialExtent, SpatialExtent)
            crs = spatialExtent.crs()
            return SpatialPoint(crs, spatialExtent.center())
    
    
        def __init__(self, crs, *args):
    
            if not isinstance(crs, QgsCoordinateReferenceSystem):
                crs = QgsCoordinateReferenceSystem(crs)
    
            assert isinstance(crs, QgsCoordinateReferenceSystem)
            super(SpatialPoint, self).__init__(*args)
            self.mCrs = crs
    
    
        def __hash__(self):
            return hash(str(self))
    
    
        def setCrs(self, crs):
            assert isinstance(crs, QgsCoordinateReferenceSystem)
            self.mCrs = crs
    
        def crs(self):
            return self.mCrs
    
    
        def toPixelPosition(self, rasterDataSource, allowOutOfRaster=False):
            """
            Returns the pixel position of this SpatialPoint within the rasterDataSource
            :param rasterDataSource: gdal.Dataset
            :param allowOutOfRaster: set True to return out-of-raster pixel positions, e.g. QPoint(-1,0)
            :return: the pixel position as QPoint
            """
            ds = gdalDataset(rasterDataSource)
            ns, nl = ds.RasterXSize, ds.RasterYSize
            gt = ds.GetGeoTransform()
    
            pt = self.toCrs(ds.GetProjection())
            if pt is None:
                return None
    
            px = geo2px(pt, gt)
            if not allowOutOfRaster:
                if px.x() < 0 or px.x() >= ns:
                    return None
                if px.y() < 0 or px.y() >= nl:
                    return None
            return px
    
    
        def toCrs(self, crs):
            assert isinstance(crs, QgsCoordinateReferenceSystem)
    
            if self.mCrs != crs:
    
                pt = saveTransform(pt, self.mCrs, crs)
    
            return SpatialPoint(crs, pt) if pt else None
    
        def __reduce_ex__(self, protocol):
            return self.__class__, (self.crs().toWkt(), self.x(), self.y()), {}
    
        def __eq__(self, other):
            if not isinstance(other, SpatialPoint):
                return False
            return self.x() == other.x() and \
                   self.y() == other.y() and \
                   self.crs() == other.crs()
    
    
        def __copy__(self):
    
            return SpatialPoint(self.crs(), self.x(), self.y())
    
        def __str__(self):
            return self.__repr__()
    
    
        def __repr__(self):
    
    
            if self.crs().mapUnits() == QgsUnitTypes.DistanceDegrees:
                return '{:.1f} {:.1f} {}'.format(self.x(), self.y(), self.crs().authid())
            else:
                return '{:.5f} {:.5f} {:}'.format(self.x(), self.y(), self.crs().authid())
    
    
    
    
    def findParent(qObject, parentType, checkInstance = False):
        parent = qObject.parent()
        if checkInstance:
            while parent != None and not isinstance(parent, parentType):
                parent = parent.parent()
        else:
            while parent != None and type(parent) != parentType:
                parent = parent.parent()
        return parent
    
    
    def saveTransform(geom, crs1, crs2):
        assert isinstance(crs1, QgsCoordinateReferenceSystem)
        assert isinstance(crs2, QgsCoordinateReferenceSystem)
    
        result = None
        if isinstance(geom, QgsRectangle):
            if geom.isEmpty():
                return None
    
    
    
            transform = QgsCoordinateTransform()
            transform.setSourceCrs(crs1)
            transform.setDestinationCrs(crs2)
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
                rect = transform.transformBoundingBox(geom)
    
                result = SpatialExtent(crs2, rect)
            except:
    
                messageLog('Can not transform from {} to {} on rectangle {}'.format( \
    
                    crs1.description(), crs2.description(), str(geom)))
    
    
        elif isinstance(geom, QgsPointXY):
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            transform = QgsCoordinateTransform()
    
            transform.setSourceCrs(crs1)
            transform.setDestinationCrs(crs2)
    
                result = SpatialPoint(crs2, pt)
            except:
    
                messageLog('Can not transform from {} to {} on QgsPointXY {}'.format( \
    
                    crs1.description(), crs2.description(), str(geom)))
        return result
    
    
    
    def gdalDataset(pathOrDataset, eAccess=gdal.GA_ReadOnly):
        """
    
        :param pathOrDataset: path or gdal.Dataset
        :return: gdal.Dataset
        """
        if not isinstance(pathOrDataset, gdal.Dataset):
            pathOrDataset = gdal.Open(pathOrDataset, eAccess)
        assert isinstance(pathOrDataset, gdal.Dataset)
        return pathOrDataset
    
    
    def geo2pxF(geo, gt):
        """
        Returns the pixel position related to a Geo-Coordinate in floating point precision.
        :param geo: Geo-Coordinate as QgsPoint
        :param gt: GDAL Geo-Transformation tuple, as described in http://www.gdal.org/gdal_datamodel.html
        :return: pixel position as QPointF
        """
    
        assert isinstance(geo, QgsPointXY)
    
        # see http://www.gdal.org/gdal_datamodel.html
        px = (geo.x() - gt[0]) / gt[1]  # x pixel
        py = (geo.y() - gt[3]) / gt[5]  # y pixel
        return QPointF(px,py)
    
    
    def createGeoTransform(gsd, ul_x, ul_y):
        """
        Create a GDAL Affine GeoTransform vector for north-up images.
        See http://www.gdal.org/gdal_datamodel.html for details
        :param gsd: ground-sampling-distance / pixel-size
        :param ul_x: upper-left X
        :param ul_y: upper-left Y
        :return: (tuple)
        """
        if isinstance(gsd, tuple) or isinstance(gsd, list):
            gt1 = gsd[0] #pixel width
            gt5 = gsd[1] #pixel height
        else:
            gsd = float(gsd)
            gt1 = gt5 = gsd #pixel width == pixel height
    
        gt0 = ul_x
        gt3 = ul_y
    
        gt2 = gt4 = 0
    
        return (gt0, gt1, gt2, gt3, gt4, gt5)
    
    def geo2px(geo, gt):
        """
        Returns the pixel position related to a Geo-Coordinate as integer number.
        Floating-point coordinate are casted to integer coordinate, e.g. the pixel coordinate (0.815, 23.42) is returned as (0,23)
    
        :param geo: Geo-Coordinate as QgsPointXY
    
        :param gt: GDAL Geo-Transformation tuple, as described in http://www.gdal.org/gdal_datamodel.html
        :return: pixel position as QPpint
        """
        px = geo2pxF(geo, gt)
        return QPoint(int(px.x()), int(px.y()))
    
    
        """
        Converts a pixel coordinate into a geo-coordinate
    
        :param px: QPoint() with pixel coordinates
        :param gt: geo-transformation
        :param pxCenter: True to return geo-coordinate of pixel center, False to return upper-left edge
    
        #see http://www.gdal.org/gdal_datamodel.html
    
        gx = gt[0] + px.x()*gt[1]+px.y()*gt[2]
        gy = gt[3] + px.x()*gt[4]+px.y()*gt[5]
    
    
        if pxCenter:
            p2 = px2geo(QPoint(px.x()+1, px.y()+1), gt, pxCenter=False)
    
            gx = 0.5*(gx + p2.x())
            gy = 0.5*(gy + p2.y())
    
        return QgsPointXY(gx, gy)
    
    class SpatialExtent(QgsRectangle):
        """
        Object to keep QgsRectangle and QgsCoordinateReferenceSystem together
        """
        @staticmethod
        def fromMapCanvas(mapCanvas, fullExtent=False):
            assert isinstance(mapCanvas, QgsMapCanvas)
    
            if fullExtent:
                extent = mapCanvas.fullExtent()
            else:
                extent = mapCanvas.extent()
            crs = mapCanvas.mapSettings().destinationCrs()
            return SpatialExtent(crs, extent)
    
    
        @staticmethod
        def world():
            crs = QgsCoordinateReferenceSystem('EPSG:4326')
            ext = QgsRectangle(-180,-90,180,90)
            return SpatialExtent(crs, ext)
    
    
    
        @staticmethod
        def fromRasterSource(pathSrc):
            ds = gdalDataset(pathSrc)
            assert isinstance(ds, gdal.Dataset)
            ns, nl = ds.RasterXSize, ds.RasterYSize
            gt = ds.GetGeoTransform()
            crs = QgsCoordinateReferenceSystem(ds.GetProjection())
    
            xValues = []
            yValues = []
            for x in [0, ns]:
                for y in [0, nl]:
                    px = px2geo(QPoint(x,y), gt)
                    xValues.append(px.x())
                    yValues.append(px.y())
    
            return SpatialExtent(crs, min(xValues), min(yValues),
                                      max(xValues), max(yValues))
    
    
    
    
    
    
        @staticmethod
        def fromLayer(mapLayer):
            assert isinstance(mapLayer, QgsMapLayer)
            extent = mapLayer.extent()
            crs = mapLayer.crs()
            return SpatialExtent(crs, extent)
    
        def __init__(self, crs, *args):
    
            if not isinstance(crs, QgsCoordinateReferenceSystem):
                crs = QgsCoordinateReferenceSystem(crs)
    
            assert isinstance(crs, QgsCoordinateReferenceSystem)
            super(SpatialExtent, self).__init__(*args)
            self.mCrs = crs
    
        def setCrs(self, crs):
            assert isinstance(crs, QgsCoordinateReferenceSystem)
            self.mCrs = crs
    
        def crs(self):
            return self.mCrs
    
        def toCrs(self, crs):
            assert isinstance(crs, QgsCoordinateReferenceSystem)
            box = QgsRectangle(self)
            if self.mCrs != crs:
    
                box = saveTransform(box, self.mCrs, crs)
            return SpatialExtent(crs, box) if box else None
    
        def spatialCenter(self):
            return SpatialPoint(self.crs(), self.center())
    
    
        def combineExtentWith(self, *args):
            if args is None:
                return
            elif isinstance(args[0], SpatialExtent):
    
                ext = args[0]
                extent2 = ext.toCrs(self.crs())
    
                self.combineExtentWith(QgsRectangle(extent2))
            else:
                super(SpatialExtent, self).combineExtentWith(*args)
    
            return self
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        def setCenter(self, centerPoint:SpatialPoint, crs=None):
            """
            Sets the center. Can be used to move the SpatialExtent
            :param centerPoint:
            :param crs: QgsCoordinateReferenceSystem
            :return:
            """
            if isinstance(centerPoint, SpatialPoint):
                crs = centerPoint.crs()
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            if isinstance(crs, QgsCoordinateReferenceSystem) and crs != self.crs():
                trans = QgsCoordinateTransform()
                trans.setSourceCrs(crs)
                trans.setDestinationCrs(self.crs())
    
                centerPoint = trans.transform(centerPoint)
    
            delta = centerPoint - self.center()
            self.setXMaximum(self.xMaximum() + delta.x())
            self.setXMinimum(self.xMinimum() + delta.x())
            self.setYMaximum(self.yMaximum() + delta.y())
            self.setYMinimum(self.yMinimum() + delta.y())
    
            return self
    
        def __cmp__(self, other):
            if other is None: return 1
            s = ""
    
    
            return QgsPointXY(*self.upperRight())
    
            return QgsPointXY(*self.upperLeft())
    
            return QgsPointXY(*self.lowerRight())
    
            return QgsPointXY(*self.lowerLeft())
    
        def upperRight(self):
            return self.xMaximum(), self.yMaximum()
    
        def upperLeft(self):
            return self.xMinimum(), self.yMaximum()
    
        def lowerRight(self):
            return self.xMaximum(), self.yMinimum()
    
        def lowerLeft(self):
            return self.xMinimum(), self.yMinimum()
    
    
    
            if not isinstance(other, SpatialExtent):
                return False
    
            return self.toString() == other.toString()
    
        def __sub__(self, other):
            raise NotImplementedError()
    
        def __mul__(self, other):
            raise NotImplementedError()
    
        def __copy__(self):
            return SpatialExtent(self.crs(), QgsRectangle(self))
    
        def __reduce_ex__(self, protocol):
            return self.__class__, (self.crs().toWkt(),
                                    self.xMinimum(), self.yMinimum(),
                                    self.xMaximum(), self.yMaximum()
                                    ), {}
    
        def __str__(self):
            return self.__repr__()
    
    
        def __repr__(self):
    
            return '{} {} {}'.format(self.upperLeft(), self.lowerRight(), self.crs().authid())
    
    
        """
        Normalizes string, converts to lowercase, removes non-alpha characters,
        and converts spaces to hyphens.
        see https://stackoverflow.com/questions/295135/turn-a-string-into-a-valid-filename
        :return: path
        """
    
        text = re.sub(r'[^\w\s.-]', '', text).strip().lower()
        text = re.sub(r'[-\s]+', '_', text)
        return re.sub(r'[ ]+]','',text)
    
    def value2str(value, sep=' '):
        """
        Converts a value into a string
        :param value:
        :param sep:
        :return:
        """
        if isinstance(value, list):
            value = sep.join([str(v) for v in value])
        elif isinstance(value, np.array):
            value = value2str(value.astype(list), sep=sep)
        elif value is None:
            value = ''
        else:
            value = str(value)
        return value
    
    
    
    # works in Python 2 & 3
    class _Singleton(type):
        """ A metaclass that creates a Singleton base class when called. """
        _instances = {}
        def __call__(cls, *args, **kwargs):
            if cls not in cls._instances:
                cls._instances[cls] = super(_Singleton, cls).__call__(*args, **kwargs)
            return cls._instances[cls]
    
    class Singleton(_Singleton('SingletonMeta', (object,), {})): pass
    
    """
    #work, but require metaclass pattern 
    
    class Singleton(type):
        _instances = {}
    
        def __call__(cls, *args, **kwds):
            if cls not in cls._instances:
                cls._instances[cls] = super(Singleton, cls).__call__(*args,**kwds)
            return cls._instances[cls]
    
    class KeepRefs(object):
        __refs__ = defaultdict(list)
        def __init__(self):
            self.__refs__[self.__class__].append(weakref.ref(self))
    
        @classmethod
        def instances(cls):
            for inst_ref in cls.__refs__[cls]:
                inst = inst_ref()
                if inst is not None:
    
    def defaultBands(dataset):
        """
        Returns a list of 3 default bands
        :param dataset:
        :return:
        """
        if isinstance(dataset, str):
            return defaultBands(gdal.Open(dataset))
        elif isinstance(dataset, QgsRasterDataProvider):
            return defaultBands(dataset.dataSourceUri())
        elif isinstance(dataset, QgsRasterLayer):
            return defaultBands(dataset.source())
        elif isinstance(dataset, gdal.Dataset):
    
            db = dataset.GetMetadataItem(str('default_bands'), str('ENVI'))
            if db != None:
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
                db = [int(n) for n in re.findall(r'\d+')]
    
                return db
            db = [0, 0, 0]
            cis = [gdal.GCI_RedBand, gdal.GCI_GreenBand, gdal.GCI_BlueBand]
            for b in range(dataset.RasterCount):
                band = dataset.GetRasterBand(b + 1)
                assert isinstance(band, gdal.Band)
                ci = band.GetColorInterpretation()
                if ci in cis:
                    db[cis.index(ci)] = b
            if db != [0, 0, 0]:
                return db
    
            rl = QgsRasterLayer(dataset.GetFileList()[0])
            defaultRenderer = rl.renderer()
            if isinstance(defaultRenderer, QgsRasterRenderer):
                db = defaultRenderer.usesBands()
                if len(db) == 0:
                    return [0, 1, 2]
                if len(db) > 3:
                    db = db[0:3]
                db = [b-1 for b in db]
            return db
    
        else:
            raise Exception()
    
    ######### Lookup  tables
    METRIC_EXPONENTS = {
        "nm": -9, "um": -6, u"µm": -6, "mm": -3, "cm": -2, "dm": -1, "m": 0, "hm": 2, "km": 3
    }
    # add synonyms
    METRIC_EXPONENTS['nanometers'] = METRIC_EXPONENTS['nm']
    METRIC_EXPONENTS['micrometers'] = METRIC_EXPONENTS['um']
    METRIC_EXPONENTS['millimeters'] = METRIC_EXPONENTS['mm']
    METRIC_EXPONENTS['centimeters'] = METRIC_EXPONENTS['cm']
    METRIC_EXPONENTS['decimeters'] = METRIC_EXPONENTS['dm']
    METRIC_EXPONENTS['meters'] = METRIC_EXPONENTS['m']
    METRIC_EXPONENTS['hectometers'] = METRIC_EXPONENTS['hm']
    METRIC_EXPONENTS['kilometers'] = METRIC_EXPONENTS['km']
    
    
    LUT_WAVELENGTH = dict({'B': 480,
                           'G': 570,
                           'R': 660,
                           'NIR': 850,
                           'SWIR': 1650,
                           'SWIR1': 1650,
                           'SWIR2': 2150
                           })
    
    def convertMetricUnit(value, u1, u2):
        """converts value, given in unit u1, to u2"""
        assert u1 in METRIC_EXPONENTS.keys()
        assert u2 in METRIC_EXPONENTS.keys()
    
        e1 = METRIC_EXPONENTS[u1]
        e2 = METRIC_EXPONENTS[u2]
    
        return value * 10 ** (e1 - e2)
    
    def bandClosestToWavelength(dataset, wl, wl_unit='nm'):
        """
    
        Returns the band index (!) of an image dataset closest to wavelength `wl`.
    
        :param dataset: str | gdal.Dataset
        :param wl: wavelength to search the closed band for
        :param wl_unit: unit of wavelength. Default = nm
        :return: band index | 0 of wavelength information is not provided
        """
        if isinstance(wl, str):
            assert wl.upper() in LUT_WAVELENGTH.keys(), wl
            return bandClosestToWavelength(dataset, LUT_WAVELENGTH[wl.upper()], wl_unit='nm')
        else:
            try:
                wl = float(wl)
                ds_wl, ds_wlu = parseWavelength(dataset)
    
                if ds_wl is None or ds_wlu is None:
                    return 0
    
    
                if ds_wlu != wl_unit:
                    wl = convertMetricUnit(wl, wl_unit, ds_wlu)
                return int(np.argmin(np.abs(ds_wl - wl)))
            except:
                pass
        return 0
    
    
    def cloneRenderer(renderer):
    
        assert isinstance(renderer, QgsRasterRenderer)
        cloned = renderer.clone()
    
        #handle specific issues if cloning is not exactly the same
        if isinstance(cloned, QgsSingleBandPseudoColorRenderer):
            cloned.setClassificationMin(renderer.classificationMin())
            cloned.setClassificationMax(renderer.classificationMax())
    
        return cloned
    
    
    
    def parseWavelength(dataset):
        """
        Returns the wavelength + wavelength unit of a dataset
        :param dataset:
        :return: (wl, wl_u) or (None, None), if not existing
        """
    
        wl = None
        wlu = None
    
        if isinstance(dataset, str):
            return parseWavelength(gdal.Open(dataset))
        elif isinstance(dataset, QgsRasterDataProvider):
            return parseWavelength(dataset.dataSourceUri())