temporalprofiles.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 qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
from qgis.PyQt.QtWidgets import *
import numpy as np
from osgeo import ogr, osr, gdal
from .externals import pyqtgraph as pg
from .externals.pyqtgraph import functions as fn, AxisItem
from .externals.qps.plotstyling.plotstyling import PlotStyle

from .timeseries import TimeSeries, TimeSeriesDate, SensorInstrument, TimeSeriesSource
from .pixelloader import PixelLoader, PixelLoaderTask
from .utils import *
from .externals.qps.speclib.spectrallibraries import createQgsField


LABEL_EXPRESSION_2D = 'DN or Index'
LABEL_TIME = 'Date'
DEBUG = False
OPENGL_AVAILABLE = False
DEFAULT_SAVE_PATH = None
DEFAULT_CRS = QgsCoordinateReferenceSystem('EPSG:4326')

FN_ID = 'fid'
FN_X = 'x'
FN_Y = 'y'
FN_NAME = 'name'

FN_DOY = 'DOY'
FN_DTG = 'DTG'
FN_IS_NODATA ='is_nodata'
FN_GEO_X = 'geo_x'
FN_GEO_Y = 'geo_y'
FN_PX_X = 'px_x'
FN_PX_Y = 'px_y'

#FN_N_TOTAL = 'n'
#FN_N_NODATA = 'no_data'
#FN_N_LOADED = 'loaded'
#FN_N_LOADED_PERCENT = 'percent'

regBandKey = re.compile(r"(?<!\w)b\d+(?!\w)", re.IGNORECASE)
regBandKeyExact = re.compile(r'^' + regBandKey.pattern + '$', re.IGNORECASE)



try:
    import OpenGL

    OPENGL_AVAILABLE = True

except:
    pass



def temporalProfileFeatureFields(sensor: SensorInstrument, singleBandOnly=False) -> QgsFields:
    """
    Returns the fields of a single temporal profile
    :return:
    """
    assert isinstance(sensor, SensorInstrument)
    fields = QgsFields()
    fields.append(createQgsField(FN_DTG, '2011-09-12', comment='Date-time-group'))
    fields.append(createQgsField(FN_DOY, 42, comment='Day-of-year'))
    fields.append(createQgsField(FN_GEO_X, 12.1233, comment='geo-coordinate x/east value'))
    fields.append(createQgsField(FN_GEO_Y, 12.1233, comment='geo-coordinate y/north value'))
    fields.append(createQgsField(FN_PX_X, 42, comment='pixel-coordinate x index'))
    fields.append(createQgsField(FN_PX_Y, 24, comment='pixel-coordinate y index'))


    for b in range(sensor.nb):
        bandKey = bandIndex2bandKey(b)
        fields.append(createQgsField(bandKey, 1.0, comment='value band {}'.format(b+1)))

    return fields

def sensorExampleQgsFeature(sensor:SensorInstrument, singleBandOnly=False)->QgsFeature:
    """
    Returns an exemplary QgsFeature with value for a specific sensor
    :param sensor: SensorInstrument
    :param singleBandOnly:
    :return:
    """
    # populate with exemplary band values (generally stored as floats)

    fields = temporalProfileFeatureFields(sensor)
    f = QgsFeature(fields)
    pt = QgsPointXY(12.34567, 12.34567)
    f.setGeometry(QgsGeometry.fromPointXY(pt))
    f.setAttribute(FN_GEO_X, pt.x())
    f.setAttribute(FN_GEO_Y, pt.y())
    f.setAttribute(FN_PX_X, 1)
    f.setAttribute(FN_PX_Y, 1)
    dtg = datetime.date.today()
    doy = dateDOY(dtg)
    f.setAttribute(FN_DTG, str(dtg))
    f.setAttribute(FN_DOY, doy)

    for b in range(sensor.nb):
        bandKey = bandIndex2bandKey(b)
        f.setAttribute(bandKey, 1.0)

    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 bandIndex2bandKey(i : int):
    assert i >= 0