main.py

# -*- coding: utf-8 -*-
"""
/***************************************************************************
 EnMAPBox
                                 A QGIS plugin
 EnMAP-Box V3
                              -------------------
        begin                : 2015-08-20
        git sha              : $Format:%H$
        copyright            : (C) 2015 by HU-Berlin
        email                : bj@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.                                   *
 *                                                                         *
 ***************************************************************************/
"""

# Import the code for the dialog
import os, sys, re, fnmatch, collections, copy, traceback, six
from qgis.core import *
#os.environ['PATH'] += os.pathsep + r'C:\OSGeo4W64\bin'

from osgeo import gdal, ogr, osr, gdal_array

DEBUG = True
import qgis.analysis
try:
    from qgis.gui import *
    import qgis
    import qgis_add_ins
    qgis_available = True

    #import console.console_output
    #console.show_console()
    #sys.stdout = console.console_output.writeOut()
    #sys.stderr = console.console_output.writeOut()

except:
    print('Can not find QGIS instance')
    qgis_available = False

import numpy as np

import multiprocessing, site
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.uic.Compiler.qtproxies import QtGui, QtCore
import sys
import code
import codecs

#abbreviations
from timeseriesviewer import jp, mkdir, DIR_SITE_PACKAGES, file_search

site.addsitedir(DIR_SITE_PACKAGES)

#I don't know why, but this is required to run this in QGIS
#todo: still required?
path = os.path.abspath(jp(sys.exec_prefix, '../../bin/pythonw.exe'))
if os.path.exists(path):
    multiprocessing.set_executable(path)
    sys.argv = [ None ]

#ensure that required non-standard modules are available

import pyqtgraph as pg
from timeseriesviewer.ui import widgets

regLandsatSceneID = re.compile(r"L[EMCT][1234578]{1}[12]\d{12}[a-zA-Z]{3}\d{2}")


class TimeSeriesTableModel(QAbstractTableModel):
    columnames = ['date','sensor','ns','nl','nb','image','mask']

    def __init__(self, TS, parent=None, *args):
        super(QAbstractTableModel, self).__init__()
        assert isinstance(TS, TimeSeries)
        self.TS = TS

    def rowCount(self, parent = QModelIndex()):
        return len(self.TS)

    def columnCount(self, parent = QModelIndex()):
        return len(self.columnames)

    def removeRows(self, row, count , parent=QModelIndex()):
        self.beginRemoveRows(parent, row, row+count-1)
        toRemove = self._data[row:row+count]
        for i in toRemove:
            self._data.remove(i)

        self.endRemoveRows()

    def getDateFromIndex(self, index):
        if index.isValid():
            i = index.row()
            if i >= 0 and i < len(self.TS):
                return self.TS.getTSDs()[i]
        return None

    def getTimeSeriesDatumFromIndex(self, index):

        if index.isValid():
            i = index.row()
            if i >= 0 and i < len(self.TS):
                date = self.TS.getTSDs()[i]
                return self.TS.data[date]

        return None



    def data(self, index, role = Qt.DisplayRole):
        if role is None or not index.isValid():
            return None


        value = None
        ic_name = self.columnames[index.column()]
        TSD = self.getTimeSeriesDatumFromIndex(index)
        keys = list(TSD.__dict__.keys())
        if role == Qt.DisplayRole or role == Qt.ToolTipRole:
            if ic_name == 'name':
                value = os.path.basename(TSD.pathImg)
            elif ic_name == 'sensor':
                if role == Qt.ToolTipRole:
                    value = TSD.sensor.getDescription()
                else:
                    value = str(TSD.sensor)
            elif ic_name == 'date':
                value = '{}'.format(TSD.date)
            elif ic_name == 'image':
                value = TSD.pathImg
            elif ic_name == 'mask':
                value = TSD.pathMsk
            elif ic_name in keys:
                value = TSD.__dict__[ic_name]
            else:
                s = ""
        elif role == Qt.BackgroundColorRole:
            value = None
        elif role == Qt.UserRole:
            value = TSD

        return value

    #def flags(self, index):
    #    return Qt.ItemIsEnabled

    def flags(self, index):
        if index.isValid():
            item = self.getTimeSeriesDatumFromIndex(index)
            cname = self.columnames[index.column()]
            if cname.startswith('d'): #relative values can be edited
                flags = Qt.ItemIsEnabled | Qt.ItemIsSelectable | Qt.ItemIsEditable
            else:
                flags = Qt.ItemIsEnabled | Qt.ItemIsSelectable
            return flags
            #return item.qt_flags(index.column())
        return None

    def headerData(self, col, orientation, role):
        if Qt is None:
            return None
        if orientation == Qt.Horizontal and role == Qt.DisplayRole:
            return self.columnames[col]
        elif orientation == Qt.Vertical and role == Qt.DisplayRole:
            return col
        return None

