sensecarbon_tsv.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
from qgis.core import *
#os.environ['PATH'] += os.pathsep + r'C:\OSGeo4W64\bin'

from osgeo import gdal, ogr, osr, gdal_array

DEBUG = True

try:
    from qgis.gui import *
    import qgis
    import qgis_add_ins
    qgis_available = True
except:
    qgis_available = False

import numpy as np
import six
import multiprocessing
from PyQt4.QtCore import *
from PyQt4.QtGui import *


#abbreviations
mkdir = lambda p: os.makedirs(p, exist_ok=True)
jp = os.path.join

def expand_python_path(path):
    if path not in sys.path:
        sys.path.append(path)



#I don't know why, but this is required to run this in QGIS
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
PLUGIN_DIR = os.path.dirname(__file__)
LIB_DIR = jp(PLUGIN_DIR, 'libs')
expand_python_path(PLUGIN_DIR)
expand_python_path(LIB_DIR)
try:
    import pyqtgraph
except:
    expand_python_path(jp(LIB_DIR,'pyqtgraph'))

import pyqtgraph as pg

import tsv_widgets
from sensecarbon_tsv_gui import *



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

def file_search(rootdir, wildcard, recursive=False, ignoreCase=False):
    assert rootdir is not None
    if not os.path.isdir(rootdir):
        six.print_("Path is not a directory:{}".format(rootdir), file=sys.stderr)

    results = []

    for root, dirs, files in os.walk(rootdir):
        for file in files:
            if (ignoreCase and fnmatch.fnmatch(file.lower(), wildcard.lower())) \
                    or fnmatch.fnmatch(file, wildcard):
                results.append(os.path.join(root, file))
        if not recursive:
            break
    return results



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 Qt is None or index.isValid() == False:
            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.): 'L7 ETM+' \
                  ,(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(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.initSensor)

        self.Sensors = self.TS.Sensors

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





    def initSensor(self, sensor, recommended_bands=None):
        assert type(sensor) is SensorConfiguration
        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 = tsv_widgets.ImageChipViewSettings(sensor)

            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):
        assert type(sensor) is SensorConfiguration
        dsRef = gdal.Open(self.Sensors[sensor][0])
        return [dsRef.GetRasterBand(b).ComputeRasterMinMax() for b in bands]


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

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

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


    def useMaskValues(self):

        #todo:
        return False


class SensorConfiguration(object):
    """
    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 = '{} b x {} 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):
    def __init__(self, parent=None, iface=None, TSD=None, bands=None):
        super(ImageChipLabel, self).__init__(parent)
        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 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):
        import time
        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 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 getDS(ds):
    if type(ds) is not gdal.Dataset:
        ds = gdal.Open(ds)
    return ds

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):

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

        dsImg = gdal.Open(pathImg)
        assert dsImg

        date = getImageDate(dsImg)
        assert date is not None
        self.date = date.astype(str)

        self.ns = dsImg.RasterXSize
        self.nl = dsImg.RasterYSize
        self.nb = dsImg.RasterCount

        self.srs_wkt = dsImg.GetProjection()
        self.gt = list(dsImg.GetGeoTransform())

        refBand = dsImg.GetRasterBand(1)
        self.etype = refBand.DataType

        self.nodata = refBand.GetNoDataValue()

        self.bandnames = list()
        for b in range(self.nb):
            name = dsImg.GetRasterBand(b+1).GetDescription()
            if name is None or name == '':
                name = 'Band {}'.format(b+1)
            self.bandnames.append(name)

        self.wavelength = None
        domains = dsImg.GetMetadataDomainList()
        if domains:
            for domain in domains:
                md = dsImg.GetMetadata_Dict(domain)
                if 'wavelength' in md.keys():
                    wl = md['wavelength']
                    wl = re.split('[;,{}]', wl)
                    wl = [float(w) for w in wl]
                    assert len(wl) == self.nb
                    self.wavelength = wl
                    break

        self.sensor = SensorConfiguration(self.nb, self.gt[1], self.gt[5], self.bandnames, self.wavelength)


        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):
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.srs_wkt)
        return srs

    def getBoundingBox(self, srs=None):
        ext = list()


        for px in [0,self.ns]:
            for py in [0, self.nl]:
                ext.append(px2Coordinate(self.gt, px, py))



        if srs is not None:
            assert type(srs) is osr.SpatialReference
            my_srs = self.getSpatialReference()
            if not my_srs.IsSame(srs):
                #todo: consider srs differences
                trans = osr.CoordinateTransformation(my_srs, srs)
                ext = trans.TransformPoints(ext)
                ext = [(e[0], e[1]) for e in ext]

        return ext


    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]):

        ds = gdal.Open(self.pathImg, gdal.GA_ReadOnly)

        assert len(px_x) == 2 and px_x[0] <= px_x[1]
        assert len(px_y) == 2 and px_y[0] <= px_y[1]

        ns = px_x[1]-px_x[0]+1
        nl = px_y[1]-px_y[0]+1
        assert ns >= 0
        assert nl >= 0

        src_ns = ds.RasterXSize
        src_nl = ds.RasterYSize


        chipdata = dict()



        #pixel indices in source image
        x0 = max([0, px_x[0]])
        y0 = max([0, px_y[0]])
        x1 = min([src_ns, px_x[1]])
        y1 = min([src_nl, px_y[1]])
        win_xsize = x1-x0+1
        win_ysize = y1-y0+1

        #pixel indices in image chip (ideally 0 and ns-1 or nl-1)
        i0 = x0 - px_x[0]
        i1 = i0 + win_xsize

        j0 = y0 - px_y[0]
        j1 = j0+ win_ysize




        templateImg = np.zeros((nl,ns))
        if self.nodata:
            templateImg *= self.nodata

        templateImg = templateImg.astype(self.getdtype())
        templateMsk = np.ones((nl,ns), dtype='bool')

        if win_xsize < 1 or win_ysize < 1:
            six.print_('Selected image chip is out of raster image {}'.format(self.pathImg), file=sys.stderr)
            for i, b in enumerate(bands):
                chipdata[b] = np.copy(templateImg)

        else:
            for i, b in enumerate(bands):
                band = ds.GetRasterBand(b)
                data = np.copy(templateImg)
                data[j0:j1,i0:i1] = band.ReadAsArray(xoff=x0, yoff=y0, win_xsize=win_xsize,win_ysize=win_ysize)