Skip to content
Snippets Groups Projects
sensecarbon_tsv.py 40.5 KiB
Newer Older
unknown's avatar
unknown committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
# -*- 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.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon
# Initialize Qt resources from file resources.py
import resources
try:
    from qgis.core import *
    from qgis.gui import *
    qgis_available = True
except:
    qgis_available = False

# Import the code for the dialog
from sensecarbon_tsv_gui import SenseCarbon_TSVGui
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import skimage
from osgeo import gdal, ogr, osr, gdal_array
import numpy as np
import pickle
import six
import os, sys, re, fnmatch, collections, copy
import multiprocessing
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()


    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:
            value = self._data[index.row()]

        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 TSItemModel(QAbstractItemModel):

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

    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)

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


class TimeSeries(QObject):

    #define signals
    datumAdded = pyqtSignal(list, name='datumAdded')
    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

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

        self.Pool = None
        self.nb = None
        self.bandnames = list()


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



    def getMaxExtent(self):

        x = []
        y = []

        for TSD in self.data.values():
            bb = TSD.getBoundingBox(self.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 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}

            self.Pool.apply_async(PFunc_TimeSeries_getSpatialChip, \
                             args=args, kwds=kwds, \
                             callback=self._callback_spatialchips, error_callback=self._callback_error)

        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):
            self.addMask(file, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value)
            self.progress.emit(0,i+1,l)
        self.progress.emit(0,0,l)

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

        if date in self.data.keys():
            TSD = self.data[date]
            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 removeDate(self, date):
        self.data.pop(date, None)
        if len(self.data) == 0:
            self.nb = None
            self.bandnames = None
            self.srs = None

    def addFile(self, pathImg, pathMsk=None, _quiet=True):
        print('Add image {}...'.format(pathImg))
        TSD = TimeSeriesDatum(pathImg, pathMsk=pathMsk)


        if self.nb is None:
            self.nb = TSD.nb
            self.bandnames = TSD.bandnames
            self.srs = TSD.getSpatialReference()

        else:
            assert self.nb == TSD.nb

        self.data[TSD.getDate()] = TSD

        if not _quiet:
            self._sortTimeSeriesData()
            self.datumAdded.emit(self.bandnames[:])

    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,l)
        self.datumAdded.emit(self.bandnames[:])



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

        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)

        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:
            my_srs = self.getSpatialReference()
            if not my_srs.IsSame(srs):
                #todo: consider srs differences
                raise NotImplemented()
                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)

        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.ones((nl,ns)) * 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', 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}')