class TimeSeriesItemModel(QAbstractItemModel):

    def __init__(self, TS):
        QAbstractItemModel.__init__(self)
        #self.rootItem = TreeItem[]
        assert type(TS) is TimeSeries
        self.TS = TS

    def index(self, row, column, parent = QModelIndex()):
        if not parent.isValid():
            parentItem = self.rootItem
        else:
            parentItem = parent.internalPointer()
        childItem = parentItem.child(row)
        if childItem:
            return self.createIndex(row, column, childItem)
        else:
            return QModelIndex()

    def setData(self, index, value, role = Qt.EditRole):
        if role == Qt.EditRole:
            row = index.row()

            return False
        return False

    def data(self, index, role=Qt.DisplayRole):
        data = None
        if role == Qt.DisplayRole or role == Qt.EditRole:
            data = 'sampletext'


        return data

    def flags(self, QModelIndex):
        return Qt.ItemIsSelectable

    def rowCount(self, index=QModelIndex()):
        return len(self.TS)

    #---------------------------------------------------------------------------
    def columnCount(self, index=QModelIndex()):
        return 1

LUT_SENSORNAMES = {(6, 30., 30.): 'Landsat Legacy' \
                  , (7,30.,30.): 'L8 OLI' \
                  , (4,10.,10.): 'S2 MSI 10m' \
                  , (6,20.,20.): 'S2 MSI 20m' \
                  , (3,30.,30.): 'S2 MSI 60m' \
                  , (3,30.,30.): 'S2 MSI 60m' \
                  , (5,5.,5.): 'RE 5m' \
                   }


class BandView(QObject):

    removeView = pyqtSignal(object)

    def __init__(self, TS, recommended_bands=None):
        assert type(TS) is TimeSeries
        self.representation = collections.OrderedDict()
        self.TS = TS
        self.TS.sensorAdded.connect(self.checkSensors)
        self.TS.changed.connect(self.checkSensors)

        self.Sensors = self.TS.Sensors

        import copy
        for sensor in self.Sensors:
            self.initSensor(copy.deepcopy(sensor), recommended_bands=recommended_bands)



    def checkSensors(self):
        represented_sensors = set(self.representation.keys())
        ts_sensors = set(self.TS.Sensors.keys())

        to_add = ts_sensors - represented_sensors
        to_remove = represented_sensors - ts_sensors
        for S in to_remove:
            self.representation[S].close()
            self.representation.pop(S)
        for S in to_add:
            self.initSensor(S)


    def initSensor(self, sensor, recommended_bands=None):
        """
        :param sensor:
        :param recommended_bands:
        :return:
        """

        assert type(sensor) is SensorInstrument
        if sensor not in self.representation.keys():
            #self.bandMappings[sensor] = ((0, 0, 5000), (1, 0, 5000), (2, 0, 5000))
            #x = imagechipviewsettings_widget.ImageChipViewSettings(sensor)
            #x = tsv_widgets.BandViewSettings(sensor)
            x = widgets.ImageChipViewSettings(sensor)
            #x.removeView.connect(lambda : self.removeView.emit(self))

            if recommended_bands is not None:
                assert min(recommended_bands) > 0
                if max(recommended_bands) < sensor.nb:
                    x.setBands(recommended_bands)
            x.create()
            self.representation[sensor] = x


    def getSensorStats(self, sensor, bands):
        """

        :param sensor:
        :param bands:
        :return:
        """
        assert type(sensor) is SensorInstrument
        dsRef = gdal.Open(self.Sensors[sensor][0])
        return [dsRef.GetRasterBand(b).ComputeRasterMinMax() for b in bands]


    def getRanges(self, sensor):
        return self.getWidget(sensor).getRGBSettings()[1]

    def getBands(self, sensor):
        return self.getWidget(sensor).getRGBSettings()[0]


    def getRGBSettings(self, sensor):
        return self.getWidget(sensor).getRGBSettings()

    def getWidget(self, sensor):
        assert type(sensor) is SensorInstrument
        return self.representation[sensor]



    def useMaskValues(self):

        #todo:
        return False


