timeseries.py

from __future__ import absolute_import
import six, sys, os, gc, re, collections, site, inspect, time, traceback, copy

import bisect, datetime
from osgeo import gdal, ogr

from qgis import *
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from PyQt4.QtXml import *

import numpy as np

from timeseriesviewer import DIR_REPO, DIR_EXAMPLES, dprint, jp, findAbsolutePath

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

METRIC_EXPONENTS = {
    "nm":-9,"um": -6, "mm":-3, "cm":-2, "dm":-1, "m": 0,"hm":2, "km":3
}
#add synonyms
METRIC_EXPONENTS['nanometers'] = METRIC_EXPONENTS['nm']
METRIC_EXPONENTS['micrometers'] = METRIC_EXPONENTS['um']
METRIC_EXPONENTS['millimeters'] = METRIC_EXPONENTS['mm']
METRIC_EXPONENTS['centimeters'] = METRIC_EXPONENTS['cm']
METRIC_EXPONENTS['decimeters'] = METRIC_EXPONENTS['dm']
METRIC_EXPONENTS['meters'] = METRIC_EXPONENTS['m']
METRIC_EXPONENTS['hectometers'] = METRIC_EXPONENTS['hm']
METRIC_EXPONENTS['kilometers'] = METRIC_EXPONENTS['km']

def convertMetricUnit(value, u1, u2):
    assert u1 in METRIC_EXPONENTS.keys()
    assert u2 in METRIC_EXPONENTS.keys()

    e1 = METRIC_EXPONENTS[u1]
    e2 = METRIC_EXPONENTS[u2]

    return value * 10**(e1-e2)

def verifyVRT(pathVRT):

    ds = gdal.Open(pathVRT)
    if ds is None:
        return False
    s = ""

    return True



class SensorInstrument(QObject):

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


    LUT_Wavelengths = dict({'B':480,
                            'G':570,
                            'R':660,
                            'nIR':850,
                            'swIR':1650,
                            'swIR1':1650,
                            'swIR2':2150
                            })
    """
    Describes a Sensor Configuration
    """
    def __init__(self, refLyr, sensor_name=None):
        super(SensorInstrument, self).__init__()
        assert isinstance(refLyr, QgsRasterLayer)
        assert refLyr.isValid()
        #QgsMapLayerRegistry.instance().addMapLayer(refLyr)
        self.nb = refLyr.bandCount()
        self.bandDataType = refLyr.dataProvider().dataType(1)
        self.refUri = refLyr.dataProvider().dataSourceUri()
        r = refLyr.renderer()


        """
        dom = QDomDocument()
        root = dom.createElement('root')
        refLyr.renderer().writeXML(dom, root)
        dom.appendChild(root)
        self.renderXML = dom.toString()
        """

        #todo: better band names
        self.bandNames = [refLyr.bandName(i) for i in range(1, self.nb + 1)]
        #self.refLyr = refLyr

        self.TS = None
        px_size_x = refLyr.rasterUnitsPerPixelX()
        px_size_y = refLyr.rasterUnitsPerPixelY()

        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

        #find wavelength
        wl, wlu = parseWavelength(refLyr)
        self.wavelengths = np.asarray(wl)
        self.wavelengthUnits = wlu

        if sensor_name is None:
            id = (self.nb, self.px_size_x, self.px_size_y)
            sensor_name = '{}bands@{}m'.format(self.nb, self.px_size_x)

        self.sensorName = sensor_name

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

    def dataType(self, p_int):
        return self.bandDataType

    def bandClosestToWavelength(self, wl, wl_unit='nm'):
        """
        Returns the band index (>=0) of the band closest to wavelength wl
        :param wl:
        :param wl_unit:
        :return:
        """
        if not self.wavelengthsDefined():
            return None

        if wl in SensorInstrument.LUT_Wavelengths.keys():
            wl_unit = 'nm'
            wl = SensorInstrument.LUT_Wavelengths[wl]

        wl = float(wl)
        if self.wavelengthUnits != wl_unit:
            wl = convertMetricUnit(wl, wl_unit, self.wavelengthUnits)


        return np.argmin(np.abs(self.wavelengths - wl))




    def wavelengthsDefined(self):
        return self.wavelengths is not None and \
                self.wavelengthUnits is not None

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

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

        return '\n'.join(info)


class TSDLoaderSignals(QObject):

    sigRasterLayerLoaded = pyqtSignal(QgsRasterLayer)
    sigFinished = pyqtSignal()

class TSDLoader(QRunnable):
    """
    Runnable to load QgsRasterLayers from a parallel thread
    """
    def __init__(self, tsd_paths):
        super(TSDLoader, self).__init__()
        self.signals = TSDLoaderSignals()

        self.paths = list()
        for p in tsd_paths:
            if not (isinstance(p, tuple) or isinstance(p, list)):
                p = [p]
            self.paths.append(p)


    def run(self):
        lyrs = []
        for path in self.paths:
            TSD = TimeSeriesDatum(path)
            lyr = QgsRasterLayer(path)
            if lyr:
                lyrs.append(lyr)
                self.signals.sigRasterLayerLoaded.emit(lyr)
                dprint('{} loaded'.format(path))
            else:
                dprint('Failed to load {}'.format(path))
        self.signals.sigFinished.emit()
        #return lyrs


