Skip to content
Snippets Groups Projects
sensecarbon_tsv.py 51.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • unknown's avatar
    unknown committed
    # -*- 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.                                   *
     *                                                                         *
     ***************************************************************************/
    """
    
    unknown's avatar
    unknown committed
    
    
    # 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
    
    
    
    
    unknown's avatar
    unknown committed
    try:
        from qgis.gui import *
    
    unknown's avatar
    unknown committed
        import qgis
    
    unknown's avatar
    unknown committed
        import qgis_add_ins
    
    unknown's avatar
    unknown committed
        qgis_available = True
    except:
        qgis_available = False
    
    import numpy as np
    import pickle
    
    unknown's avatar
    unknown committed
    import six
    import multiprocessing
    
    
    
    #I don't know why but this is required to run this in QGIS
    
    
    unknown's avatar
    unknown committed
    path = os.path.abspath(os.path.join(sys.exec_prefix, '../../bin/pythonw.exe'))
    if os.path.exists(path):
        multiprocessing.set_executable(path)
        sys.argv = [ None ]
    
    unknown's avatar
    unknown committed
    
    
    pluginDir = os.path.dirname(__file__)
    sys.path.append(pluginDir)
    sys.path.append(os.path.join(pluginDir, 'qimage2ndarray'))
    
    import qimage2ndarray
    
    
    unknown's avatar
    unknown committed
    from PyQt4.QtCore import *
    from PyQt4.QtGui import *
    from sensecarbon_tsv_gui import SenseCarbon_TSVGui
    
    
    unknown's avatar
    unknown committed
    DEBUG = True
    
    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','name','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()
    
    
    unknown's avatar
    unknown committed
        def getDateFromIndex(self, index):
            if index.isValid():
                i = index.row()
                if i >= 0 and i < len(self.TS):
                    return self.TS.getDates()[i]
            return None
    
    unknown's avatar
    unknown committed
    
        def getTimeSeriesDatumFromIndex(self, index):
    
            if index.isValid():
                i = index.row()
                if i >= 0 and i < len(self.TS):
                    date = self.TS.getDates()[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 == '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:
    
    unknown's avatar
    unknown committed
    
            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
    
    
    unknown's avatar
    unknown committed
    class TimeSeriesItemModel(QAbstractItemModel):
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
        def __init__(self, TS):
    
    unknown's avatar
    unknown committed
            QAbstractItemModel.__init__(self)
            #self.rootItem = TreeItem[]
    
    unknown's avatar
    unknown committed
            assert type(TS) is TimeSeries
            self.TS = TS
    
    unknown's avatar
    unknown committed
    
        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()):
    
    unknown's avatar
    unknown committed
            return len(self.TS)
    
    unknown's avatar
    unknown committed
    
        #---------------------------------------------------------------------------
        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' \
                      }
    
    unknown's avatar
    unknown committed
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    class SensorConfiguration(object):
    
        """
        Describes a Sensor Configuration
        """
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        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
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            if band_names is not None:
                assert len(band_names) == nb
    
            else:
                band_names = ['Band {}'.format(b+1) for b in range(nb)]
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            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)
    
            print(sensor_name)
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            self.sensor_name = sensor_name
    
    
            self.hashvalue = hash(','.join(self.band_names))
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        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
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        def __repr__(self):
    
            return self.sensor_name
    
    unknown's avatar
    unknown committed
    class TimeSeries(QObject):
    
        #define signals
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    
        #fire when a new TSD is added
        datumAdded = pyqtSignal(name='datumAdded')
    
        #fire when a new sensor configuration is added
        sensorAdded = pyqtSignal(object, name='sensorAdded')
    
    
    
    unknown's avatar
    unknown committed
        changed = pyqtSignal()
    
    unknown's avatar
    unknown committed
        progress = pyqtSignal(int,int,int, name='progress')
        chipLoaded = pyqtSignal(object, name='chiploaded')
        closed = pyqtSignal()
        error = pyqtSignal(object)
    
        data = collections.OrderedDict()
    
        CHIP_BUFFER=dict()
    
        shape = None
    
    
        SensorConfigurations = collections.OrderedDict()
    
    unknown's avatar
    unknown committed
        def __init__(self, imageFiles=None, maskFiles=None):
            QObject.__init__(self)
    
            self.Pool = None
    
            if imageFiles is not None:
                self.addFiles(imageFiles)
            if maskFiles is not None:
                self.addMasks(maskFiles)
    
    
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
        def getMaxExtent(self, srs=None):
    
    unknown's avatar
    unknown committed
    
            x = []
            y = []
    
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            if srs is None:
                TSD = self.data[self.getDates()[0]]
                srs = TSD.getSpatialReference()
    
    
    unknown's avatar
    unknown committed
            for TSD in self.data.values():
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
                bb = TSD.getBoundingBox(srs)
    
    unknown's avatar
    unknown committed
                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 getDates(self):
            return list(self.data.keys())
    
    
    
        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, dates=None, bands=[4,5,3], 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)
    
            if dates is None:
                dates = self.getDates()
    
            self._callback_progress_max = len(dates)
            self._callback_progress_done = 0
    
            for date in dates:
    
                TSD = copy.deepcopy(self.data[date])
                #TSD = pickle.dumps(self.data[date])
                args = (TSD, bbWkt, srsWkt)
                kwds = {'bands':bands}
    
    
    unknown's avatar
    unknown committed
                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)
    
    unknown's avatar
    unknown committed
    
            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)
            assert l > 0
    
            self.progress.emit(0,0,l)
            for i, file in enumerate(files):
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    
                try:
                    self.addMask(file, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value, _quiet=True)
                except:
                    pass
    
    
    unknown's avatar
    unknown committed
                self.progress.emit(0,i+1,l)
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
            self.progress.emit(0,0,l)
    
    unknown's avatar
    unknown committed
            self.changed()
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
        def addMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True, _quiet=False):
    
    unknown's avatar
    unknown committed
            print('Add mask {}...'.format(pathMsk))
            ds = getDS(pathMsk)
            date = getImageDate(ds)
    
            if date in self.data.keys():
                TSD = self.data[date]
    
    unknown's avatar
    unknown committed
    
                if not _quiet:
                    self.changed.emit()
    
    
    unknown's avatar
    unknown committed
                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
    
    
    unknown's avatar
    unknown committed
        def removeAll(self):
            self.clear()
    
        def clear(self):
            dates = list(self.data.keys())
            for date in dates:
                self.removeDate(date, _quiet=True)
            self.changed.emit()
    
    
        def removeDates(self, dates):
            for date in dates:
                self.removeDate(date, _quiet=True)
            self.changed.emit()
    
        def removeDate(self, date, _quiet=False):
    
    unknown's avatar
    unknown committed
            assert type(date) is np.datetime64
    
    unknown's avatar
    unknown committed
            self.data.pop(date, None)
            if len(self.data) == 0:
                self.nb = None
                self.bandnames = None
                self.srs = None
    
    unknown's avatar
    unknown committed
            if not _quiet:
                self.changed.emit()
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
    
        def addFile(self, pathImg, pathMsk=None, _quiet=False):
    
    unknown's avatar
    unknown committed
            print(pathImg)
    
    unknown's avatar
    unknown committed
            print('Add image {}...'.format(pathImg))
    
    
            try:
                TSD = TimeSeriesDatum(pathImg, pathMsk=pathMsk)
    
    unknown's avatar
    unknown committed
    
    
                sensor = SensorConfiguration(TSD.nb, TSD.gt[1], TSD.gt[5], TSD.bandnames, TSD.wavelength)
                if sensor not in self.SensorConfigurations.keys():
                    self.SensorConfigurations[sensor] = list()
                self.SensorConfigurations[sensor].append(TSD)
    
    unknown's avatar
    unknown committed
    
    
                self.data[TSD.getDate()] = TSD
    
                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
    
    unknown's avatar
    unknown committed
    
    
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
        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)
    
    unknown's avatar
    unknown committed
                self.progress.emit(0,i+1,l)
    
            self._sortTimeSeriesData()
            self.progress.emit(0,0,l)
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
            self.datumAdded.emit()
    
    unknown's avatar
    unknown committed
    
    
    
        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]))
                info.append('  Bands: {} [{}]'.format(self.nb,','.join(self.bandnames)))
            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.getDate(), 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 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())
    
    
    unknown's avatar
    unknown committed
            refBand = dsImg.GetRasterBand(1)
            self.etype = refBand.DataType
    
            self.nodata = refBand.GetNoDataValue()
    
    
            self.bandnames = list()
    
    unknown's avatar
    unknown committed
            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()
            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
    
    unknown's avatar
    unknown committed
    
            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:
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
                assert type(srs) is osr.SpatialReference
    
    unknown's avatar
    unknown committed
                my_srs = self.getSpatialReference()
                if not my_srs.IsSame(srs):
                    #todo: consider srs differences
    
    Benjamin Jakimow's avatar
    Benjamin Jakimow committed
    
                    raise Exception('differences SRS in bounding box request')
    
    unknown's avatar
    unknown committed
                    pass
    
            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 = geometry
                srs_bb = g.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 = ogr.CreateGeometryFromWkt(geometry, srs_bb)
    
            assert srs_bb is not None and g is not None
            assert g.GetGeometryName() == 'POLYGON'
    
            if not srs_img.IsSame(srs_bb):
                g = g.Clone()
                g.TransformTo(srs_img)
            cx0,cx1,cy0,cy1 = g.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)
    
    unknown's avatar
    unknown committed
    
            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
    
    
            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
    
    
    unknown's avatar
    unknown committed
            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', 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)
                    chipdata[b] = data
                    nodatavalue = band.GetNoDataValue()
                    if nodatavalue is not None:
                        templateMsk[j0:j1,i0:i1] = np.logical_and(templateMsk[j0:j1,i0:i1], data[j0:j1,i0:i1] != nodatavalue)
    
                if self.pathMsk:
                    ds = gdal.Open(self.pathMsk)
                    tmp = ds.GetRasterBand(1).ReadAsArray(xoff=x0, yoff=y0, \
                                win_xsize=win_xsize,win_ysize=win_ysize) == 0
    
                    templateMsk[j0:j1,i0:i1] = np.logical_and(templateMsk[j0:j1,i0:i1], tmp)
    
            chipdata['mask'] = templateMsk
            return chipdata
    
        def __repr__(self):
            info = []
            info.append('TS Datum {}'.format(self.date))
            for key in ['pathImg', 'pathMsk']:
                info.append('  {}:{}'.format(key,self.__dict__[key]))
            return '\n'.join(info)
    
    regYYYYDOY = re.compile(r'(19|20)\d{5}')
    regYYYYMMDD = re.compile(r'(19|20)\d{2}-\d{2}-\d{2}')
    
    unknown's avatar
    unknown committed
    regYYYY = re.compile(r'(19|20)\d{2}')
    
    unknown's avatar
    unknown committed
    def parseAcquisitionDate(text):
        match = regYYYYMMDD.search(text)
        if match:
            return np.datetime64(match.group())
    
    unknown's avatar
    unknown committed
        match = regYYYYDOY.search(text)
    
    unknown's avatar
    unknown committed
        if match:
            return getDateTime64FromYYYYDOY(match.group())
    
    unknown's avatar
    unknown committed
        match = regYYYY.search(text)
        if match:
            return np.datetime64(match.group())
    
    unknown's avatar
    unknown committed
        return None
    
    
    
    unknown's avatar
    unknown committed
    def getDateTime64FromYYYYDOY(yyyydoy):
        return getDateTime64FromDOY(yyyydoy[0:4], yyyydoy[4:7])
    
    def getDateTime64FromDOY(year, doy):
            if type(year) is str:
                year = int(year)
            if type(doy) is str:
                doy = int(doy)
            return np.datetime64('{:04d}-01-01'.format(year)) + np.timedelta64(doy-1, 'D')
    
    
    class PictureTest(QMainWindow):
    
    
        def __init__(self, parent=None, qImage=None):
    
    unknown's avatar
    unknown committed
            super(PictureTest,self).__init__(parent)
            self.setWindowTitle("Show Image with pyqt")
            self.imageLabel=QLabel()
            self.imageLabel.setSizePolicy(QSizePolicy.Ignored,QSizePolicy.Ignored)
            self.setCentralWidget(self.imageLabel)
    
            self.cv_img = None
    
    
            if qImage:
                self.addImage(qImage)
    
    unknown's avatar
    unknown committed
    
    
        def addImage(self, qImage):
            pxmap = QPixmap.fromImage(qImage)
            self.addPixmap(pxmap)
    
        def addPixmap(self, pixmap):
            pxmap = pixmap.scaled(self.imageLabel.size(), Qt.KeepAspectRatio)
            self.imageLabel.setPixmap(pxmap)
            self.imageLabel.adjustSize()
            self.imageLabel.update()
    
        def addNumpy(self, data):
    
    unknown's avatar
    unknown committed
    
    
            img = Array2Image(data)
    
    unknown's avatar
    unknown committed
    
            #self.resize(img.width(), img.height())
    
    def getChip3d(chips, rgb_idx, ranges):
        assert len(rgb_idx) == 3 and len(rgb_idx) == len(ranges)
        for i in rgb_idx:
            assert i in chips.keys()
    
        nl, ns = chips[rgb_idx[0]].shape
        a3d = np.ndarray((3,nl,ns), dtype='float')
    
        for i, rgb_i in enumerate(rgb_idx):
            range = ranges[i]
            data = chips[rgb_i].astype('float')
            data -= range[0]
            data *= 255./range[1]
            a3d[i,:] = data
    
        np.clip(a3d, 0, 255, out=a3d)
    
        return a3d.astype('uint8')
    
    def Array2Image(d3d):
        nb, nl, ns = d3d.shape
        byteperline = nb
        d3d = d3d.transpose([1,2,0]).copy()
    
        return QImage(d3d.data, ns, nl, QImage.Format_RGB888)
    
    
    class VerticalLabel(QLabel):
        def __init__(self, text):
            super(self.__class__, self).__init__()
            self.text = text
    
        def paintEvent(self, event):
            painter = QPainter(self)
            painter.setPen(Qt.black)
            painter.translate(20, 100)
            painter.rotate(-90)
            if self.text:
                painter.drawText(0, 0, self.text)
            painter.end()
    
    unknown's avatar
    unknown committed
    
    
        def minimumSizeHint(self):
            size = QLabel.minimumSizeHint(self)
            return QSize(size.height(), size.width())
    
        def sizeHint(self):
            size = QLabel.sizeHint(self)
            return QSize(size.height(), size.width())
    
    unknown's avatar
    unknown committed
    
    
    unknown's avatar
    unknown committed
    class ImageChipBuffer(object):
    
        data = dict()
        BBox = None
        SRS = None
    
        def __init__(self):
    
            pass
    
        def hasDataCube(self, date):
            return date in self.data.keys()
    
        def getMissingBands(self, date, bands):
            missing = set(bands)
            if date in self.data.keys():
    
    unknown's avatar
    unknown committed
                missing = missing - set(self.data[date].keys())
    
    unknown's avatar
    unknown committed
    
            return missing
    
        def addDataCube(self, date, chipData):
            assert self.BBox is not None, 'Please initialize the bounding box first.'
    
            assert type(date) == np.datetime64
            assert isinstance(chipData, dict)
            if date not in self.data.keys():
                self.data[date] = dict()
            self.data[date].update(chipData)
    
        def getDataCube(self, date):
            return self.data.get(date)
    
    
        def getChipRGB(self, date, view):
    
    unknown's avatar
    unknown committed
            bands = view.getBands()
            band_ranges = view.getRanges()
            assert len(bands) == 3 and len(bands) == len(band_ranges)
            assert date in self.data, 'Date {} is not in buffer'.format(date)
            chipData = self.data[date]
            for b in bands:
                assert b in chipData.keys()
    
            nl, ns = chipData[bands[0]].shape
            rgb_data = np.ndarray((3,nl,ns), dtype='float')
    
            for i, b in enumerate(bands):
                range = band_ranges[i]
                data = chipData[b].astype('float')
                data -= range[0]
                data *= 255./range[1]
                rgb_data[i,:] = data
    
            np.clip(rgb_data, 0, 255, out=rgb_data)
            rgb_data = rgb_data.astype('uint8')
    
            if view.useMaskValues():
                rgb = view.getMaskColor()
                is_masked = np.where(np.logical_not(chipData['mask']))
                for i, c in enumerate(rgb):
                    rgb_data[i, is_masked[0], is_masked[1]] = c
    
    
    unknown's avatar
    unknown committed
    
    
        def getChipImage(self, date, view):
            rgb = self.getChipRGB(date, view)
            nb, nl, ns = rgb.shape
            rgb = rgb.transpose([1,2,0]).copy('C')
            return QImage(rgb.data, ns, nl, QImage.Format_RGB888)
    
    unknown's avatar
    unknown committed
    
        def clear(self):
            self.data.clear()
    
        def setBoundingBox(self, BBox):
            assert type(BBox) is ogr.Geometry
            SRS = BBox.GetSpatialReference()
            assert SRS is not None
            if self.BBox is None or not self.BBox.Equals(BBox) or not self.SRS.IsSame(SRS):
                self.data.clear()
                self.BBox = BBox
                self.SRS = SRS
    
        def __repr__(self):
            info = ['Chipbuffer']
            info.append('Bounding Box: {}'.format(self.bbBoxWkt))
            info.append('Chips: {}'.format(len(self.data)))
            return '\n'.join(info)
    
    
    class SenseCarbon_TSV:
        """QGIS Plugin Implementation."""