class SensorInstrument(object):

    INSTRUMENTS = dict()

    @staticmethod
    def fromRasterLayer(lyr):

        assert isinstance(lyr, QgsRasterLayer)
        nb = lyr.bandCount()
        sx = lyr.rasterUnitsPerPixelX()
        sy = lyr.rasterUnitsPerPixelY()

        bandNames = [lyr.bandName(i) for i in range(1, nb+1)]

        return SensorInstrument(nb, sx, sy,
                                band_names=bandNames,
                                wavelengths=None,
                                sensor_name=None)

    def fromGDALDataSet(self, ds):
        assert isinstance(ds, gdal.Dataset)
        nb = ds.RasterCount


    """
    Describes a Sensor Configuration
    """
    def __init__(self,nb, px_size_x,px_size_y, band_names=None, wavelengths=None, sensor_name=None):
        assert nb >= 1

        self.TS = None
        self.nb = nb
        self.px_size_x = float(abs(px_size_x))
        self.px_size_y = float(abs(px_size_y))

        assert self.px_size_x > 0
        assert self.px_size_y > 0

        if band_names is not None:
            assert len(band_names) == nb
        else:
            band_names = ['Band {}'.format(b+1) for b in range(nb)]

        self.band_names = band_names

        if wavelengths is not None:
            assert len(wavelengths) == nb

        self.wavelengths = wavelengths

        if sensor_name is None:
            id = (self.nb, self.px_size_x, self.px_size_y)
            if id in LUT_SENSORNAMES.keys():
                sensor_name = LUT_SENSORNAMES[id]
            else:
                sensor_name = '{}band@{}m'.format(self.nb, self.px_size_x)


        self.sensor_name = sensor_name

        self.hashvalue = hash(','.join(self.band_names))

    def __eq__(self, other):
        return self.nb == other.nb and \
               self.px_size_x == other.px_size_x and \
               self.px_size_y == other.px_size_y

    def __hash__(self):
        return self.hashvalue

    def __repr__(self):
        return self.sensor_name

    def getDescription(self):
        info = []
        info.append(self.sensor_name)
        info.append('{} Bands'.format(self.nb))
        info.append('Band\tName\tWavelength')
        for b in range(self.nb):
            if self.wavelengths:
                wl = str(self.wavelengths[b])
            else:
                wl = 'unknown'
            info.append('{}\t{}\t{}'.format(b+1, self.band_names[b],wl))

        return '\n'.join(info)



class ImageChipLabel(QLabel):

    clicked = pyqtSignal(object, object)


    def __init__(self, time_series_viewer=None, iface=None, TSD=None, bands=None):
        super(ImageChipLabel, self).__init__(time_series_viewer)
        self.TSV = time_series_viewer
        self.TSD = TSD
        self.bn = os.path.basename(self.TSD.pathImg)

        self.iface=iface
        self.bands=bands
        self.setContextMenuPolicy(Qt.DefaultContextMenu)
        self.setFrameShape(QFrame.StyledPanel)
        self.setSizePolicy(QSizePolicy.Fixed, QSizePolicy.Fixed)

        tt = ['Date: {}'.format(TSD.date) \
             ,'Name: {}'.format(self.bn) \
             ,'RGB:  {}'.format(','.join([str(b) for b in bands]))]

        self.setToolTip(list2str(tt))

    def mouseReleaseEvent(self, event):
	    self.clicked.emit(self, event)



    def contextMenuEvent(self, event):
        menu = QMenu()
        #add general options

        action = menu.addAction('Copy to clipboard')
        action.triggered.connect(lambda : QApplication.clipboard().setPixmap(self.pixmap()))


        #add QGIS specific options
        if self.iface:
            action = menu.addAction('Add {} to QGIS layers'.format(self.bn))
            action.triggered.connect(lambda : qgis_add_ins.add_QgsRasterLayer(self.iface, self.TSD.pathImg, self.bands))

        menu.exec_(event.globalPos())