def verifyInputImage(path, vrtInspection=''):

    if not os.path.exists(path):
        print('{}Image does not exist: '.format(vrtInspection, path))
        return False
    ds = gdal.Open(path)
    if not ds:
        print('{}GDAL unable to open: '.format(vrtInspection, path))
        return False
    if ds.GetDriver().ShortName == 'VRT':
        vrtInspection = 'Inspection {}\n'.format(path)
        validSrc = [verifyInputImage(p, vrtInspection=vrtInspection) for p in set(ds.GetFileList()) - set([path])]
        return all(validSrc)
    else:
        return True


class TimeSeriesDatum(QObject):
    @staticmethod
    def createFromPath(path):
        """
        Creates a valid TSD or returns None if this is impossible
        :param path:
        :return:
        """
        p = findAbsolutePath(path)
        if verifyInputImage(p):
            return TimeSeriesDatum(None, p)
        else:
            return None



    """
    Collects all data sets related to one sensor
    """
    sigVisibilityChanged = pyqtSignal(bool)
    sigRemoveMe = pyqtSignal()



    def __init__(self, timeSeries, pathImg, pathMsk=None):
        super(TimeSeriesDatum,self).__init__()
        assert os.path.exists(pathImg)
        self.timeSeries = timeSeries
        self.pathImg = pathImg
        self.pathMsk = None

        assert os.path.exists(pathImg)
        self.lyrImg = QgsRasterLayer(pathImg, os.path.basename(pathImg), False)
        assert self.lyrImg.isValid()
        self.uriImg = self.lyrImg.dataProvider().dataSourceUri()

        self.crs = self.lyrImg.dataProvider().crs()
        self.sensor = SensorInstrument(self.lyrImg)

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

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

        self.mVisibility = True

        if pathMsk:
            self.setMask(pathMsk)


    def rank(self):
        return self.timeSeries.index(self)

    def setVisibility(self, b):
        old = self.mVisibility
        self.mVisibility = b
        if old != self.mVisibility:
            self.sigVisibilityChanged.emit(b)

    def isVisible(self):
        return self.mVisibility


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


    def getSpatialReference(self):
        return self.crs

    def spatialExtent(self):
        from timeseriesviewer.main import SpatialExtent
        extent = SpatialExtent(self.lyrImg.crs(), self.lyrImg.extent())
        return extent

    def __repr__(self):
        return 'TS Datum {} {}'.format(self.date, str(self.sensor))

    def __cmp__(self, other):
        return cmp(str((self.date, self.sensor)), str((other.date, other.sensor)))

    #def __eq__(self, other):
    #    return self.date == other.date and self.sensor == other.sensor

    def __lt__(self, other):
        if self.date < other.date:
            return True
        else:
            return self.sensor.sensorName < other.sensor.sensorName

    def __hash__(self):
        return hash((self.date,self.sensor.sensorName))


