temporalprofiles2d.py

# -*- coding: utf-8 -*-
"""
/***************************************************************************
                              EO Time Series Viewer
                              -------------------
        begin                : 2017-08-04
        git sha              : $Format:%H$
        copyright            : (C) 2017 by HU-Berlin
        email                : benjamin.jakimow@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.                                   *
 *                                                                         *
 ***************************************************************************/
"""
# noinspection PyPep8Naming

import os, sys, pickle, datetime, re, collections

from collections import OrderedDict
from qgis.gui import *
from qgis.core import *
from PyQt5.Qt import *
from PyQt5.QtCore import *
from PyQt5.QtXml import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
import numpy as np
import pyqtgraph as pg
from pyqtgraph import functions as fn
from pyqtgraph import AxisItem
from osgeo import ogr, osr, gdal
from timeseriesviewer.timeseries import TimeSeries, TimeSeriesDatum, SensorInstrument
from timeseriesviewer.plotstyling import PlotStyle
from timeseriesviewer.pixelloader import PixelLoader, PixelLoaderTask
from timeseriesviewer.utils import *
from timeseriesviewer.models import OptionListModel, Option
LABEL_EXPRESSION_2D = 'DN or Index'
LABEL_TIME = 'Date'
DEBUG = False
OPENGL_AVAILABLE = False

try:
    import OpenGL

    OPENGL_AVAILABLE = True

except:
    pass


def qgsFieldFromKeyValue(fieldName, value):
    t = type(value)
    if t in [int, float] or np.isreal(value):

        fLen  = 0
        fPrec = 0
        fComm = ''
        fType = ''
        f = QgsField(fieldName, QVariant.Double, 'double', 40, 5)
    else:
        f = QgsField(fieldName, QVariant.String, 'text', 40, 5)
    return f


def sensorExampleQgsFeature(sensor, singleBandOnly=False):
    # populate with exemplary band values (generally stored as floats)

    if sensor is None:
        singleBandOnly = True

    fieldValues = collections.OrderedDict()
    if singleBandOnly:
        fieldValues['b'] = 1.0
    else:
        assert isinstance(sensor, SensorInstrument)
        for b in range(sensor.nb):
            fn = bandIndex2bandKey(b)
            fieldValues[fn] = 1.0

    date = datetime.date.today()
    doy = dateDOY(date)
    fieldValues['doy'] = doy
    fieldValues['date'] = str(date)


    fields = QgsFields()
    for k, v in fieldValues.items():
        fields.append(qgsFieldFromKeyValue(k,v))
    f = QgsFeature(fields)
    for k, v in fieldValues.items():
        f.setAttribute(k, v)
    return f


def dateDOY(date):
    if isinstance(date, np.datetime64):
        date = date.astype(datetime.date)
    return date.timetuple().tm_yday

def daysPerYear(year):

    if isinstance(year, np.datetime64):
        year = year.astype(datetime.date)
    if isinstance(year, datetime.date):
        year = year.timetuple().tm_year

    return dateDOY(datetime.date(year=year, month=12, day=31))

def date2num(d):
    #kindly taken from https://stackoverflow.com/questions/6451655/python-how-to-convert-datetime-dates-to-decimal-years
    if isinstance(d, np.datetime64):
        d = d.astype(datetime.datetime)

    if isinstance(d, QDate):
        d = datetime.date(d.year(), d.month(), d.day())

    assert isinstance(d, datetime.date)

    yearDuration = daysPerYear(d)
    yearElapsed = d.timetuple().tm_yday
    fraction = float(yearElapsed) / float(yearDuration)
    if fraction == 1.0:
        fraction = 0.9999999
    return float(d.year) + fraction

def num2date(n, dt64=True, qDate=False):
    n = float(n)
    if n < 1:
        n += 1

    year = int(n)
    fraction = n - year
    yearDuration = daysPerYear(year)
    yearElapsed = fraction * yearDuration

    import math
    doy = round(yearElapsed)
    if doy < 1:
        doy = 1
    try:
        date = datetime.date(year, 1, 1) + datetime.timedelta(days=doy-1)
    except:
        s = ""
    if qDate:
        return QDate(date.year, date.month, date.day)
    if dt64:
        return np.datetime64(date)
    else:
        return date



    #return np.datetime64('{:04}-01-01'.format(year), 'D') + np.timedelta64(int(yearElapsed), 'D')


def saveTemporalProfiles(profiles, path, mode='all', sep=',', loadMissingValues=False):
    if path is None or len(path) == 0:
        return

    assert mode in ['coordinate','all']

    nbMax = 0
    for sensor in profiles[0].timeSeries().sensors():
        assert isinstance(sensor, SensorInstrument)
        nbMax = max(nbMax, sensor.nb)

    ext = os.path.splitext(path)[1].lower()

    assert isinstance(ext, str)
    if ext.startswith('.'):
        ext = ext[1:]

    if loadMissingValues:
        for p in profiles:
            assert isinstance(p, TemporalProfile)
            p.loadMissingData()


    if ext == 'csv':
        #write a flat list of profiles
        lines = ['Temporal Profiles']
        lines.append(sep.join(['pid', 'name', 'date', 'img', 'location', 'band values']))

        for nP, p in enumerate(profiles):
            assert isinstance(p, TemporalProfile)
            lines.append('Profile {} "{}": {}'.format(p.mID, p.name(), p.mCoordinate))
            assert isinstance(p, TemporalProfile)
            c = p.coordinate()
            for i, tsd in enumerate(p.mTimeSeries):
                assert isinstance(tsd, TimeSeriesDatum)

                data = p.data(tsd)
                line = [p.id(), p.name(), data['date'], tsd.pathImg, c.asWkt()]
                for b in range(tsd.sensor.nb):