class TimeSeries(QObject):
    datumAdded = pyqtSignal(name='datumAdded')

    #fire when a new sensor configuration is added
    sensorAdded = pyqtSignal(object, name='sensorAdded')


    changed = pyqtSignal()
    progress = pyqtSignal(int,int,int, name='progress')
    chipLoaded = pyqtSignal(object, name='chiploaded')
    closed = pyqtSignal()
    error = pyqtSignal(object)

    def __init__(self, imageFiles=None, maskFiles=None):
        QObject.__init__(self)

        #define signals

        #fire when a new TSD is added


        self.data = collections.OrderedDict()

        self.CHIP_BUFFER=dict()

        self.shape = None

        self.Sensors = collections.OrderedDict()

        self.Pool = None

        if imageFiles is not None:
            self.addFiles(imageFiles)
        if maskFiles is not None:
            self.addMasks(maskFiles)

    _sep = ';'


    def loadFromFile(self, path):

        images = []
        masks = []
        with open(path, 'r') as f:
            lines = f.readlines()
            for l in lines:
                if re.match('^[ ]*[;#&]', l):
                    continue

                parts = re.split('[\n'+TimeSeries._sep+']', l)
                parts = [p for p in parts if p != '']
                images.append(parts[0])
                if len(parts) > 1:
                    masks.append(parts[1])

        self.addFiles(images)
        self.addMasks(masks)


    def saveToFile(self, path):
        if path is None or len(path) == 0:
            return

        lines = []
        lines.append('#Time series definition file: {}'.format(np.datetime64('now').astype(str)))
        lines.append('#<image path>[;<mask path>]')
        for TSD in self.data.values():

            line = TSD.pathImg
            if TSD.pathMsk is not None:
                line += TimeSeries._sep + TSD.pathMsk

            lines.append(line)

        lines = [l+'\n' for l in lines]

        print('Write {}'.format(path))
        with open(path, 'w') as f:
            f.writelines(lines)


    def getSRS(self):
        if len(self.data) == 0:
            return 0
        else:
            TSD = self.data[self.getTSDs()[0]]
            return TSD.getSpatialReference()

    def getWKT(self):
        srs = self.getSRS()
        return srs.ExportToWkt()

    def getSceneCenter(self, srs=None):

        if srs is None:
            srs = self.getSRS()

        bbs = list()
        for S, TSDs in self.Sensors.items():
            x = []
            y = []
            for TSD in TSDs:
                bb = TSD.getBoundingBox(srs)
                x.extend([c[0] for c in bb])
                y.extend([c[1] for c in bb])

        return None
        pass

    def getMaxExtent(self, srs=None):

        x = []
        y = []

        if srs is None:
            srs = self.getSRS()

        for TSD in self.data.values():
            bb = TSD.getBoundingBox(srs)
            x.extend([c[0] for c in bb])
            y.extend([c[1] for c in bb])

        return (min(x), min(y), max(x), max(y))

    def getObservationDates(self):
        return [tsd.getDate() for tsd in self.data.keys()]

    def getTSDs(self, date_of_interest=None):
        if date_of_interest:
            tsds = [tsd for tsd in self.data.keys() if tsd.getDate() == date_of_interest]
        else:
            tsds = list(self.data.keys())
        return tsds

    def _callback_error(self, error):
        six.print_(error, file=sys.stderr)
        self.error.emit(error)
        self._callback_progress()

    def _callback_spatialchips(self, results):
        self.chipLoaded.emit(results)
        self._callback_progress()

    def _callback_progress(self):
        self._callback_progress_done += 1
        self.progress.emit(0, self._callback_progress_done, self._callback_progress_max)

        if self._callback_progress_done >= self._callback_progress_max:
            self._callback_progress_done = 0
            self._callback_progress_max = 0
            self.progress.emit(0,0,1)

    def getSpatialChips_parallel(self, bbWkt, srsWkt, TSD_band_list, ncpu=1):
        assert type(bbWkt) is str and type(srsWkt) is str

        import multiprocessing
        if self.Pool is not None:
            self.Pool.terminate()

        self.Pool = multiprocessing.Pool(processes=ncpu)


        self._callback_progress_max = len(TSD_band_list)
        self._callback_progress_done = 0

        for T in TSD_band_list:

            TSD = copy.deepcopy(T[0])
            bands = T[1]
            #TSD = pickle.dumps(self.data[date])
            args = (TSD, bbWkt, srsWkt)
            kwds = {'bands':bands}

            if six.PY3:
                self.Pool.apply_async(PFunc_TimeSeries_getSpatialChip, \
                                 args=args, kwds=kwds, \
                                 callback=self._callback_spatialchips, error_callback=self._callback_error)
            else:
                self.Pool.apply_async(PFunc_TimeSeries_getSpatialChip, \
                                      args, kwds, self._callback_spatialchips)

        s = ""

        pass

    def getImageChips(self, xy, size=50, bands=[4,5,6], dates=None):
        chipCollection = collections.OrderedDict()

        if dates is None:
            dates = self.data.keys()
        for date in dates:
            TSD = self.data[date]
            chipCollection[date] = TSD.readImageChip(xy, size=size, bands=bands)

        return chipCollection

    def addMasks(self, files, raise_errors=True, mask_value=0, exclude_mask_value=True):
        assert isinstance(files, list)
        l = len(files)

        self.progress.emit(0,0,l)
        for i, file in enumerate(files):

            try:
                self.addMask(file, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value, _quiet=True)
            except:
                pass

            self.progress.emit(0,i+1,l)

        self.progress.emit(0,0,1)
        self.changed.emit()

    def addMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True, _quiet=False):
        print('Add mask {}...'.format(pathMsk))
        ds = getDS(pathMsk)
        date = getImageDate(ds)

        if date in self.data.keys():
            TSD = self.data[date]

            if not _quiet:
                self.changed.emit()

            return TSD.setMask(pathMsk, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value)
        else:
            info = 'TimeSeries does not contain date {} {}'.format(date, pathMsk)
            if raise_errors:
                raise Exception(info)
            else:
                six.print_(info, file=sys.stderr)
            return False

    def removeAll(self):
        self.clear()

    def clear(self):
        self.Sensors.clear()
        self.data.clear()
        self.changed.emit()


    def removeDates(self, TSDs):
        for TSD in TSDs:
            self.removeTSD(TSD, _quiet=True)
        self.changed.emit()

    def removeTSD(self, TSD, _quiet=False):

        assert type(TSD) is TimeSeriesDatum
        S = TSD.sensor
        self.Sensors[S].remove(TSD)
        self.data.pop(TSD, None)
        if len(self.Sensors[S]) == 0:
            self.Sensors.pop(S)
        if not _quiet:
            self.changed.emit()


    def addFile(self, pathImg, pathMsk=None, _quiet=False):

        six.print_(pathImg)
        six.print_('Add image {}...'.format(pathImg))

        try:
            sensorAdded = False
            TSD = TimeSeriesDatum(pathImg, pathMsk=pathMsk)
            existingSensors = list(self.Sensors.keys())
            if TSD.sensor not in existingSensors:
                self.Sensors[TSD.sensor] = list()
                sensorAdded = True
            else:
                TSD.sensor = existingSensors[existingSensors.index(TSD.sensor)]

            if TSD in self.data.keys():
                six.print_('Time series datum already added: {}'.format(str(TSD)), file=sys.stderr)
            else:
                self.Sensors[TSD.sensor].append(TSD)
                self.data[TSD] = TSD

            if sensorAdded:
                self.sensorAdded.emit(TSD.sensor)

            if not _quiet:
                self._sortTimeSeriesData()
                self.changed.emit()
                self.datumAdded.emit()
        except:
            exc_type, exc_value, exc_traceback = sys.exc_info()
            traceback.print_exception(exc_type, exc_value, exc_traceback, limit=2)
            six.print_('Unable to add {}'.format(file), file=sys.stderr)
            pass



    def addFiles(self, files):
        assert isinstance(files, list)
        l = len(files)
        assert l > 0

        self.progress.emit(0,0,l)
        for i, file in enumerate(files):
            self.addFile(file, _quiet=True)
            self.progress.emit(0,i+1,l)

        self._sortTimeSeriesData()
        self.progress.emit(0,0,1)
        self.datumAdded.emit()
        self.changed.emit()



    def _sortTimeSeriesData(self):
        self.data = collections.OrderedDict(sorted(self.data.items(), key=lambda t:t[0]))

    def __len__(self):
        return len(self.data)

    def __repr__(self):
        info = []
        info.append('TimeSeries:')
        l = len(self)
        info.append('  Scenes: {}'.format(l))

        if l > 0:
            keys = list(self.data.keys())
            info.append('  Range: {} to {}'.format(keys[0], keys[-1]))
        return '\n'.join(info)