def parseAcquisitionDate(text):
    match = regYYYYMMDD.search(text)
    if match:
        return np.datetime64(match.group())
    match = regYYYYDOY
    if match:
        return getDateTime64FromYYYYDOY(match.group())
    return None


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

    def addNumpy(self, data):


        if False:
            nb, nl, ns = data.shape
            byteperline = nb
            data = data.transpose([1,2,0]).copy()

            img = QImage(data.data, ns, nl, QImage.Format_RGB888)
            self.imageLabel.setPixmap(QPixmap.fromImage(img))
        img = Array2Image(data)
        pxMap = QPixmap.fromImage(img)
        pxMap = pxMap.scaled(self.imageLabel.size(), Qt.KeepAspectRatio)
        self.imageLabel.setPixmap(pxMap)

        #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 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():
            missing = missing - self.data[date].keys()

        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 getChipImage(self, date, view):
        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

        rgb_data = rgb_data.transpose([1,2,0]).copy()
        return QImage(rgb_data.data, ns, nl, QImage.Format_RGB888)


    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."""

    def __init__(self, iface):
        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # Save reference to the QGIS interface
        self.iface = iface
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = 'placeholder'
        #locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'EnMAPBox_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)

            if qVersion() > '4.3.3':
                QCoreApplication.installTranslator(self.translator)

        # Create the dialog (after translation) and keep reference
        self.dlg = SenseCarbon_TSVGui()
        D = self.dlg
        self.TS = TimeSeries()
        self.TS.datumAdded.connect(self.ua_datumAdded)
        self.TS.progress.connect(self.ua_TSprogress)
        self.TS.chipLoaded.connect(self.ua_showPxCoordinate_addChips)
        D.tableView_TimeSeries.setModel(TimeSeriesTableModel(self.TS))
        D.tableView_TimeSeries.horizontalHeader().setResizeMode(QHeaderView.ResizeToContents)

        self.VIEWS = list()
        self.ImageChipBuffer = ImageChipBuffer()
        self.CHIPWIDGETS = collections.OrderedDict()

        self.ValidatorPxX = QIntValidator(0,99999)
        self.ValidatorPxY = QIntValidator(0,99999)
        D.btn_showPxCoordinate.clicked.connect(lambda: self.ua_showPxCoordinate_start())
        D.btn_selectByCoordinate.clicked.connect(self.ua_selectByCoordinate)
        D.btn_selectByCoordinate.clicked.connect(self.ua_selectByRectangle)
        D.btn_addBandView.clicked.connect(lambda :self.ua_addView())
        D.btn_addTSImages.clicked.connect(lambda :self.ua_addTSImages())
        D.btn_addTSMasks.clicked.connect(lambda :self.ua_addTSMasks())
        D.btn_removeTSD.clicked.connect(lambda : self.ua_removeTSD(None))

        D.spinBox_ncpu.setRange(0, multiprocessing.cpu_count())

        # Declare instance attributes
        self.actions = []
        #self.menu = self.tr(u'&EnMAP-Box')
        # TODO: We are going to let the user set this up in a future iteration
        #self.toolbar = self.iface.addToolBar(u'EnMAPBox')
        #self.toolbar.setObjectName(u'EnMAPBox')

        self.CPV = self.dlg.scrollAreaWidgetContents.layout()
        self.check_enabled()
        s = ""

    def ua_selectByRectangle(self):

        pass

    def ua_selectByCoordinate(self):

        pass

    def ua_TSprogress(self, v_min, v, v_max):
        assert v_min <= v and v <= v_max
        P = self.dlg.progressBar
        if P.minimum() != v_min or P.maximum() != v_max:
            P.setRange(v_min, v_max)
        P.setValue(v)

    def ua_datumAdded(self):
        cb_centerdate = self.dlg.cb_centerdate
        cb_centerdate.clear()
        if len(self.TS) > 0:
            if self.dlg.spinBox_coordinate_x.value() == 0.0 and \
               self.dlg.spinBox_coordinate_y.value() == 0.0:
                xmin, ymin, xmax, ymax = self.TS.getMaxExtent()
                self.dlg.spinBox_coordinate_x.setRange(xmin, xmax)
                self.dlg.spinBox_coordinate_y.setRange(ymin, ymax)
                self.dlg.spinBox_coordinate_x.setValue(0.5*(xmin+xmax))
                self.dlg.spinBox_coordinate_y.setValue(0.5*(ymin+ymax))
                s =""

                for date in self.TS.getDates():
                    cb_centerdate.addItem(date.astype('str'), date)
            s = ""
        self.dlg.tableView_TimeSeries.resizeColumnsToContents()

    def check_enabled(self):
        D = self.dlg
        D.btn_showPxCoordinate.setEnabled(len(self.VIEWS) > 0 and len(self.TS) > 0)
        D.btn_selectByCoordinate.setEnabled(qgis_available)
        D.btn_selectByRectangle.setEnabled(qgis_available)


    # noinspection PyMethodMayBeStatic
    def tr(self, message):
        """Get the translation for a string using Qt translation API.

        We implement this ourselves since we do not inherit QObject.

        :param message: String for translation.
        :type message: str, QString

        :returns: Translated version of message.
        :rtype: QString
        """
        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('EnMAPBox', message)


    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=True,
        status_tip=None,
        whats_this=None,
        parent=None):
        """Add a toolbar icon to the toolbar.

        :param icon_path: Path to the icon for this action. Can be a resource
            path (e.g. ':/plugins/foo/bar.png') or a normal file system path.
        :type icon_path: str

        :param text: Text that should be shown in menu items for this action.
        :type text: str

        :param callback: Function to be called when the action is triggered.
        :type callback: function

        :param enabled_flag: A flag indicating if the action should be enabled
            by default. Defaults to True.
        :type enabled_flag: bool

        :param add_to_menu: Flag indicating whether the action should also
            be added to the menu. Defaults to True.
        :type add_to_menu: bool

        :param add_to_toolbar: Flag indicating whether the action should also
            be added to the toolbar. Defaults to True.
        :type add_to_toolbar: bool

        :param status_tip: Optional text to show in a popup when mouse pointer
            hovers over the action.
        :type status_tip: str

        :param parent: Parent widget for the new action. Defaults None.
        :type parent: QWidget

        :param whats_this: Optional text to show in the status bar when the
            mouse pointer hovers over the action.

        :returns: The action that was created. Note that the action is also
            added to self.actions list.
        :rtype: QAction
        """

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            self.toolbar.addAction(action)

        if add_to_menu:
            self.iface.addPluginToMenu(
                self.menu,
                action)

        self.actions.append(action)