class TimeSeries(QObject):

    sigTimeSeriesDatesAdded = pyqtSignal(list)
    sigTimeSeriesDatesRemoved = pyqtSignal(list)

    sigSensorAdded = pyqtSignal(SensorInstrument)
    sigSensorRemoved = pyqtSignal(SensorInstrument)

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

        #define signals

        #fire when a new TSD is added


        #self.data = collections.OrderedDict()
        self.data = list()

        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, n_max=None):

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

        if n_max:
            n_max = min([len(images), n_max])
            self.addFiles(images[0:n_max])
        else:
            self.addFiles(images)
        #self.addMasks(masks)


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

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

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

            lines.append(line)

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

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

    def getMaxSpatialExtent(self, crs=None):
        if len(self.data) == 0:
            return None

        extent = self.data[0].spatialExtent()
        if len(self.data) > 1:
            for TSD in self.data[1:]:
                extent.combineExtentWith(TSD.spatialExtent())

        return extent


    def tsdFromPath(self, path):
        for tsd in self.data:
            if tsd.pathImg == path:
                return tsd
        return False

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

    def getTSD(self, pathOfInterest):
        for tsd in self.data:
            if tsd.pathImg == pathOfInterest:
                return tsd
        return None

    def getTSDs(self, dateOfInterest=None, sensorOfInterest=None):
        tsds = self.data[:]
        if dateOfInterest:
            tsds = [tsd for tsd in tsds if tsd.getDate() == dateOfInterest]
        if sensorOfInterest:
            tsds = [tsd for tsd in tsds if tsd.sensor == sensorOfInterest]
        return tsds

    def clear(self):
        self.Sensors.clear()
        dates = self.data[:]
        del self.data[:]
        self.sigTimeSeriesDatesRemoved.emit(dates)


    def removeDates(self, TSDs):
        removed = list()
        for TSD in TSDs:
            assert type(TSD) is TimeSeriesDatum
            self.data.remove(TSD)
            TSD.timeSeries = None
            removed.append(TSD)

            S = TSD.sensor
            #self.Sensors[S].remove(TSD)
            if len(self.Sensors[S]) == 0:
                self.Sensors.pop(S)
                self.sigSensorRemoved(S)

        self.sigTimeSeriesDatesRemoved.emit(removed)


    def addTimeSeriesDates(self, timeSeriesDates):
        assert isinstance(timeSeriesDates, list)
        added = list()
        for TSD in timeSeriesDates:
            try:
                sensorAdded = False
                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:
                    six.print_('Time series datum already added: {}'.format(str(TSD)), file=sys.stderr)
                else:
                    self.Sensors[TSD.sensor].append(TSD)
                    #insert sorted

                    bisect.insort(self.data, TSD)
                    TSD.timeSeries = self
                    TSD.sigRemoveMe.connect(lambda : self.removeDates([TSD]))
                    added.append(TSD)
                if sensorAdded:
                    self.sigSensorAdded.emit(TSD.sensor)

            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

        if len(added) > 0:
            self.sigTimeSeriesDatesAdded.emit(added)


    def addFiles(self, files):
        assert isinstance(files, list)
        for i, file in enumerate(files):
            tsd = TimeSeriesDatum.createFromPath(file)
            if tsd is None:
                dprint('Unable to add: {}'.format(file), file=sys.stderr)
            else:
                self.addTimeSeriesDates([tsd])

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

    def __iter__(self):
        return iter(self.data)

    def __getitem__(self, slice):
        return self.data[slice]

    def __delitem__(self, slice):
        self.removeDates(slice)

    def __contains__(self, item):
        return item in self.data

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


        return '\n'.join(info)

regAcqDate = re.compile(r'acquisition[ ]*(time|date|day)', re.I)
regLandsatSceneID = re.compile(r"L[EMCT][1234578]{1}[12]\d{12}[a-zA-Z]{3}\d{2}")

def getImageDate(lyr):
    assert isinstance(lyr, QgsRasterLayer)
    mdLines = str(lyr.metadata()).splitlines()
    date = None
    #find date in metadata
    for line in [l for l in mdLines if regAcqDate.search(l)]:
        date = parseAcquisitionDate(line)
        if date:
            return date
    #find date in filename
    dn, fn = os.path.split(str(lyr.dataProvider().dataSourceUri()))
    date = parseAcquisitionDate(fn)
    if date: return date

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

    return date


def 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


regYYYYDOY = re.compile(r'(19|20)\d{5}')
regYYYYMMDD = re.compile(r'(19|20)\d{2}-\d{2}-\d{2}')
regYYYY = re.compile(r'(19|20)\d{2}')




def parseWavelength(lyr):
    wl = None
    wlu = None
    assert isinstance(lyr, QgsRasterLayer)
    md = [l.split('=') for l in str(lyr.metadata()).splitlines() if 'wavelength' in l.lower()]
    #see http://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html for supported wavelength units
    regWLU = re.compile('((micro|nano|centi)meters)|(um|nm|mm|cm|m|GHz|MHz)')
    for kv in md:
        key, value = kv
        key = key.lower()
        if key == 'center wavelength':
            tmp = re.findall('\d*\.\d+|\d+', value) #find floats
            if len(tmp) == 0:
                tmp = re.findall('\d+', value) #find integers
            if len(tmp) == lyr.bandCount():
                wl = [float(w) for w in tmp]

        if key == 'wavelength units':
            match = regWLU.search(value)
            if match:
                wlu = match.group()

            names = ['nanometers','micrometers','millimeters','centimeters','decimenters']
            si   = ['nm','um','mm','cm','dm']
            if wlu in names:
                wlu = si[names.index(wlu)]

    return wl, wlu



def parseAcquisitionDate(text):
    match = regLandsatSceneID.search(text)
    if match:
        id = match.group()
        return getDateTime64FromYYYYDOY(id[9:16])
    match = regYYYYMMDD.search(text)
    if match:
        return np.datetime64(match.group())
    match = regYYYYDOY.search(text)
    if match:
        return getDateTime64FromYYYYDOY(match.group())
    match = regYYYY.search(text)
    if match:
        return np.datetime64(match.group())
    return None



def getDateTime64FromYYYYDOY(yyyydoy):
    return getDateFromDOY(yyyydoy[0:4], yyyydoy[4:7])

def getDOYfromDate(dt):

    return (dt.astype('datetime64[D]') - dt.astype('datetime64[Y]')).astype(int)+1

def getDateFromDOY(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')


if __name__ == '__main__':

    assert getDOYfromDate(np.datetime64('2014-01-01')) == 1
    assert getDOYfromDate(np.datetime64('2017-12-31')) == 365

    print convertMetricUnit(100, 'cm', 'm')
    print convertMetricUnit(1, 'm', 'um')