def PFunc_TimeSeries_getSpatialChip(TSD, bbWkt, srsWkt , bands=[4,5,3]):

    chipdata = TSD.readSpatialChip(bbWkt, srs=srsWkt, bands=bands)

    return TSD, chipdata

def px2Coordinate(gt, pxX, pxY, upper_left=True):
    cx = gt[0] + pxX*gt[1] + pxY*gt[2]
    cy = gt[3] + pxX*gt[4] + pxY*gt[5]
    if not upper_left:
        cx += gt[1]*0.5
        cy += gt[5]*0.5
    return cx, cy

def coordinate2px(gt, cx, cy):
    px = int((cx - gt[0]) / gt[1])
    py = int((cy - gt[3]) / gt[5])
    return px, py


def getBoundingBoxPolygon(points, srs=None):
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for point in points:
        ring.AddPoint(point[0], point[1])
    bb = ogr.Geometry(ogr.wkbPolygon)
    bb.AddGeometry(ring)
    if srs:
        bb.AssignSpatialReference(srs)
    return bb



def transformGeometry(geom, crsSrc, crsDst, trans=None):
    if trans is None:
        assert isinstance(crsSrc, QgsCoordinateReferenceSystem)
        assert isinstance(crsDst, QgsCoordinateReferenceSystem)
        return transformGeometry(geom, None, None, trans=QgsCoordinateTransform(crsSrc, crsDst))
    else:
        assert isinstance(trans, QgsCoordinateTransform)
        return trans.transform(geom)




def getDS(ds):
    if type(ds) is not gdal.Dataset:
        ds = gdal.Open(ds)
    return ds

regAcqDate = re.compile(r'acquisition[ ]*(time|date|day)', re.I)
def getImageDate2(lyr):
    assert isinstance(lyr, QgsRasterLayer)
    mdLines = str(lyr.metadata()).splitlines()
    date = None
    #find date in metadata
    for line in [l for l in mdLines if regAcqDate.search(l)]:
        date = parseAcquisitionDate(line)
        if date:
            return date
    #find date in filename
    dn, fn = os.path.split(str(lyr.dataProvider().dataSourceUri()))
    date = parseAcquisitionDate(fn)
    if date: return date

    #find date in file directory path
    date = parseAcquisitionDate(date)

    return date

def getBandNames(lyr):
    assert isinstance(lyr, QgsRasterLayer)
    dp = lyr.dataProvider()
    assert isinstance(dp, QgsRasterDataProvider)
    if str(dp.name()) == 'gdal':
        s = ""
    else:
        return lyr

def getImageDate(ds):
    if type(ds) is str:
        ds = gdal.Open(ds)

    path = ds.GetFileList()[0]
    to_check = [os.path.basename(path), os.path.dirname(path)]

    regAcqDate = re.compile(r'acquisition (time|date|day)', re.I)
    for key, value in ds.GetMetadata_Dict().items():
        if regAcqDate.search(key):
            to_check.insert(0, value)

    for text in to_check:
        date = parseAcquisitionDate(text)
        if date:
            return date

    raise Exception('Can not identify acquisition date of {}'.format(path))


class TimeSeriesDatum(object):
    """
    Collects all data sets related to one sensor
    """

    def __init__(self, pathImg, pathMsk=None):
        self.pathImg = pathImg
        self.pathMsk = None

        self.lyrImg = QgsRasterLayer(pathImg, os.path.basename(pathImg), False)
        self.crs = self.lyrImg.dataProvider().crs()
        self.sensor = SensorInstrument.fromRasterLayer(self.lyrImg)

        self.date = getImageDate2(self.lyrImg)
        assert self.date is not None, 'Unable to find acquisition date of {}'.format(pathImg)

        self.ns = self.lyrImg.width()
        self.nl = self.lyrImg.height()
        self.nb = self.lyrImg.bandCount()
        self.srs_wkt = str(self.crs.toWkt())


        if pathMsk:
            self.setMask(pathMsk)

    def getdtype(self):
        return gdal_array.GDALTypeCodeToNumericTypeCode(self.etype)

    def getDate(self):
        return np.datetime64(self.date)


    def getSpatialReference(self):
        return self.crs
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.srs_wkt)
        return srs


    def getBoundingBox(self, srs=None):

        bbox = self.lyrImg.extent()
        if srs:
            assert isinstance(srs, QgsCoordinateReferenceSystem)
            bbox = transformGeometry(bbox, self.crs, srs)
        return bbox


    def setMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True):
        dsMsk = gdal.Open(pathMsk)
        mskDate = getImageDate(dsMsk)


        errors = list()
        if mskDate and mskDate != self.getDate():
            errors.append('Mask date differs from image date')
        if self.ns != dsMsk.RasterXSize or self.nl != dsMsk.RasterYSize:
            errors.append('Spatial size differs')
        if dsMsk.RasterCount != 1:
            errors.append('Mask has > 1 bands')

        projImg = self.getSpatialReference()
        projMsk = osr.SpatialReference()
        projMsk.ImportFromWkt(dsMsk.GetProjection())

        if not projImg.IsSame(projMsk):
            errors.append('Spatial Reference differs from image')
        if self.gt != list(dsMsk.GetGeoTransform()):
            errors.append('Geotransformation differs from image')

        if len(errors) > 0:
            errors.insert(0, 'pathImg:{} \npathMsk:{}'.format(self.pathImg, pathMsk))
            errors = '\n'.join(errors)
            if raise_errors:
                raise Exception(errors)
            else:
                six.print_(errors, file=sys.stderr)
                return False
        else:
            self.pathMsk = pathMsk
            self.mask_value = mask_value
            self.exclude_mask_value = exclude_mask_value

            return True

    def readSpatialChip(self, geometry, srs=None, bands=[4,5,3]):

        srs_img = osr.SpatialReference()
        srs_img.ImportFromWkt(self.srs_wkt)


        if type(geometry) is ogr.Geometry:
            g_bb = geometry
            srs_bb = g_bb.GetSpatialReference()
        else:
            assert srs is not None and type(srs) in [str, osr.SpatialReference]
            if type(srs) is str:
                srs_bb = osr.SpatialReference()
                srs_bb.ImportFromWkt(srs)
            else:
                srs_bb = srs.Clone()
            g_bb = ogr.CreateGeometryFromWkt(geometry, srs_bb)

        assert srs_bb is not None and g_bb is not None
        assert g_bb.GetGeometryName() == 'POLYGON'

        if not srs_img.IsSame(srs_bb):
            g_bb = g_bb.Clone()
            g_bb.TransformTo(srs_img)

        cx0,cx1,cy0,cy1 = g_bb.GetEnvelope()

        ul_px = coordinate2px(self.gt, min([cx0, cx1]), max([cy0, cy1]))
        lr_px = coordinate2px(self.gt, max([cx0, cx1]), min([cy0, cy1]))
        lr_px = [c+1 for c in lr_px] #+1

        return self.readImageChip([ul_px[0], lr_px[0]], [ul_px[1], lr_px[1]], bands=bands)

    def readImageChip(self, px_x, px_y, bands=[4,5,3]):