diff --git a/CHANGES.txt b/CHANGES.txt
index 7c1936e51fac49f88a2699d1e40689c94824a29e..66ef74f4c5bcbb6a60bf828bb3cb051a729e7608 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -9,8 +9,6 @@
- Width of plot lines now scale-independent (issue #64, QPen.setCosmetic(True))
- adding fields to spectral library (issue #61)
-
-
2018-06-04
increased version to 0.6
SpectralLibrary Module
diff --git a/test/test_mapcanvas.py b/test/test_mapcanvas.py
index cebeb5143ec0ac54681705d4cef605e2e7165dd4..e7c7ee4d738d60d452ed117b740633c0bcf34f6d 100644
--- a/test/test_mapcanvas.py
+++ b/test/test_mapcanvas.py
@@ -25,7 +25,8 @@ import unittest, tempfile
from timeseriesviewer.mapcanvas import *
-QGIS_APP = initQgisApplication()
+resourceDir = os.path.join(DIR_REPO,'qgisresources')
+QGIS_APP = initQgisApplication(qgisResourceDir=resourceDir)
class testclassDialogTest(unittest.TestCase):
@@ -47,7 +48,6 @@ class testclassDialogTest(unittest.TestCase):
self.assertIsInstance(m.mLayerModel, MapCanvasLayerModel)
- m.saveAsImage()
if __name__ == "__main__":
diff --git a/test/test_stackedbandinput.py b/test/test_stackedbandinput.py
new file mode 100644
index 0000000000000000000000000000000000000000..b2c68f8db01d4ea5305a8ed02bfbb7878ffa7cef
--- /dev/null
+++ b/test/test_stackedbandinput.py
@@ -0,0 +1,219 @@
+# -*- coding: utf-8 -*-
+
+"""
+***************************************************************************
+
+ ---------------------
+ Date :
+ Copyright : (C) 2017 by Benjamin Jakimow
+ Email : benjamin jakimow at geo dot hu-berlin dot 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
+
+from timeseriesviewer.utils import initQgisApplication, nextColor
+from PyQt5.QtGui import *
+from PyQt5.QtCore import *
+import unittest, tempfile
+
+from timeseriesviewer.stackedbandinput import *
+from example.Images import Img_2014_06_16_LE72270652014167CUB00_BOA, Img_2014_05_07_LC82270652014127LGN00_BOA
+resourceDir = os.path.join(DIR_REPO,'qgisresources')
+QGIS_APP = initQgisApplication(qgisResourceDir=resourceDir)
+
+
+class testclassDialogTest(unittest.TestCase):
+ """Test rerources work."""
+
+ @classmethod
+ def setUpClass(cls):
+ pass
+ #cls.srcDir = r'F:\Temp\EOTSV_Dev\DF'
+ #cls.stackFiles = file_search(cls.srcDir, '*.tif')
+
+ def setUp(self):
+ """Runs before each test."""
+ pass
+
+ def tearDown(self):
+ """Runs after each test."""
+ pass
+
+ def createTestDatasets(self):
+
+ vsiDir = '/vsimem/tmp'
+ from timeseriesviewer.temporalprofiles2d import date2num
+ ns = 50
+ nl = 100
+
+ r1 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(16, 'D'), dtype=np.datetime64)
+ r2 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(8, 'D'), dtype=np.datetime64)
+ drv = gdal.GetDriverByName('ENVI')
+
+ crs = osr.SpatialReference()
+ crs.ImportFromEPSG(32633)
+
+ assert isinstance(drv, gdal.Driver)
+ datasets = []
+ for i, r in enumerate([r1, r2]):
+ p = '{}tmpstack{}.bsq'.format(vsiDir, i+1)
+ nb = len(r)
+ ds = drv.Create(p, ns, nl, nb, eType=gdal.GDT_Float32)
+ assert isinstance(ds, gdal.Dataset)
+
+ ds.SetProjection(crs.ExportToWkt())
+
+ dateString = ','.join([str(d) for d in r])
+ dateString = '{{{}}}'.format(dateString)
+ ds.SetMetadataItem('wavelength', dateString, 'ENVI')
+
+ for b, date in enumerate(r):
+ decimalYear = date2num(date)
+
+ band = ds.GetRasterBand(b+1)
+ assert isinstance(band, gdal.Band)
+ band.Fill(decimalYear)
+ ds.FlushCache()
+ datasets.append(p)
+
+
+ if i == 0:
+ #create a classification image stack
+
+ nc = nb
+ data = np.ones((nb, ns, nl), dtype=np.uint8)
+ classNames = ['unclassified']
+ colorTable = gdal.ColorTable()
+ colorTable.SetColorEntry(0, (0,0,0))
+ assert isinstance(colorTable, gdal.ColorTable)
+ color = QColor('green')
+
+ for j, date in enumerate(r):
+ c = j+1
+ data[j, j:-1, 0:j] = c
+ classNames.append('Class {}'.format(date))
+ colorTable.SetColorEntry(c, color.getRgb())
+ color = nextColor(color)
+
+ p = '{}tmpClassificationStack.bsq'.format(vsiDir)
+ ds = gdal_array.SaveArray(data,p, format='ENVI', prototype=datasets[0])
+ ds.GetRasterBand(1).SetColorTable(colorTable)
+ ds.GetRasterBand(1).SetCategoryNames(classNames)
+
+ assert isinstance(ds, gdal.Dataset)
+ ds.SetMetadataItem('wavelength', dateString, 'ENVI')
+ ds.FlushCache()
+ datasets.append(ds)
+ return datasets
+
+
+ def test_inputmodel(self):
+ testData = self.createTestDatasets()
+ m = InputStackTableModel()
+ m.insertSources(testData)
+ self.assertTrue(len(m) == len(testData))
+
+ dTotal, dIntersecton = m.dateInfo()
+
+ self.assertTrue(len(dTotal) > 0)
+ self.assertTrue(len(dIntersecton) > 0)
+ self.assertTrue(len(dTotal) > len(dIntersecton))
+
+
+ def test_outputmodel(self):
+
+ m = OutputImageModel()
+ m.setOutputDir('/vsimem/dub')
+ m.setOutputPrefix('myPrefix')
+
+ testData = self.createTestDatasets()
+
+ stackInfo = InputStackInfo(testData[0])
+
+ self.assertTrue(len(stackInfo) > 0)
+
+ stackInfos = [InputStackInfo(f) for f in testData]
+
+ dates = set()
+
+ for s in stackInfos:
+ dates.update(s.dates())
+
+ m.setMultiStackSources(stackInfos, list(dates))
+
+ self.assertTrue(len(m) > 0)
+
+ outInfo = m.mOutputImages[0]
+ self.assertIsInstance(outInfo, OutputVRTDescription)
+
+ xml = m.vrtXML(outInfo)
+ self.assertIsInstance(xml, str)
+
+ eTree = m.vrtXML(outInfo, asElementTree=True)
+ self.assertIsInstance(eTree, ElementTree.Element)
+
+
+ def test_dateparsing(self):
+
+ dsDates = gdal.OpenEx(self.stackFiles[1], allowed_drivers=['ENVI'])
+
+ #dsDates = gdal.Open(self.stackFiles[1])
+ dates = datesFromDataset(dsDates)
+ self.assertEqual(len(dates), dsDates.RasterCount)
+
+ dsNoDates = gdal.OpenEx(self.stackFiles[0], allowed_drivers=['ENVI'])
+ dates = datesFromDataset(dsNoDates)
+ self.assertEqual(len(dates), 0)
+
+ s = ""
+
+ def test_dialog(self):
+ d = StackedBandInputDialog()
+ d.addSources(self.createTestDatasets())
+
+ r = d.exec_()
+
+ self.assertTrue(r in [QDialog.Rejected, QDialog.Accepted])
+ if r == QDialog.Accepted:
+ images = d.writtenFiles()
+ self.assertTrue(len(images) > 0)
+
+ for p in images:
+ ds = gdal.Open(p)
+ self.assertIsInstance(ds, gdal.Dataset)
+ lyr = QgsRasterLayer(p)
+ self.assertIsInstance(lyr, QgsRasterLayer)
+ self.assertTrue(lyr.isValid())
+ else:
+ self.assertTrue(len(d.writtenFiles()) == 0)
+
+ #QGIS_APP.exec_()
+ pass
+
+ def test_withTSV(self):
+
+
+
+ testImages = self.createTestDatasets()
+ from timeseriesviewer.main import TimeSeriesViewer
+ TSV = TimeSeriesViewer(None)
+ TSV.show()
+
+ d = StackedBandInputDialog()
+ d.addSources(self.createTestDatasets())
+ writtenFiles = d.saveImages()
+ self.assertTrue(len(writtenFiles) > 0)
+ TSV.loadImageFiles(writtenFiles)
+
+ self.assertTrue(len(TSV.TS) == len(writtenFiles))
+
+ QGIS_APP.exec_()
+if __name__ == "__main__":
+ unittest.main()
diff --git a/test/test_timeseries.py b/test/test_timeseries.py
index e470194e2a6300e5dc7478954a82328c8909bf8b..e6160d408e92698eb294158a6e7c0758b7d62f65 100644
--- a/test/test_timeseries.py
+++ b/test/test_timeseries.py
@@ -5,12 +5,57 @@ import os
import unittest
import example
import example.Images
-from osgeo import gdal
+from osgeo import gdal, ogr, osr
from timeseriesviewer.utils import file_search
-from timeseriesviewer.timeseries import TimeSeries, SensorInstrument, TimeSeriesDatum
+from timeseriesviewer.timeseries import *
class TestInit(unittest.TestCase):
+
+ def createTestDatasets(self):
+
+ vsiDir = '/vsimem/tmp'
+ from timeseriesviewer.temporalprofiles2d import date2num
+ ns = 50
+ nl = 100
+
+ r1 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(16, 'D'), dtype=np.datetime64)
+ r2 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(8, 'D'), dtype=np.datetime64)
+ drv = gdal.GetDriverByName('ENVI')
+
+ crs = osr.SpatialReference()
+ crs.ImportFromEPSG(32633)
+
+ assert isinstance(drv, gdal.Driver)
+ datasets = []
+
+ for i, r in enumerate([r1, r2]):
+ p = '{}tmpstack{}.bsq'.format(vsiDir, i+1)
+
+ ds = drv.Create(p, ns, nl, len(r), eType=gdal.GDT_Float32)
+ assert isinstance(ds, gdal.Dataset)
+
+ ds.SetProjection(crs.ExportToWkt())
+
+ dateString = ','.join([str(d) for d in r])
+ dateString = '{{{}}}'.format(dateString)
+ ds.SetMetadataItem('wavelength', dateString, 'ENVI')
+
+ for b, date in enumerate(r):
+ decimalYear = date2num(date)
+
+ band = ds.GetRasterBand(b+1)
+ assert isinstance(band, gdal.Band)
+ band.Fill(decimalYear)
+ ds.FlushCache()
+ datasets.append(p)
+
+
+
+
+
+ return datasets
+
def test_timeseriesdatum(self):
file = example.Images.Img_2014_03_20_LC82270652014079LGN00_BOA
@@ -19,6 +64,8 @@ class TestInit(unittest.TestCase):
self.assertIsInstance(tsd, TimeSeriesDatum)
self.assertEqual(tsd.nb, 6)
+
+
def test_timeseries(self):
files = file_search(os.path.dirname(example.__file__), '*.tif', recursive=True)
@@ -48,7 +95,7 @@ class TestInit(unittest.TestCase):
dsRE = gdal.Open(pathRE)
assert isinstance(dsRE, gdal.Dataset)
- tsdRE = TS.tsdFromPath(pathRE)
+ tsdRE = TS.getTSD(pathRE)
self.assertIsInstance(tsdRE, TimeSeriesDatum)
sRE = tsdRE.sensor
self.assertIsInstance(sRE, SensorInstrument)
diff --git a/timeseriesviewer/tests.py b/test/tests.py
similarity index 77%
rename from timeseriesviewer/tests.py
rename to test/tests.py
index 74d02ed8f8e994c1d696200ebcc8c37ffadd6a40..ffeca468f2ba773233363f351a2e2763d178b110 100644
--- a/timeseriesviewer/tests.py
+++ b/test/tests.py
@@ -28,16 +28,62 @@ from timeseriesviewer import *
from timeseriesviewer.timeseries import *
from timeseriesviewer import DIR_EXAMPLES
from timeseriesviewer.utils import file_search
+from osgeo import ogr, osr, gdal, gdal_array
+
+
+
class TestObjects():
+ """
+ Creates objects to be used for testing. It is prefered to generate objects in-memory.
+ """
@staticmethod
- def timeSeries():
+ def createTimeSeries():
TS = TimeSeries()
files = file_search(DIR_EXAMPLES, '*.bsq', recursive=True)
TS.addFiles(files)
return TS
+ @staticmethod
+ def createTimeSeriesStacks():
+ vsiDir = '/vsimem/tmp'
+ from timeseriesviewer.temporalprofiles2d import date2num
+ ns = 50
+ nl = 100
+
+ r1 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(16, 'D'), dtype=np.datetime64)
+ r2 = np.arange('2000-01-01', '2005-06-14', step=np.timedelta64(8, 'D'), dtype=np.datetime64)
+ drv = gdal.GetDriverByName('ENVI')
+
+ crs = osr.SpatialReference()
+ crs.ImportFromEPSG(32633)
+
+ assert isinstance(drv, gdal.Driver)
+ datasets = []
+ for i, r in enumerate([r1, r2]):
+ p = '{}tmpstack{}.bsq'.format(vsiDir, i + 1)
+
+ ds = drv.Create(p, ns, nl, len(r), eType=gdal.GDT_Float32)
+ assert isinstance(ds, gdal.Dataset)
+
+ ds.SetProjection(crs.ExportToWkt())
+
+ dateString = ','.join([str(d) for d in r])
+ dateString = '{{{}}}'.format(dateString)
+ ds.SetMetadataItem('wavelength', dateString, 'ENVI')
+
+ for b, date in enumerate(r):
+ decimalYear = date2num(date)
+
+ band = ds.GetRasterBand(b + 1)
+ assert isinstance(band, gdal.Band)
+ band.Fill(decimalYear)
+ ds.FlushCache()
+ datasets.append(p)
+
+ return datasets
+
@staticmethod
def spectralProfiles(n):
"""
diff --git a/timeseriesviewer/dateparser.py b/timeseriesviewer/dateparser.py
index 110de79f58a5090da1a0c5ea8c370aa3b068371c..b55da28783e1d7048d8da4d82f193035ff736db5 100644
--- a/timeseriesviewer/dateparser.py
+++ b/timeseriesviewer/dateparser.py
@@ -1,9 +1,10 @@
-import os, re, logging
+import os, re, logging, datetime
from osgeo import gdal
import numpy as np
from qgis import *
-logger = logging.getLogger(__file__)
+from qgis.PyQt.QtCore import QDate
+
#regular expression. compile them only once
@@ -19,6 +20,7 @@ regYYYYMMDD = re.compile(r'(?P<year>(19|20)\d\d)(?P<hyphen>-?)(?P<month>1[0-2]|0
regMissingHypen = re.compile(r'^\d{8}')
regYYYYMM = re.compile(r'([0-9]{4})-(1[0-2]|0[1-9])')
regYYYYDOY = re.compile(r'(?P<year>(19|20)\d\d)-?(?P<day>36[0-6]|3[0-5][0-9]|[12][0-9]{2}|0[1-9][0-9]|00[1-9])')
+regDecimalYear = re.compile(r'(?P<year>(19|20)\d\d)\.(?P<datefraction>\d\d\d)')
def matchOrNone(regex, text):
match = regex.search(text)
@@ -27,7 +29,52 @@ def matchOrNone(regex, text):
else:
return None
-def extractDateTimeGroup(text):
+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 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
+
+
+def extractDateTimeGroup(text:str)->np.datetime64:
+ """
+ Extracts a date-time-group from a text string
+ :param text: a string
+ :return: numpy.datetime64 in case of success, or None
+ """
match = regISODate1.search(text)
if match:
matchedText = match.group()
@@ -47,6 +94,14 @@ def extractDateTimeGroup(text):
if match:
return np.datetime64(match.group())
+ match = regDecimalYear.search(text)
+ if match:
+ year = float(match.group('year'))
+ df = float(match.group('datefraction'))
+ num = match.group()
+ return num2date(num)
+
+
return None
def datetime64FromYYYYMMDD(yyyymmdd):
diff --git a/timeseriesviewer/main.py b/timeseriesviewer/main.py
index 7bf76adcd4685ad835417455a33fd4d29fadf49c..e090f3e7276c6fdb42f85cace7c48d474665e50b 100644
--- a/timeseriesviewer/main.py
+++ b/timeseriesviewer/main.py
@@ -463,7 +463,7 @@ class TimeSeriesViewer(QgisInterface, QObject):
self.ui.actionClearTS.triggered.connect(self.clearTimeSeries)
self.ui.actionSaveTS.triggered.connect(self.saveTimeSeriesDefinition)
self.ui.actionAddTSExample.triggered.connect(self.loadExampleTimeSeries)
-
+ self.ui.actionLoadTimeSeriesStack.triggered.connect(self.loadTimeSeriesStack)
self.ui.actionShowCrosshair.toggled.connect(self.spatialTemporalVis.setShowCrosshair)
#connect buttons with actions
@@ -674,6 +674,15 @@ class TimeSeriesViewer(QgisInterface, QObject):
if path is not None:
s.setValue('FILE_TS_DEFINITION', path)
+ def loadTimeSeriesStack(self):
+
+ from timeseriesviewer.stackedbandinput import StackedBandInputDialog
+
+ d = StackedBandInputDialog(parent=self.ui)
+ if d.exec_() == QDialog.Accepted:
+ C = d.saveImages()
+ self.addTimeSeriesImages(writtenFiles)
+
def loadExampleTimeSeries(self):
diff --git a/timeseriesviewer/mapvisualization.py b/timeseriesviewer/mapvisualization.py
index e6e32efff49d2d2830ecdd3a4872791047df6023..53324cda18d0d221ab21a645485d53239cac66ba 100644
--- a/timeseriesviewer/mapvisualization.py
+++ b/timeseriesviewer/mapvisualization.py
@@ -42,7 +42,6 @@ from timeseriesviewer.crosshair import CrosshairStyle
-#dummyPath = os.path.abspath(os.path.join(os.path.dirname(__file__), 'dummy_gdal.vrt'))
#assert os.path.isfile(dummyPath)
#lyr = QgsRasterLayer(dummyPath)
#assert lyr.isValid()
@@ -1532,21 +1531,65 @@ RENDER_CLASSES['renderer-v2'] = {
}
+
+
def rendererToXml(renderer):
+ """
+ Returns a renderer XML representation
+ :param renderer: QgsRasterRender | QgsFeatureRenderer
+ :return: QDomDocument
+ """
doc = QDomDocument()
+ err = ''
if isinstance(renderer, QgsRasterRenderer):
- path = os.path.join(os.path.dirname(__file__), 'dummy_gdal.vrt')
+ #create a dummy raster layer
+ import uuid
+ from timeseriesviewer.virtualrasters import write_vsimem, read_vsimem
+ xml = """<VRTDataset rasterXSize="1" rasterYSize="1">
+ <GeoTransform> 0.0000000000000000e+00, 1.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.0000000000000000e+00</GeoTransform>
+ <VRTRasterBand dataType="Float32" band="1">
+ <Metadata>
+ <MDI key="STATISTICS_MAXIMUM">0</MDI>
+ <MDI key="STATISTICS_MEAN">0</MDI>
+ <MDI key="STATISTICS_MINIMUM">0</MDI>
+ <MDI key="STATISTICS_STDDEV">0</MDI>
+ </Metadata>
+ <Description>Band 1</Description>
+ <Histograms>
+ <HistItem>
+ <HistMin>0</HistMin>
+ <HistMax>0</HistMax>
+ <BucketCount>1</BucketCount>
+ <IncludeOutOfRange>0</IncludeOutOfRange>
+ <Approximate>0</Approximate>
+ <HistCounts>0</HistCounts>
+ </HistItem>
+ </Histograms>
+ </VRTRasterBand>
+ </VRTDataset>
+ """
+ path = '/vsimem/{}.vrt'.format(uuid.uuid4())
+ drv = gdal.GetDriverByName('VRT')
+ assert isinstance(drv, gdal.Driver)
+ write_vsimem(path, xml)
lyr = QgsRasterLayer(path)
+ assert lyr.isValid()
lyr.setRenderer(renderer.clone())
+ lyr.exportNamedStyle(doc, err)
+ #remove dummy raster layer
+ lyr = None
+ drv.Delete(path)
elif isinstance(renderer, QgsFeatureRenderer):
#todo: distinguish vector type from requested renderer
lyr = QgsVectorLayer('Point?crs=epsg:4326&field=id:integer', 'dummy', 'memory')
lyr.setRenderer(renderer.clone())
+ lyr.exportNamedStyle(doc, err)
+ lyr = None
else:
raise NotImplementedError()
- err = ''
- lyr.exportNamedStyle(doc, err)
+
+
return doc
diff --git a/timeseriesviewer/pixelloader.py b/timeseriesviewer/pixelloader.py
index ccf65d3d4deac6edbd03728ab109bb3ffc04b234..dd828b011dba2f8ee5351423e9072bb076b84b98 100644
--- a/timeseriesviewer/pixelloader.py
+++ b/timeseriesviewer/pixelloader.py
@@ -67,7 +67,7 @@ class PixelLoaderTask(object):
def fromDump(byte_object):
return pickle.loads(byte_object)
- def __init__(self, source, geometries, bandIndices=None, **kwargs):
+ def __init__(self, source:str, geometries, bandIndices=None, **kwargs):
"""
:param jobId: jobId number as given by the calling PixelLoader
@@ -83,8 +83,6 @@ class PixelLoaderTask(object):
for geometry in geometries:
assert type(geometry) in [SpatialExtent, SpatialPoint]
- assert os.path.isfile(source)
-
#assert isinstance(source, str) or isinstance(source, unicode)
self.sourcePath = source
diff --git a/timeseriesviewer/profilevisualization.py b/timeseriesviewer/profilevisualization.py
index 88bb882618c0d69f2fb498935a9d702099b90943..08d885f629dae0df4b8634e419d16924d74e51ca 100644
--- a/timeseriesviewer/profilevisualization.py
+++ b/timeseriesviewer/profilevisualization.py
@@ -30,7 +30,7 @@ from qgis.PyQt.QtGui import *
from timeseriesviewer import jp, SETTINGS
from timeseriesviewer.timeseries import *
-from timeseriesviewer.utils import SpatialExtent, SpatialPoint, px2geo, loadUI
+from timeseriesviewer.utils import SpatialExtent, SpatialPoint, px2geo, loadUI, nextColor
from timeseriesviewer.plotstyling import PlotStyle, PlotStyleButton
from timeseriesviewer.pixelloader import PixelLoader, PixelLoaderTask
from timeseriesviewer.sensorvisualization import SensorListModel
@@ -1284,35 +1284,6 @@ class ProfileViewDockUI(QgsDockWidget, loadUI('profileviewdock.ui')):
w.update()
self.setWindowTitle(title)
-NEXT_COLOR_HUE_DELTA_CON = 10
-NEXT_COLOR_HUE_DELTA_CAT = 100
-def nextColor(color, mode='cat'):
- """
- Returns another color
- :param color: the previous color
- :param mode: 'cat' - for categorical color jump (next color looks pretty different to previous)
- 'con' - for continuous color jump (next color looks similar to previous)
- :return:
- """
- assert mode in ['cat','con']
- assert isinstance(color, QColor)
- hue, sat, value, alpha = color.getHsl()
- if mode == 'cat':
- hue += NEXT_COLOR_HUE_DELTA_CAT
- elif mode == 'con':
- hue += NEXT_COLOR_HUE_DELTA_CON
- if sat == 0:
- sat = 255
- value = 128
- alpha = 255
- s = ""
- while hue > 360:
- hue -= 360
-
- return QColor.fromHsl(hue, sat, value, alpha)
-
-
-
class SpectralTemporalVisualization(QObject):
diff --git a/timeseriesviewer/spectrallibraries.py b/timeseriesviewer/spectrallibraries.py
index b421f716909bd26e052c31e7b088739fb718ee45..7e99bcb5a848b662fd0dcb5d222be975c8d1ca41 100644
--- a/timeseriesviewer/spectrallibraries.py
+++ b/timeseriesviewer/spectrallibraries.py
@@ -35,10 +35,10 @@ from pyqtgraph.graphicsItems.PlotDataItem import PlotDataItem
from pyqtgraph.graphicsItems.PlotItem import PlotItem
import pyqtgraph.functions as fn
import numpy as np
-from osgeo import gdal, gdal_array
+from osgeo import gdal, gdal_array, ogr
from timeseriesviewer.utils import *
-#from timeseriesviewer.virtualrasters import *
+from timeseriesviewer.virtualrasters import describeRawFile
from timeseriesviewer.models import *
from timeseriesviewer.plotstyling import PlotStyle, PlotStyleDialog, MARKERSYMBOLS2QGIS_SYMBOLS
import timeseriesviewer.mimedata as mimedata
@@ -1863,7 +1863,7 @@ class EnviSpectralLibraryIO(AbstractSpectralLibraryIO):
xSize = int(hdr['samples'])
ySize = int(hdr['lines'])
bands = int(hdr['bands'])
- byteOrder = 'MSB' if hdr['byte order'] == 0 else 'LSB'
+ byteOrder = 'LSB' if int(hdr['byte order']) == 0 else 'MSB'
if pathVrt is None:
id = uuid.UUID()
@@ -2556,36 +2556,6 @@ class SpectralLibraryPlotWidget(PlotWidget):
#if self.mModel.rowCount() > 0:
# self.onRowsInserted(self.mModel.index(0,0), 0, self.mModel.rowCount())
- def onMouseMoved2D(self, evt):
- pos = evt[0] ## using signal proxy turns original arguments into a tuple
-
- plotItem = self.getPlotItem()
- if plotItem.sceneBoundingRect().contains(pos):
- vb = plotItem.vb
- assert isinstance(vb, DateTimeViewBox)
- mousePoint = vb.mapSceneToView(pos)
- x = mousePoint.x()
- if x >= 0:
- y = mousePoint.y()
- date = num2date(x)
- doy = dateDOY(date)
- plotItem.vb.updateCurrentDate(num2date(x, dt64=True))
- self.mInfoLabelCursor.setText('DN {:0.2f}\nDate {}\nDOY {}'.format(
- mousePoint.y(), date, doy),
- color=self.mInfoColor)
-
- s = self.size()
- pos = QPointF(s.width(), 0)
- self.mInfoLabelCursor.setVisible(vb.mActionShowCursorValues.isChecked())
- self.mInfoLabelCursor.setPos(pos)
-
- b = vb.mActionShowCrosshair.isChecked()
- self.mCrosshairLineH.setVisible(b)
- self.mCrosshairLineV.setVisible(b)
- self.mCrosshairLineH.pen.setColor(self.mInfoColor)
- self.mCrosshairLineV.pen.setColor(self.mInfoColor)
- self.mCrosshairLineV.setPos(mousePoint.x())
- self.mCrosshairLineH.setPos(mousePoint.y())
def speclib(self):
diff --git a/timeseriesviewer/stackedbandinput.py b/timeseriesviewer/stackedbandinput.py
new file mode 100644
index 0000000000000000000000000000000000000000..8cc4e547afba109af751da9566b4d9cb9560ba70
--- /dev/null
+++ b/timeseriesviewer/stackedbandinput.py
@@ -0,0 +1,832 @@
+# -*- coding: utf-8 -*-
+# noinspection PyPep8Naming
+"""
+***************************************************************************
+ stackedbandinput.py
+
+ Sometimes time-series-data is written out as stacked band images, having one observation per band.
+ This module helps to use such data as EOTS input.
+ ---------------------
+ Date : June 2018
+ Copyright : (C) 2018 by Benjamin Jakimow
+ 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. *
+* *
+***************************************************************************
+"""
+
+import os, re, tempfile, pickle, copy, shutil, locale, uuid, csv, io
+from xml.etree import ElementTree
+from collections import OrderedDict
+from qgis.core import *
+from qgis.gui import *
+from qgis.utils import qgsfunction
+from qgis.PyQt.QtCore import *
+from qgis.PyQt.QtGui import *
+from qgis.PyQt.QtWidgets import *
+from qgis.core import QgsField, QgsFields, QgsFeature, QgsMapLayer, QgsVectorLayer, QgsConditionalStyle
+from qgis.gui import QgsMapCanvas, QgsDockWidget
+from pyqtgraph.widgets.PlotWidget import PlotWidget
+from pyqtgraph.graphicsItems.PlotDataItem import PlotDataItem
+from pyqtgraph.graphicsItems.PlotItem import PlotItem
+import pyqtgraph.functions as fn
+import numpy as np
+from osgeo import gdal, gdal_array
+import numpy as np
+from timeseriesviewer.utils import *
+from timeseriesviewer.virtualrasters import *
+from timeseriesviewer.models import *
+from timeseriesviewer.dateparser import *
+from timeseriesviewer.plotstyling import PlotStyle, PlotStyleDialog, MARKERSYMBOLS2QGIS_SYMBOLS
+import timeseriesviewer.mimedata as mimedata
+
+
+def datesFromDataset(dataset:gdal.Dataset)->list:
+
+ nb = dataset.RasterCount
+
+ def checkDates(dateList):
+ if not len(dateList) == nb:
+ return False
+ for d in dateList:
+ if not isinstance(d, np.datetime64):
+ return False
+ return True
+
+ searchedKeys = []
+ searchedKeys.append(re.compile('aquisition[ ]*dates$', re.I))
+ searchedKeys.append(re.compile('observation[ ]*dates$', re.I))
+ searchedKeys.append(re.compile('wavelength$', re.I))
+
+ #1. Check Metadata
+ for domain in dataset.GetMetadataDomainList():
+ domainData = dataset.GetMetadata_Dict(domain)
+ assert isinstance(domainData, dict)
+
+ for key, values in domainData.items():
+ for regex in searchedKeys:
+ if regex.search(key.strip()):
+ values = re.sub('[{}]', '', values)
+ values = values.split(',')
+ dateValues = [extractDateTimeGroup(t) for t in values]
+ if checkDates(dateValues):
+ return dateValues
+
+
+ #2. Check Band Names
+ bandDates = [extractDateTimeGroup(dataset.GetRasterBand(b+1).GetDescription()) for b in range(nb)]
+ bandDates = [b for b in bandDates if isinstance(b, np.datetime64)]
+ if checkDates(bandDates):
+ return bandDates
+
+ return []
+
+class InputStackInfo(object):
+
+ def __init__(self, dataset):
+ if isinstance(dataset, str):
+ #test ENVI header first
+ basename = os.path.splitext(dataset)[0]
+ ds = None
+ if os.path.isfile(basename+'.hdr'):
+ ds = gdal.OpenEx(dataset, allowed_drivers=['ENVI'])
+ if not isinstance(ds, gdal.Dataset):
+ ds = gdal.Open(dataset)
+ if not isinstance(ds, gdal.Dataset):
+ raise Exception('Unable to open {}'.format(dataset))
+
+ dataset = ds
+ del ds
+
+ assert isinstance(dataset, gdal.Dataset)
+
+ self.mMetadataDomains = dataset.GetMetadataDomainList()
+ self.mMetaData = OrderedDict()
+
+ for domain in self.mMetadataDomains:
+ self.mMetaData[domain] = dataset.GetMetadata_Dict(domain)
+
+ self.ns = dataset.RasterXSize
+ self.nl = dataset.RasterYSize
+ self.nb = dataset.RasterCount
+
+ self.wkt = dataset.GetProjection()
+ self.gt = dataset.GetGeoTransform()
+
+ self.colorTable = dataset.GetRasterBand(1).GetColorTable()
+ self.classNames = dataset.GetRasterBand(1).GetCategoryNames()
+
+ self.path = dataset.GetFileList()[0]
+
+ self.outputBandName = os.path.basename(self.path)
+ if len(self.outputBandName) == 0:
+ self.outputBandName = ''
+
+ self.bandnames = []
+ self.nodatavalues = []
+ for b in range(self.nb):
+ band = dataset.GetRasterBand(b+1)
+ assert isinstance(band, gdal.Band)
+ self.bandnames.append(band.GetDescription())
+ self.nodatavalues.append(band.GetNoDataValue())
+
+
+ self.mDates = datesFromDataset(dataset)
+
+ def __len__(self):
+ return len(self.mDates)
+
+ def dates(self)->list:
+ """Returns a list of dates"""
+ return self.mDates
+
+
+ def structure(self):
+ return (self.ns, self.nl, self.nb, self.gt, self.wkt)
+
+ def wavelength(self):
+ return self.mMetaData[''].get('wavelength')
+
+ def setWavelength(self, wl):
+ self.mMetaData['']['wavelength'] = wl
+
+
+class OutputVRTDescription(object):
+ """
+ Descrbies an output VRT
+ """
+
+ def __init__(self, path:str, date:np.datetime64):
+ super(OutputVRTDescription, self).__init__()
+ self.mPath = path
+ self.mDate = date
+
+
+ def setPath(self, path:str):
+ self.mPath = path
+
+
+
+class InputStackTableModel(QAbstractTableModel):
+
+
+
+ def __init__(self, parent=None):
+
+ super(InputStackTableModel, self).__init__(parent)
+ self.mStackImages = []
+
+ self.cn_source = 'Source'
+ self.cn_dates = 'Dates'
+ self.cn_crs = 'GT + CRS'
+ self.cn_ns = 'ns'
+ self.cn_nl = 'nl'
+ self.cn_nb = 'nb'
+ self.cn_name = 'Band Name'
+ self.cn_wl = 'Wavelength'
+ self.mColumnNames = [self.cn_source, self.cn_dates, self.cn_name, self.cn_wl, self.cn_ns, self.cn_nl, self.cn_nb, self.cn_crs]
+
+ self.mColumnTooltips = {}
+ self.mColumnTooltips[self.cn_source] = 'Stack source uri / file path'
+ self.mColumnTooltips[self.cn_crs] = 'Geo-Transformation + Coordinate Reference System'
+ self.mColumnTooltips[self.cn_ns] = 'Number of samples / pixel in horizontal direction'
+ self.mColumnTooltips[self.cn_nl] = 'Number of lines / pixel in vertical direction'
+ self.mColumnTooltips[self.cn_nb] = 'Number of bands'
+ self.mColumnTooltips[self.cn_name] = 'Prefix of band name in output image'
+ self.mColumnTooltips[self.cn_wl] = 'Wavelength in output image'
+ self.mColumnTooltips[self.cn_dates] = 'Identified dates'
+
+ def __len__(self):
+ return len(self.mStackImages)
+
+ def __iter__(self):
+ return iter(self.mStackImages)
+
+ def columnName(self, i) -> str:
+ if isinstance(i, QModelIndex):
+ i = i.column()
+ return self.mColumnNames[i]
+
+
+
+ def dateInfo(self):
+ """
+ Returns a list with all extracted dates and a list of date in common between all datasets
+ :return: [all dates], [dates in common]
+ """
+ if len(self) == 0:
+ return [],[]
+ datesTotal = set()
+ datesInCommon = None
+ for i, f in enumerate(self.mStackImages):
+ assert isinstance(f, InputStackInfo)
+
+ dates = f.dates()
+ if datesInCommon is None:
+ datesInCommon = set(dates)
+ else:
+ datesInCommon = datesInCommon.intersection(dates)
+
+ datesTotal = datesTotal.union(f.dates())
+
+ return sorted(list(datesTotal)), sorted(list(datesInCommon))
+
+ def flags(self, index):
+ if index.isValid():
+ columnName = self.columnName(index)
+ flags = Qt.ItemIsEnabled | Qt.ItemIsSelectable
+ if columnName in [self.cn_name, self.cn_wl]: #allow check state
+ flags = flags | Qt.ItemIsEditable
+
+ 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:
+ cname = self.mColumnNames[col]
+ if role == Qt.DisplayRole:
+ return cname
+ elif role == Qt.ToolTipRole:
+ return self.mColumnTooltips.get(cname)
+ elif orientation == Qt.Vertical and role == Qt.DisplayRole:
+ return col
+ return None
+
+
+ def rowCount(self, parent=None):
+ return len(self.mStackImages)
+
+ def columnCount(self, parent: QModelIndex):
+ return len(self.mColumnNames)
+
+ def insertSources(self, paths, i=None):
+ """
+ Inserts new datasources
+ :param paths: [list-of-datasources]
+ :param i: index where to add the first datasource.
+ """
+
+ if i == None:
+ i = self.rowCount()
+
+ if not isinstance(paths, list):
+ paths = [paths]
+
+ infos = [InputStackInfo(p) for p in paths]
+ if len(infos) > 0:
+
+ self.beginInsertRows(QModelIndex(), i, i+len(infos)-1)
+ for j, info in enumerate(infos):
+ assert isinstance(info, InputStackInfo)
+ if len(info.outputBandName) == 0:
+ info.outputBandName = 'Band {}'.format(i+j+1)
+ self.mStackImages.insert(i+j, info)
+ self.endInsertRows()
+
+ def removeSources(self, stackInfos:list):
+
+ for stackInfo in stackInfos:
+ assert stackInfo in self.mStackImages
+
+ for stackInfo in stackInfos:
+ assert isinstance(stackInfo, InputStackInfo)
+
+ idx = self.info2index(stackInfo)
+
+ self.beginRemoveRows(QModelIndex(), idx.row(), idx.row())
+ self.mStackImages.remove(stackInfo)
+ self.endRemoveRows()
+
+ def isValid(self):
+ l = len(self.mStackImages)
+ if l == 0:
+ return False
+ ref = self.mStackImages[0]
+ assert isinstance(ref, InputStackInfo)
+
+ #all input stacks need to have the same characteristic
+ for stackInfo in self.mStackImages[1:]:
+ assert isinstance(stackInfo, InputStackInfo)
+ if not ref.dates() == stackInfo.dates():
+ return False
+ if not ref.structure() == stackInfo.structure():
+ return False
+ return True
+
+
+ def index2info(self, index:QModelIndex) -> InputStackInfo:
+ return self.mStackImages[index.row()]
+
+ def info2index(self, info:InputStackInfo) -> QModelIndex:
+ r = self.mStackImages.index(info)
+ return self.createIndex(r,0, info)
+
+ def data(self, index: QModelIndex, role: int):
+ if not index.isValid():
+ return None
+
+ info = self.mStackImages[index.row()]
+ assert isinstance(info, InputStackInfo)
+ cname = self.columnName(index)
+
+ if role in [Qt.DisplayRole, Qt.ToolTipRole]:
+ if cname == self.cn_source:
+ return info.path
+ if cname == self.cn_dates:
+ dates = info.dates()
+ if role == Qt.DisplayRole:
+ return len(dates)
+ if role == Qt.ToolTipRole:
+ if len(dates) == 0:
+ return 'No dates identified. Can not use this image as input'
+ else:
+ if len(dates) > 11:
+ dates = dates[0:10] + ['...']
+ return '\n'.join([str(d) for d in dates])
+
+
+ if cname == self.cn_ns:
+ return info.ns
+ if cname == self.cn_nl:
+ return info.nl
+ if cname == self.cn_nb:
+ return info.nb
+ if cname == self.cn_crs:
+ return '{} {}'.format(info.gt, info.wkt)
+ elif cname == self.cn_wl:
+ return info.mMetaData[''].get('wavelength')
+ elif cname == self.cn_name:
+ return info.outputBandName
+
+ if role == Qt.EditRole:
+ if cname == self.cn_wl:
+ return info.mMetaData[''].get('wavelength')
+ elif cname == self.cn_name:
+ return info.outputBandName
+
+ if role == Qt.BackgroundColorRole:
+ if cname in [self.cn_name, self.cn_wl]:
+ return QColor('yellow')
+
+ return None
+
+ def setData(self, index: QModelIndex, value, role: int):
+
+ if not index.isValid():
+ return None
+
+ info = self.index2info(index)
+ cname = self.columnName(index)
+
+ changed = False
+ if role == Qt.EditRole:
+ if cname == self.cn_name:
+ if isinstance(value, str) and len(value) > 0:
+ info.outputBandName = value
+ changed = True
+ elif cname == self.cn_wl:
+ info.setWavelength(value)
+ changed = True
+ if changed:
+ self.dataChanged.emit(index, index)
+ return changed
+
+class OutputImageModel(QAbstractTableModel):
+
+ def __init__(self, parent=None):
+ super(OutputImageModel, self).__init__(parent)
+
+
+ self.cn_uri = 'Path'
+ self.cn_date = 'Date'
+ self.mOutputImages = []
+
+ self.mColumnNames = [self.cn_date, self.cn_uri]
+ self.mColumnTooltips = {}
+ self.mColumnTooltips[self.cn_uri] = 'Output location'
+ self.masterVRT_DateLookup = {}
+ self.masterVRT_SourceBandTemplates = {}
+ self.masterVRT_InputStacks = None
+ self.masterVRT_XML = None
+ self.mOutputDir = '/vsimem/'
+ self.mOutputPrefix = 'date'
+
+
+
+ def headerData(self, col, orientation, role):
+ if Qt is None:
+ return None
+ if orientation == Qt.Horizontal:
+ cname = self.mColumnNames[col]
+ if role == Qt.DisplayRole:
+ return cname
+ elif role == Qt.ToolTipRole:
+ return self.mColumnTooltips.get(cname)
+ elif orientation == Qt.Vertical and role == Qt.DisplayRole:
+ return col
+ return None
+
+ def createVRTUri(self, date:np.datetime64):
+
+ path = os.path.join(self.mOutputDir, self.mOutputPrefix)
+ path = '{}{}.vrt'.format(path, date)
+
+ return path
+
+ def clearOutputs(self):
+ self.beginRemoveRows(QModelIndex(), 0, self.rowCount() - 1)
+ self.mOutputImages = []
+ self.endRemoveRows()
+
+ def setMultiStackSources(self, listOfInputStacks:list, dates:list):
+
+ self.clearOutputs()
+
+ if listOfInputStacks is None or len(listOfInputStacks) == 0:
+ return
+ if dates is None or len(dates) == 0:
+ return
+ for s in listOfInputStacks:
+ assert isinstance(s, InputStackInfo)
+ dates = sorted(dates)
+
+ listOfInputStacks = [s for s in listOfInputStacks if len(s) > 0]
+ numberOfOutputVRTBands = len(listOfInputStacks)
+ self.masterVRT_DateLookup.clear()
+ self.masterVRT_InputStacks = listOfInputStacks
+ self.masterVRT_SourceBandTemplates.clear()
+ #dates = set()
+ #for s in listOfInputStacks:
+ # for d in s.dates():
+ # dates.add(d)
+ #dates = sorted(list(dates))
+
+ #create a LUT to get the stack indices for a related date (not each stack might contain a band for each date)
+
+ for stackIndex, s in enumerate(listOfInputStacks):
+ for bandIndex, bandDate in enumerate(s.dates()):
+ if bandDate not in self.masterVRT_DateLookup.keys():
+ self.masterVRT_DateLookup[bandDate] = []
+ self.masterVRT_DateLookup[bandDate].append((stackIndex, bandIndex))
+
+ #create VRT Template XML
+ VRT = VRTRaster()
+ wavelength = []
+ for stackIndex, stack in enumerate(listOfInputStacks):
+ assert isinstance(stack, InputStackInfo)
+ vrtBand = VRTRasterBand()
+ vrtBand.setName(stack.outputBandName)
+ vrtSrc = VRTRasterInputSourceBand(stack.path, 0)
+ vrtBand.addSource(vrtSrc)
+ wavelength.append(stack.wavelength())
+ VRT.addVirtualBand(vrtBand)
+
+
+
+ pathVSITmp = '/vsimem/temp.vrt'
+ dsVRT = VRT.saveVRT(pathVSITmp)
+ dsVRT.SetMetadataItem('acquisition date', 'XML_REPLACE_DATE')
+
+ if None not in wavelength:
+ dsVRT.SetMetadataItem('wavelength', ','.join(str(wl) for wl in wavelength))
+ dsVRT.SetMetadataItem('wavelength units', 'Nanometers')
+
+ for stackIndex, stack in enumerate(listOfInputStacks):
+ band = dsVRT.GetRasterBand(stackIndex+1)
+ assert isinstance(band, gdal.Band)
+ assert isinstance(stack, InputStackInfo)
+ if stack.colorTable:
+ band.SetColorTable(stack.colorTable)
+ if stack.classNames:
+ band.SetCategoryNames(stack.classNames)
+
+ dsVRT.FlushCache()
+ drv = dsVRT.GetDriver()
+ masterVRT_XML = read_vsimem(pathVSITmp).decode('utf-8')
+ drv.Delete(pathVSITmp)
+ outputVRTs = []
+
+
+ eTree = ElementTree.fromstring(masterVRT_XML)
+ for iBand, elemBand in enumerate(eTree.findall('VRTRasterBand')):
+ sourceElements = elemBand.findall('ComplexSource') + elemBand.findall('SimpleSource')
+ assert len(sourceElements) == 1
+ self.masterVRT_SourceBandTemplates[iBand] = copy.deepcopy(sourceElements[0])
+ elemBand.remove(sourceElements[0])
+
+ for date in dates:
+ assert isinstance(date, np.datetime64)
+ path = self.createVRTUri(date)
+ outputDescription = OutputVRTDescription(path, date)
+ outputVRTs.append(outputDescription)
+
+ self.masterVRT_XML = eTree
+
+
+ self.beginInsertRows(QModelIndex(), 0, len(outputVRTs)-1)
+ self.mOutputImages = outputVRTs[:]
+ self.endInsertRows()
+
+ def setOutputDir(self, path:str):
+ self.mOutputDir = path
+ self.updateOutputURIs()
+
+ def setOutputPrefix(self, basename:str):
+ self.mOutputPrefix = basename
+ self.updateOutputURIs()
+
+ def updateOutputURIs(self):
+ c = self.mColumnNames.index(self.cn_uri)
+ ul = self.createIndex(0, c)
+ lr = self.createIndex(self.rowCount()-1, c)
+
+ for outputVRT in self:
+ assert isinstance(outputVRT, OutputVRTDescription)
+ outputVRT.setPath(self.createVRTUri(outputVRT.mDate))
+ self.dataChanged.emit(ul, lr)
+
+
+ def __len__(self):
+ return len(self.mOutputImages)
+
+ def __iter__(self):
+ return iter(self.mOutputImages)
+
+ def rowCount(self, parent=None) -> int:
+ return len(self.mOutputImages)
+
+ def columnCount(self, parent=None) -> int:
+ return len(self.mColumnNames)
+
+ def columnName(self, i) -> str:
+ if isinstance(i, QModelIndex):
+ i = i.column()
+ return self.mColumnNames[i]
+
+ def columnIndex(self, columnName:str)-> QModelIndex:
+ c = self.mColumnNames.index(columnName)
+ return self.createIndex(0, c)
+
+ def index2vrt(self, index:QModelIndex) -> OutputVRTDescription:
+ return self.mOutputImages[index.row()]
+
+ def vrt2index(self, vrt:OutputVRTDescription) -> QModelIndex:
+ i = self.mOutputImages[vrt]
+ return self.createIndex(i, 0, vrt)
+
+ def data(self, index: QModelIndex, role: int):
+
+ if not index.isValid():
+ return None
+
+ cname = self.columnName(index)
+ vrt = self.index2vrt(index)
+ if role in [Qt.DisplayRole, Qt.ToolTipRole]:
+ if cname == self.cn_uri:
+ return vrt.mPath
+ if cname == self.cn_date:
+ return str(vrt.mDate)
+
+ def vrtXML(self, outputDefinition:OutputVRTDescription, asElementTree=False) -> str:
+ """
+ Create the VRT XML related to an outputDefinition
+ :param outputDefinition:
+ :return: str
+ """
+
+ # re.search(tmpXml, '<MDI key='>')
+
+ # xml = copy.deepcopy(eTree)
+ if self.masterVRT_XML is None:
+ return None
+ #xmlTree = ElementTree.fromstring(self.masterVRT_XML)
+ xmlTree = copy.deepcopy(self.masterVRT_XML)
+
+ # set metadata
+ for elem in xmlTree.findall('Metadata/MDI'):
+ if elem.attrib['key'] == 'acquisition date':
+ elem.text = str(outputDefinition.mDate)
+
+ # insert required rasterbands
+ requiredBands = self.masterVRT_DateLookup[outputDefinition.mDate]
+
+ xmlVRTBands = xmlTree.findall('VRTRasterBand')
+
+ for t in requiredBands:
+ stackIndex, stackBandIndex = t
+
+ stackSourceXMLTemplate = copy.deepcopy(self.masterVRT_SourceBandTemplates[stackIndex])
+ stackSourceXMLTemplate.find('SourceBand').text = str(stackBandIndex+1)
+ xmlVRTBands[stackIndex].append(stackSourceXMLTemplate)
+
+ if asElementTree:
+ return xmlTree
+ else:
+ return ElementTree.tostring(xmlTree).decode('utf-8')
+
+
+
+
+
+class StackedBandInputDialog(QDialog, loadUI('stackedinputdatadialog.ui')):
+
+ def __init__(self, parent=None):
+
+ super(StackedBandInputDialog, self).__init__(parent=parent)
+ self.setupUi(self)
+ self.setWindowTitle('Stacked Time Series Data Input')
+ self.mWrittenFiles = []
+
+ self.tableModelInputStacks = InputStackTableModel()
+ self.tableModelInputStacks.rowsInserted.connect(self.updateOutputs)
+ self.tableModelInputStacks.dataChanged.connect(self.updateOutputs)
+ self.tableModelInputStacks.rowsRemoved.connect(self.updateOutputs)
+ self.tableModelInputStacks.rowsInserted.connect(self.updateInputInfo)
+ self.tableModelInputStacks.rowsRemoved.connect(self.updateInputInfo)
+ self.tableViewSourceStacks.setModel(self.tableModelInputStacks)
+
+ self.tableModelOutputImages = OutputImageModel()
+ self.tableModelOutputImages.rowsInserted.connect(self.updateOutputInfo)
+ self.tableModelOutputImages.rowsRemoved.connect(self.updateOutputInfo)
+ self.tableModelOutputImages.dataChanged.connect(self.updateOutputInfo)
+ self.tableViewOutputImages.setModel(self.tableModelOutputImages)
+ self.tableViewOutputImages.horizontalHeader().setSectionResizeMode(QHeaderView.ResizeToContents)
+
+ self.buttonGroupDateMode.buttonClicked.connect(self.updateOutputs)
+ self.buttonGroupOutputLocation.buttonClicked.connect(self.updateOutputs)
+
+ self.cbOpenInQGIS.setEnabled(isinstance(qgis.utils.iface, QgisInterface))
+ self.tbFilePrefix.textChanged.connect(self.tableModelOutputImages.setOutputPrefix)
+ self.tbFilePrefix.setText('img')
+
+ self.fileWidgetOutputDir.setStorageMode(QgsFileWidget.GetDirectory)
+ self.fileWidgetOutputDir.fileChanged.connect(self.tableModelOutputImages.setOutputDir)
+
+ sm = self.tableViewSourceStacks.selectionModel()
+ assert isinstance(sm, QItemSelectionModel)
+ sm.selectionChanged.connect(self.onSourceStackSelectionChanged)
+ self.onSourceStackSelectionChanged([],[])
+
+ sm = self.tableViewOutputImages.selectionModel()
+ assert isinstance(sm, QItemSelectionModel)
+ sm.selectionChanged.connect(self.onOutputImageSelectionChanged)
+
+ self.initActions()
+
+ def writtenFiles(self):
+ """
+ Returns the files written after pressing the "Save" button.
+ :return: [list-of-written-file-paths]
+ """
+ return self.mWrittenFiles[:]
+
+ def updateOutputs(self, *args):
+ """
+ Updates the output file information
+ """
+ self.tableModelOutputImages.clearOutputs()
+ inputStacks = self.tableModelInputStacks.mStackImages
+ datesTotal, datesIntersection = self.tableModelInputStacks.dateInfo()
+ if self.rbDatesAll.isChecked():
+ self.tableModelOutputImages.setMultiStackSources(inputStacks, datesTotal)
+ elif self.rbDatesIntersection.isChecked():
+ self.tableModelOutputImages.setMultiStackSources(inputStacks, datesIntersection)
+
+ if self.rbSaveInMemory.isChecked():
+ self.tableModelOutputImages.setOutputDir('/vsimem/')
+ elif self.rbSaveInDirectory.isChecked():
+ self.tableModelOutputImages.setOutputDir(self.fileWidgetOutputDir.filePath())
+
+ def updateInputInfo(self):
+ """
+ Updates the input file information
+ """
+
+ n = len(self.tableModelInputStacks)
+ datesTotal, datesInCommon = self.tableModelInputStacks.dateInfo()
+ info = None
+ if n > 0:
+ nAll = len(datesTotal)
+ nInt = len(datesInCommon)
+ info = '{} Input Images with {} dates in total, {} in intersection'.format(n, nAll, nInt)
+
+ self.tbInfoInputImages.setText(info)
+
+ def updateOutputInfo(self):
+
+ n = len(self.tableModelOutputImages)
+ info = None
+ if n > 0:
+ nb = len(self.tableModelOutputImages.masterVRT_InputStacks)
+ info = '{} output images with {} bands to {}'.format(n, nb, self.tableModelOutputImages.mOutputDir)
+ self.buttonBox.button(QDialogButtonBox.Save).setEnabled(n > 0)
+ self.tbInfoOutputImages.setText(info)
+
+ def initActions(self):
+ """
+ Initializes QActions and what they trigger.
+ """
+
+ self.actionAddSourceStack.triggered.connect(self.onAddSource)
+ self.actionRemoveSourceStack.triggered.connect(self.onRemoveSources)
+
+ self.btnAddSourceStack.setDefaultAction(self.actionAddSourceStack)
+ self.btnRemoveSourceStack.setDefaultAction(self.actionRemoveSourceStack)
+
+ self.buttonBox.button(QDialogButtonBox.Save).clicked.connect(self.accept)
+ self.buttonBox.button(QDialogButtonBox.Cancel).clicked.connect(self.close)
+
+ def onAddSource(self, *args):
+ """
+ Reacts on new added datasets
+ """
+ from timeseriesviewer import SETTINGS
+ defDir = SETTINGS.value('DIR_FILESEARCH')
+ filters = QgsProviderRegistry.instance().fileVectorFilters()
+ files, filter = QFileDialog.getOpenFileNames(directory=defDir, filter=filters)
+
+ if len(files) > 0:
+ self.tableModelInputStacks.insertSources(files)
+ s = ""
+
+
+ def addSources(self, paths):
+ """
+ Adds new datasources
+ :param paths: [list-of-new-datasources]
+ :return:
+ """
+ self.tableModelInputStacks.insertSources(paths)
+
+ def onRemoveSources(self, *args):
+
+ model = self.tableViewSourceStacks.selectionModel()
+ assert isinstance(model, QItemSelectionModel)
+
+ infos = [self.tableModelInputStacks.index2info(idx) for idx in model.selectedRows()]
+ self.tableModelInputStacks.removeSources(infos)
+
+ def onSourceStackSelectionChanged(self, selected, deselected):
+
+ self.actionRemoveSourceStack.setEnabled(len(selected) > 0)
+
+
+ def onOutputImageSelectionChanged(self, selected, deselected):
+
+ if len(selected) > 0:
+ idx = selected.indexes()[0]
+
+ vrtOutput = self.tableModelOutputImages.index2vrt(idx)
+ assert isinstance(vrtOutput, OutputVRTDescription)
+ xml = self.tableModelOutputImages.vrtXML(vrtOutput)
+ self.tbXMLPreview.setPlainText(xml)
+ else:
+ self.tbXMLPreview.setPlainText(None)
+ s = ""
+
+
+ def saveImages(self):
+ """
+ Write the VRT images
+ :return: [list-of-written-file-paths]
+ """
+
+
+ nTotal = len(self.tableModelOutputImages)
+ if nTotal == 0:
+ return
+
+ writtenFiles = []
+ self.progressBar.setValue(0)
+ from timeseriesviewer.virtualrasters import write_vsimem, read_vsimem
+ for i, outVRT in enumerate(self.tableModelOutputImages):
+ assert isinstance(outVRT, OutputVRTDescription)
+ xml = self.tableModelOutputImages.vrtXML(outVRT)
+
+ if outVRT.mPath.startswith('/vsimem/'):
+ write_vsimem(outVRT.mPath, xml)
+ else:
+ f = open(outVRT.mPath, 'w', encoding='utf-8')
+ f.write(xml)
+ f.flush()
+ f.close()
+
+ writtenFiles.append(outVRT.mPath)
+
+ self.progressBar.setValue(int(100. * i / nTotal))
+
+ QTimer.singleShot(500, lambda: self.progressBar.setValue(0))
+
+ if self.cbOpenInQGIS.isEnabled() and self.cbOpenInQGIS.isChecked():
+ mapLayers = [QgsRasterLayer(p) for p in writtenFiles]
+ QgsProject.instance().addMapLayers(mapLayers, addToLegend=True)
+ self.mWrittenFiles.extend(writtenFiles)
+ return writtenFiles
\ No newline at end of file
diff --git a/timeseriesviewer/temporalprofiles3dGL.py b/timeseriesviewer/temporalprofiles3dGL.py
index bb6dcd1562627a103cbbe507439794b64105c1e9..a9068258794d28d370078a3052003ea6c3fe7ce2 100644
--- a/timeseriesviewer/temporalprofiles3dGL.py
+++ b/timeseriesviewer/temporalprofiles3dGL.py
@@ -750,7 +750,6 @@ if __name__ == '__main__':
app = initQgisApplication()
from timeseriesviewer.utils import *
- from timeseriesviewer.tests import TestObjects
w = ViewWidget3D()
w.show()
diff --git a/timeseriesviewer/timeseries.py b/timeseriesviewer/timeseries.py
index 577a6c41e407486831e77333d618d7db06a11e66..c02c95dd798c57a44824a00a5c08a6009c63c21e 100644
--- a/timeseriesviewer/timeseries.py
+++ b/timeseriesviewer/timeseries.py
@@ -87,6 +87,9 @@ def getDS(pathOrDataset):
class SensorInstrument(QObject):
+ """
+ Describes a Sensor Configuration
+ """
SensorNameSettingsPrefix = 'SensorName.'
sigNameChanged = pyqtSignal(str)
@@ -98,9 +101,7 @@ class SensorInstrument(QObject):
'swIR1':1650,
'swIR2':2150
})
- """
- Describes a Sensor Configuration
- """
+
def __init__(self, pathImg, sensor_name=None):
super(SensorInstrument, self).__init__()
@@ -126,9 +127,9 @@ class SensorInstrument(QObject):
self.wavelengths = np.asarray(wl)
- self._id = '{}b{}m'.format(self.nb, self.px_size_x)
+ self.mId = '{}b{}m'.format(self.nb, self.px_size_x)
if wl is not None:
- self._id += ';'.join([str(w) for w in self.wavelengths])+str(self.wavelengthUnits)
+ self.mId += ';'.join([str(w) for w in self.wavelengths]) + str(self.wavelengthUnits)
if sensor_name is None:
@@ -141,10 +142,10 @@ class SensorInstrument(QObject):
def id(self):
- return self._id
+ return self.mId
def _sensorSettingsKey(self):
- return SensorInstrument.SensorNameSettingsPrefix+self._id
+ return SensorInstrument.SensorNameSettingsPrefix+self.mId
def setName(self, name):
self._name = name
@@ -182,6 +183,10 @@ class SensorInstrument(QObject):
def wavelengthsDefined(self):
+ """
+ Returns if wavelenght and wavelength units are defined
+ :return: bool
+ """
return self.wavelengths is not None and \
self.wavelengthUnits is not None
@@ -199,6 +204,10 @@ class SensorInstrument(QObject):
return str(self.__class__) +' ' + self.name()
def getDescription(self):
+ """
+ Returns a human-readable description
+ :return: str
+ """
info = []
info.append(self.name())
info.append('{} Bands'.format(self.nb))
@@ -213,42 +222,45 @@ class SensorInstrument(QObject):
return '\n'.join(info)
-def verifyInputImage(path, vrtInspection=''):
- if path is None or not isinstance(path, str):
- return None
- ds = gdal.Open(path)
+def verifyInputImage(datasource):
+ """
+ Checks if an image source can be uses as TimeSeriesDatum, i.e. if it can be read by gdal.Open() and
+ if we can extract an observation date as numpy.datetime64.
+ :param datasource: str with data source uri or gdal.Dataset
+ :return: bool
+ """
- if not ds:
- #logger.error('{}GDAL unable to open: '.format(vrtInspection, path))
+ if datasource is None:
+ return None
+ if isinstance(datasource, str):
+ datasource = gdal.Open(datasource)
+ if not isinstance(datasource, gdal.Dataset):
return False
- if ds.RasterCount == 0 and len(ds.GetSubDatasets()) > 0:
+ if datasource.RasterCount == 0 and len(datasource.GetSubDatasets()) > 0:
#logger.error('Can not open container {}.\nPlease specify a subdataset'.format(path))
return False
- if ds.GetDriver().ShortName == 'VRT':
- vrtInspection = 'VRT Inspection {}\n'.format(path)
- nextFiles = set(ds.GetFileList()) - set([path])
- validSrc = [verifyInputImage(p, vrtInspection=vrtInspection) for p in nextFiles]
- if not all(validSrc):
- return False
+ if datasource.GetDriver().ShortName == 'VRT':
+ files = datasource.GetFileList()
+ if len(files) > 0:
+ for f in files:
+ subDS = gdal.Open(f)
+ if not isinstance(subDS, gdal.Dataset):
+ return False
from timeseriesviewer.dateparser import parseDateFromDataSet
- date = parseDateFromDataSet(ds)
+ date = parseDateFromDataSet(datasource)
if date is None:
return False
return True
-def pixel2coord(gt, x, y):
- """Returns global coordinates from pixel x, y coords"""
- """https://scriptndebug.wordpress.com/2014/11/24/latitudelongitude-of-each-pixel-using-python-and-gdal/"""
- xoff, a, b, yoff, d, e = gt
- xp = a * x + b * y + xoff
- yp = d * x + e * y + yoff
- return (xp, yp)
-
class TimeSeriesDatum(QObject):
+ """
+ The core entity of a TimeSeries. Contains a data source, an observation date (numpy.datetime64) and
+ the SensorInstrument the data source is linked to.
+ """
@staticmethod
def createFromPath(path):
"""
@@ -267,11 +279,6 @@ class TimeSeriesDatum(QObject):
return tsd
-
-
- """
- Collects all data sets related to one sensor
- """
sigVisibilityChanged = pyqtSignal(bool)
sigRemoveMe = pyqtSignal()
@@ -296,9 +303,9 @@ class TimeSeriesDatum(QObject):
gt = ds.GetGeoTransform()
-
- UL = QgsPointXY(*pixel2coord(gt, 0, 0))
- LR = QgsPointXY(*pixel2coord(gt, self.ns, self.nl))
+ from timeseriesviewer.utils import px2geo
+ UL = QgsPointXY(*px2geo(QPoint(0,0), gt, pxCenter=False))
+ LR = QgsPointXY(*px2geo(QPoint(self.ns+1, self.nl+1), gt, pxCenter=False))
self._spatialExtent = SpatialExtent(self.crs, UL, LR)
self.srs_wkt = str(self.crs.toWkt())
@@ -307,27 +314,44 @@ class TimeSeriesDatum(QObject):
self.mVisibility = True
- def rank(self):
- return self.timeSeries.index(self)
-
def setVisibility(self, b):
+ """
+ Sets the visibility of the TimeSeriesDatum, i.e. whether linked MapCanvases will be shown to the user
+ :param b: bool
+ """
old = self.mVisibility
self.mVisibility = b
if old != self.mVisibility:
self.sigVisibilityChanged.emit(b)
def isVisible(self):
+ """
+ Returns whether the TimeSeriesDatum is visible as MapCanvas
+ :return: bool
+ """
return self.mVisibility
def getDate(self):
+ """
+ Returns the observation date
+ :return: numpy.datetime64
+ """
return np.datetime64(self.date)
def getSpatialReference(self):
+ """
+ Returns the QgsCoordinateReferenceSystem of the data source.
+ :return: QgsCoordinateReferenceSystem
+ """
return self.crs
def spatialExtent(self):
+ """
+ Returns the SpatialExtent of the data source
+ :return: SpatialExtent
+ """
return self._spatialExtent
def __repr__(self):
@@ -354,7 +378,10 @@ class TimeSeriesTableView(QTableView):
super(TimeSeriesTableView, self).__init__(parent)
def contextMenuEvent(self, event):
-
+ """
+ Creates and shows an QMenu
+ :param event:
+ """
menu = QMenu(self)
a = menu.addAction('Copy value(s)')
a.triggered.connect(self.onCopyValues)
@@ -365,6 +392,10 @@ class TimeSeriesTableView(QTableView):
menu.popup(QCursor.pos())
def onSetCheckState(self, checkState):
+ """
+ Sets a ChecState to all selected rows
+ :param checkState: Qt.CheckState
+ """
indices = self.selectionModel().selectedIndexes()
rows = sorted(list(set([i.row() for i in indices])))
model = self.model()
@@ -374,6 +405,9 @@ class TimeSeriesTableView(QTableView):
model.setData(idx, checkState, Qt.CheckStateRole)
def onCopyValues(self):
+ """
+ Copies selected cell values to the clipboard
+ """
indices = self.selectionModel().selectedIndexes()
model = self.model()
if isinstance(model, TimeSeriesTableModel):
@@ -390,12 +424,13 @@ class TimeSeriesTableView(QTableView):
QApplication.clipboard().setText(info)
s = ""
- def selectSelectedObservations(b):
- assert isinstance(b, bool)
class TimeSeriesDockUI(QgsDockWidget, loadUI('timeseriesdock.ui')):
+ """
+ QgsDockWidget that shows the TimeSeries
+ """
def __init__(self, parent=None):
super(TimeSeriesDockUI, self).__init__(parent)
self.setupUi(self)
@@ -414,6 +449,9 @@ class TimeSeriesDockUI(QgsDockWidget, loadUI('timeseriesdock.ui')):
self.setTimeSeries(None)
def setStatus(self):
+ """
+ Updates the status of the TimeSeries
+ """
from timeseriesviewer.timeseries import TimeSeries
if isinstance(self.TS, TimeSeries):
ndates = len(self.TS)
@@ -423,7 +461,13 @@ class TimeSeriesDockUI(QgsDockWidget, loadUI('timeseriesdock.ui')):
msg += ', {} to {}'.format(str(self.TS[0].date), str(self.TS[-1].date))
self.progressInfo.setText(msg)
- def setProgressInfo(self, nDone, nMax, message=None):
+ def setProgressInfo(self, nDone:int, nMax:int, message=None):
+ """
+ Sets the progress bar of the TimeSeriesDockUI
+ :param nDone: number of added data sources
+ :param nMax: total number of data source to be added
+ :param message: error / other kind of info message
+ """
if self.progressBar.maximum() != nMax:
self.progressBar.setMaximum(nMax)
self.progressBar.setValue(nDone)
@@ -433,14 +477,25 @@ class TimeSeriesDockUI(QgsDockWidget, loadUI('timeseriesdock.ui')):
QTimer.singleShot(3000, lambda: self.setStatus())
def onSelectionChanged(self, *args):
+ """
+ Slot to react on user-driven changes of the selected TimeSeriesDatum rows.
+ """
self.btnRemoveTSD.setEnabled(self.SM is not None and len(self.SM.selectedRows()) > 0)
def selectedTimeSeriesDates(self):
+ """
+ Returns the TimeSeriesDatum selected by a user.
+ :return: [list-of-TimeSeriesDatum]
+ """
if self.SM is not None:
return [self.TSM.data(idx, Qt.UserRole) for idx in self.SM.selectedRows()]
return []
def setTimeSeries(self, TS):
+ """
+ Sets the TimeSeries to be shown in the TimeSeriesDockUI
+ :param TS: TimeSeries
+ """
from timeseriesviewer.timeseries import TimeSeries
self.TS = TS
self.TSM = None
@@ -461,6 +516,9 @@ class TimeSeriesDockUI(QgsDockWidget, loadUI('timeseriesdock.ui')):
class TimeSeries(QObject):
+ """
+ The sorted list of data sources that specify the time series
+ """
sigTimeSeriesDatesAdded = pyqtSignal(list)
sigTimeSeriesDatesRemoved = pyqtSignal(list)
@@ -495,9 +553,18 @@ class TimeSeries(QObject):
_sep = ';'
def sensors(self):
+ """
+ Returns the list of sensors derived from the TimeSeries data soures
+ :return: [list-of-SensorInstruments]
+ """
return list(self.Sensors.keys())
def loadFromFile(self, path, n_max=None):
+ """
+ Loads a CSV file with source images of a TimeSeries
+ :param path: str, Path of CSV file
+ :param n_max: optional, maximum number of files to load
+ """
images = []
masks = []
@@ -522,6 +589,11 @@ class TimeSeries(QObject):
def saveToFile(self, path):
+ """
+ Saves the TimeSeries sources into a CSV file
+ :param path: str, path of CSV file
+ :return: path of CSV file
+ """
if path is None or len(path) == 0:
return None
@@ -541,7 +613,12 @@ class TimeSeries(QObject):
messageLog('Time series source images written to {}'.format(path))
return path
+
def getPixelSizes(self):
+ """
+ Returns the pixel sizes of all SensorInstruments
+ :return: [list-of-QgsRectangles]
+ """
r = []
for sensor in self.Sensors.keys():
@@ -549,7 +626,13 @@ class TimeSeries(QObject):
return r
return None
+
def getMaxSpatialExtent(self, crs=None):
+ """
+ Returns the maximum SpatialExtent of all images of the TimeSeries
+ :param crs: QgsCoordinateSystem to express the SpatialExtent coordinates.
+ :return:
+ """
if len(self.data) == 0:
return None
@@ -565,23 +648,31 @@ class TimeSeries(QObject):
extent = extent.toCrs(crs)
return extent
-
- def tsdFromPath(self, path):
- for tsd in self.data:
- if tsd.pathImg == path:
- return tsd
- return False
-
def getObservationDates(self):
+ """
+ Returns the observations dates of ecach TimeSeries data source
+ :return: [list-of-numpy.datetime64]
+ """
return [tsd.getDate() for tsd in self.data]
def getTSD(self, pathOfInterest):
+ """
+ Returns the TimeSeriesDatum related to an image source
+ :param pathOfInterest: str, image source uri
+ :return: TimeSeriesDatum
+ """
for tsd in self.data:
if tsd.pathImg == pathOfInterest:
return tsd
return None
def getTSDs(self, dateOfInterest=None, sensorOfInterest=None):
+ """
+ Returns a list of TimeSeriesDatum of the TimeSeries. By default all TimeSeriesDatum will be returned.
+ :param dateOfInterest: numpy.datetime64 to return the TimeSeriesDatum for
+ :param sensorOfInterest: SensorInstrument of interest to return the [list-of-TimeSeriesDatum] for.
+ :return: [list-of-TimeSeriesDatum]
+ """
tsds = self.data[:]
if dateOfInterest:
tsds = [tsd for tsd in tsds if tsd.getDate() == dateOfInterest]
@@ -590,16 +681,23 @@ class TimeSeries(QObject):
return tsds
def clear(self):
+ """
+ Removes all data sources from the TimeSeries (which will be empty after calling this routine).
+ """
self.removeDates(self[:])
def removeDates(self, TSDs):
+ """
+ Removes a list of TimeSeriesDatum
+ :param TSDs: [list-of-TimeSeriesDatum]
+ """
removed = list()
for TSD in TSDs:
assert type(TSD) is TimeSeriesDatum
self.data.remove(TSD)
- TSD.timeSeries = None
+ TSD.createTimeSeries = None
removed.append(TSD)
S = TSD.sensor
@@ -612,6 +710,10 @@ class TimeSeries(QObject):
def addTimeSeriesDates(self, timeSeriesDates):
+ """
+ Adds a list of TimeSeriesDatum
+ :param timeSeriesDates: [list-of-TimeSeriesDatum]
+ """
assert isinstance(timeSeriesDates, list)
added = list()
for TSD in timeSeriesDates:
@@ -649,6 +751,10 @@ class TimeSeries(QObject):
def addFiles(self, files):
+ """
+ Adds new data sources to the TimeSeries
+ :param files: [list-of-str]
+ """
if isinstance(files, str):
files = [files]
assert isinstance(files, list)
@@ -656,7 +762,7 @@ class TimeSeries(QObject):
nMax = len(files)
nDone = 0
- self.sigLoadingProgress.emit(0,nMax, 'Start loading {} files...'.format(nMax))
+ self.sigLoadingProgress.emit(0, nMax, 'Start loading {} files...'.format(nMax))
for i, file in enumerate(files):
t0 = np.datetime64('now')
diff --git a/timeseriesviewer/ui/resources.qrc b/timeseriesviewer/ui/resources.qrc
index 71a4f7f6395a9ccaa547565d128c4c4c44c5ea9b..fdf14936064b520d963ec187646c8c65358911a2 100644
--- a/timeseriesviewer/ui/resources.qrc
+++ b/timeseriesviewer/ui/resources.qrc
@@ -1,6 +1,5 @@
<RCC>
<qresource prefix="timeseriesviewer">
- <file>../../ABOUT.txt</file>
<file>icons/mIconTemporalProfileRefresh.svg</file>
<file>icons/profile2speclib_auto.svg</file>
<file>icons/mActionAddOgrLayer.svg</file>
diff --git a/timeseriesviewer/ui/stackedinputdatadialog.ui b/timeseriesviewer/ui/stackedinputdatadialog.ui
new file mode 100644
index 0000000000000000000000000000000000000000..58a7f330e303d36f53ff0b9b0b714be518347107
--- /dev/null
+++ b/timeseriesviewer/ui/stackedinputdatadialog.ui
@@ -0,0 +1,459 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>Dialog</class>
+ <widget class="QDialog" name="Dialog">
+ <property name="geometry">
+ <rect>
+ <x>0</x>
+ <y>0</y>
+ <width>569</width>
+ <height>573</height>
+ </rect>
+ </property>
+ <property name="windowTitle">
+ <string>Dialog</string>
+ </property>
+ <layout class="QVBoxLayout" name="verticalLayout">
+ <item>
+ <widget class="QgsCollapsibleGroupBox" name="gbInputStacks">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Preferred" vsizetype="Preferred">
+ <horstretch>0</horstretch>
+ <verstretch>1</verstretch>
+ </sizepolicy>
+ </property>
+ <property name="minimumSize">
+ <size>
+ <width>0</width>
+ <height>0</height>
+ </size>
+ </property>
+ <property name="title">
+ <string>Input Images (each band = one temporal observation)</string>
+ </property>
+ <property name="collapsed">
+ <bool>false</bool>
+ </property>
+ <layout class="QVBoxLayout" name="verticalLayout_4">
+ <item>
+ <layout class="QHBoxLayout" name="horizontalLayout">
+ <property name="topMargin">
+ <number>0</number>
+ </property>
+ <item>
+ <widget class="QToolButton" name="btnAddSourceStack">
+ <property name="text">
+ <string>...</string>
+ </property>
+ <property name="autoRaise">
+ <bool>true</bool>
+ </property>
+ </widget>
+ </item>
+ <item>
+ <widget class="QToolButton" name="btnRemoveSourceStack">
+ <property name="text">
+ <string>...</string>
+ </property>
+ <property name="autoRaise">
+ <bool>true</bool>
+ </property>
+ </widget>
+ </item>
+ <item>
+ <spacer name="horizontalSpacer_2">
+ <property name="orientation">
+ <enum>Qt::Horizontal</enum>
+ </property>
+ <property name="sizeHint" stdset="0">
+ <size>
+ <width>40</width>
+ <height>20</height>
+ </size>
+ </property>
+ </spacer>
+ </item>
+ </layout>
+ </item>
+ <item>
+ <widget class="QTableView" name="tableViewSourceStacks"/>
+ </item>
+ <item>
+ <widget class="QLineEdit" name="tbInfoInputImages">
+ <property name="autoFillBackground">
+ <bool>false</bool>
+ </property>
+ <property name="styleSheet">
+ <string notr="true">background-color: transparent;</string>
+ </property>
+ <property name="text">
+ <string/>
+ </property>
+ <property name="frame">
+ <bool>false</bool>
+ </property>
+ <property name="readOnly">
+ <bool>true</bool>
+ </property>
+ <property name="placeholderText">
+ <string>Add images with stacked observations to the list of input images</string>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </widget>
+ </item>
+ <item>
+ <widget class="QgsCollapsibleGroupBox" name="mGroupBox">
+ <property name="title">
+ <string>Output options</string>
+ </property>
+ <property name="collapsed">
+ <bool>false</bool>
+ </property>
+ <layout class="QGridLayout" name="gridLayout_2">
+ <item row="0" column="3">
+ <widget class="QLineEdit" name="tbFilePrefix">
+ <property name="text">
+ <string>timestack</string>
+ </property>
+ </widget>
+ </item>
+ <item row="0" column="0">
+ <widget class="QLabel" name="label">
+ <property name="text">
+ <string>Dates from</string>
+ </property>
+ </widget>
+ </item>
+ <item row="1" column="0">
+ <widget class="QLabel" name="label_5">
+ <property name="text">
+ <string>Save images to</string>
+ </property>
+ </widget>
+ </item>
+ <item row="1" column="1" colspan="3">
+ <layout class="QHBoxLayout" name="horizontalLayout_4">
+ <property name="topMargin">
+ <number>0</number>
+ </property>
+ <item>
+ <widget class="QRadioButton" name="rbSaveInMemory">
+ <property name="text">
+ <string>memory</string>
+ </property>
+ <property name="checked">
+ <bool>true</bool>
+ </property>
+ <attribute name="buttonGroup">
+ <string notr="true">buttonGroupOutputLocation</string>
+ </attribute>
+ </widget>
+ </item>
+ <item>
+ <widget class="QRadioButton" name="rbSaveInDirectory">
+ <property name="text">
+ <string>directory</string>
+ </property>
+ <attribute name="buttonGroup">
+ <string notr="true">buttonGroupOutputLocation</string>
+ </attribute>
+ </widget>
+ </item>
+ <item>
+ <widget class="QgsFileWidget" name="fileWidgetOutputDir">
+ <property name="enabled">
+ <bool>false</bool>
+ </property>
+ <property name="dialogTitle">
+ <string>Select output directory</string>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </item>
+ <item row="0" column="1">
+ <layout class="QHBoxLayout" name="horizontalLayout_2">
+ <item>
+ <widget class="QRadioButton" name="rbDatesIntersection">
+ <property name="text">
+ <string>Intersection</string>
+ </property>
+ <property name="checked">
+ <bool>true</bool>
+ </property>
+ <attribute name="buttonGroup">
+ <string notr="true">buttonGroupDateMode</string>
+ </attribute>
+ </widget>
+ </item>
+ <item>
+ <widget class="QRadioButton" name="rbDatesAll">
+ <property name="text">
+ <string>All Images</string>
+ </property>
+ <attribute name="buttonGroup">
+ <string notr="true">buttonGroupDateMode</string>
+ </attribute>
+ </widget>
+ </item>
+ <item>
+ <spacer name="horizontalSpacer_4">
+ <property name="orientation">
+ <enum>Qt::Horizontal</enum>
+ </property>
+ <property name="sizeHint" stdset="0">
+ <size>
+ <width>40</width>
+ <height>20</height>
+ </size>
+ </property>
+ </spacer>
+ </item>
+ </layout>
+ </item>
+ <item row="0" column="2">
+ <widget class="QLabel" name="label_6">
+ <property name="text">
+ <string>Prefix</string>
+ </property>
+ </widget>
+ </item>
+ <item row="2" column="1">
+ <widget class="QCheckBox" name="cbOpenInQGIS">
+ <property name="text">
+ <string>Open in QGIS</string>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </widget>
+ </item>
+ <item>
+ <widget class="QgsCollapsibleGroupBox" name="gbOutputImages">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Preferred" vsizetype="Preferred">
+ <horstretch>0</horstretch>
+ <verstretch>2</verstretch>
+ </sizepolicy>
+ </property>
+ <property name="title">
+ <string>Output Images (each image = one temporal observation)</string>
+ </property>
+ <property name="collapsed">
+ <bool>false</bool>
+ </property>
+ <layout class="QVBoxLayout" name="verticalLayout_3">
+ <item>
+ <widget class="QSplitter" name="splitter">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Expanding" vsizetype="Preferred">
+ <horstretch>0</horstretch>
+ <verstretch>0</verstretch>
+ </sizepolicy>
+ </property>
+ <property name="orientation">
+ <enum>Qt::Horizontal</enum>
+ </property>
+ <property name="childrenCollapsible">
+ <bool>true</bool>
+ </property>
+ <widget class="QTableView" name="tableViewOutputImages">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Expanding" vsizetype="Expanding">
+ <horstretch>0</horstretch>
+ <verstretch>0</verstretch>
+ </sizepolicy>
+ </property>
+ </widget>
+ <widget class="QPlainTextEdit" name="tbXMLPreview">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Expanding" vsizetype="Expanding">
+ <horstretch>2</horstretch>
+ <verstretch>0</verstretch>
+ </sizepolicy>
+ </property>
+ <property name="verticalScrollBarPolicy">
+ <enum>Qt::ScrollBarAlwaysOn</enum>
+ </property>
+ <property name="horizontalScrollBarPolicy">
+ <enum>Qt::ScrollBarAlwaysOn</enum>
+ </property>
+ <property name="sizeAdjustPolicy">
+ <enum>QAbstractScrollArea::AdjustToContents</enum>
+ </property>
+ <property name="lineWrapMode">
+ <enum>QPlainTextEdit::NoWrap</enum>
+ </property>
+ <property name="readOnly">
+ <bool>true</bool>
+ </property>
+ <property name="placeholderText">
+ <string>Select output image to show its definition here</string>
+ </property>
+ </widget>
+ </widget>
+ </item>
+ <item>
+ <widget class="QLineEdit" name="tbInfoOutputImages">
+ <property name="autoFillBackground">
+ <bool>false</bool>
+ </property>
+ <property name="styleSheet">
+ <string notr="true">background-color:transparent;</string>
+ </property>
+ <property name="text">
+ <string/>
+ </property>
+ <property name="frame">
+ <bool>false</bool>
+ </property>
+ <property name="readOnly">
+ <bool>true</bool>
+ </property>
+ <property name="placeholderText">
+ <string>No output images defined</string>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </widget>
+ </item>
+ <item>
+ <spacer name="verticalSpacer">
+ <property name="orientation">
+ <enum>Qt::Vertical</enum>
+ </property>
+ <property name="sizeHint" stdset="0">
+ <size>
+ <width>0</width>
+ <height>0</height>
+ </size>
+ </property>
+ </spacer>
+ </item>
+ <item>
+ <layout class="QHBoxLayout" name="horizontalLayout_3">
+ <property name="topMargin">
+ <number>0</number>
+ </property>
+ <item>
+ <widget class="QProgressBar" name="progressBar">
+ <property name="value">
+ <number>0</number>
+ </property>
+ </widget>
+ </item>
+ <item>
+ <widget class="QDialogButtonBox" name="buttonBox">
+ <property name="sizePolicy">
+ <sizepolicy hsizetype="Expanding" vsizetype="Preferred">
+ <horstretch>0</horstretch>
+ <verstretch>0</verstretch>
+ </sizepolicy>
+ </property>
+ <property name="orientation">
+ <enum>Qt::Horizontal</enum>
+ </property>
+ <property name="standardButtons">
+ <set>QDialogButtonBox::Cancel|QDialogButtonBox::Save</set>
+ </property>
+ <property name="centerButtons">
+ <bool>false</bool>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </item>
+ </layout>
+ <action name="actionAddSourceStack">
+ <property name="icon">
+ <iconset resource="resources.qrc">
+ <normaloff>:/timeseriesviewer/icons/mActionAddRasterLayer.svg</normaloff>:/timeseriesviewer/icons/mActionAddRasterLayer.svg</iconset>
+ </property>
+ <property name="text">
+ <string>AddSourceStack</string>
+ </property>
+ </action>
+ <action name="actionRemoveSourceStack">
+ <property name="icon">
+ <iconset resource="resources.qrc">
+ <normaloff>:/timeseriesviewer/icons/mActionRemoveTSD.svg</normaloff>:/timeseriesviewer/icons/mActionRemoveTSD.svg</iconset>
+ </property>
+ <property name="text">
+ <string>removeSourceStack</string>
+ </property>
+ </action>
+ </widget>
+ <customwidgets>
+ <customwidget>
+ <class>QgsCollapsibleGroupBox</class>
+ <extends>QGroupBox</extends>
+ <header>qgscollapsiblegroupbox.h</header>
+ <container>1</container>
+ </customwidget>
+ <customwidget>
+ <class>QgsFileWidget</class>
+ <extends>QWidget</extends>
+ <header>qgsfilewidget.h</header>
+ </customwidget>
+ </customwidgets>
+ <resources>
+ <include location="resources.qrc"/>
+ </resources>
+ <connections>
+ <connection>
+ <sender>buttonBox</sender>
+ <signal>accepted()</signal>
+ <receiver>Dialog</receiver>
+ <slot>accept()</slot>
+ <hints>
+ <hint type="sourcelabel">
+ <x>236</x>
+ <y>587</y>
+ </hint>
+ <hint type="destinationlabel">
+ <x>157</x>
+ <y>274</y>
+ </hint>
+ </hints>
+ </connection>
+ <connection>
+ <sender>buttonBox</sender>
+ <signal>rejected()</signal>
+ <receiver>Dialog</receiver>
+ <slot>reject()</slot>
+ <hints>
+ <hint type="sourcelabel">
+ <x>304</x>
+ <y>587</y>
+ </hint>
+ <hint type="destinationlabel">
+ <x>286</x>
+ <y>274</y>
+ </hint>
+ </hints>
+ </connection>
+ <connection>
+ <sender>rbSaveInDirectory</sender>
+ <signal>toggled(bool)</signal>
+ <receiver>fileWidgetOutputDir</receiver>
+ <slot>setEnabled(bool)</slot>
+ <hints>
+ <hint type="sourcelabel">
+ <x>217</x>
+ <y>238</y>
+ </hint>
+ <hint type="destinationlabel">
+ <x>375</x>
+ <y>240</y>
+ </hint>
+ </hints>
+ </connection>
+ </connections>
+ <buttongroups>
+ <buttongroup name="buttonGroupOutputLocation"/>
+ <buttongroup name="buttonGroupDateMode"/>
+ </buttongroups>
+</ui>
diff --git a/timeseriesviewer/ui/timeseriesviewer.ui b/timeseriesviewer/ui/timeseriesviewer.ui
index 337848a8954bfceba9df68d9aa3e9ec568abe95e..7e5d7708e10b62de42faa0574407e8caf96fbf24 100644
--- a/timeseriesviewer/ui/timeseriesviewer.ui
+++ b/timeseriesviewer/ui/timeseriesviewer.ui
@@ -155,6 +155,7 @@
<string>Files</string>
</property>
<addaction name="actionAddTSD"/>
+ <addaction name="actionLoadTimeSeriesStack"/>
<addaction name="actionAddVectorData"/>
<addaction name="actionLoadTS"/>
<addaction name="actionSaveTS"/>
@@ -688,6 +689,14 @@
<string>Opens the online documentation</string>
</property>
</action>
+ <action name="actionLoadTimeSeriesStack">
+ <property name="text">
+ <string>Add images from time stack</string>
+ </property>
+ <property name="toolTip">
+ <string>Load images from a stack or list of stacks witch each image band being an temporal observation</string>
+ </property>
+ </action>
</widget>
<customwidgets>
<customwidget>
diff --git a/timeseriesviewer/utils.py b/timeseriesviewer/utils.py
index f30813db35ad40caf6f30ac357125371c23a384d..b9684e6d3de260afcc62b2f0344c3efdaa0d1603 100644
--- a/timeseriesviewer/utils.py
+++ b/timeseriesviewer/utils.py
@@ -101,6 +101,35 @@ def file_search(rootdir, pattern, recursive=False, ignoreCase=False, directories
return results
+NEXT_COLOR_HUE_DELTA_CON = 10
+NEXT_COLOR_HUE_DELTA_CAT = 100
+def nextColor(color, mode='cat'):
+ """
+ Returns another color
+ :param color: the previous color
+ :param mode: 'cat' - for categorical color jump (next color looks pretty different to previous)
+ 'con' - for continuous color jump (next color looks similar to previous)
+ :return:
+ """
+ assert mode in ['cat','con']
+ assert isinstance(color, QColor)
+ hue, sat, value, alpha = color.getHsl()
+ if mode == 'cat':
+ hue += NEXT_COLOR_HUE_DELTA_CAT
+ elif mode == 'con':
+ hue += NEXT_COLOR_HUE_DELTA_CON
+ if sat == 0:
+ sat = 255
+ value = 128
+ alpha = 255
+ s = ""
+ while hue > 360:
+ hue -= 360
+
+ return QColor.fromHsl(hue, sat, value, alpha)
+
+
+
def createQgsField(name : str, exampleValue, comment:str=None):
@@ -118,12 +147,12 @@ def createQgsField(name : str, exampleValue, comment:str=None):
return QgsField(name, QVariant.Bool, 'int', len=1, comment=comment)
elif t in [int, np.int32, np.int64]:
return QgsField(name, QVariant.Int, 'int', comment=comment)
- elif t in [float, np.double, np.float, np.float64]:
+ elif t in [float, np.double, np.float16, np.float32, np.float64]:
return QgsField(name, QVariant.Double, 'double', comment=comment)
elif isinstance(exampleValue, np.ndarray):
return QgsField(name, QVariant.String, 'varchar', comment=comment)
elif isinstance(exampleValue, list):
- assert len(exampleValue)> 0, 'need at least one value in provided list'
+ assert len(exampleValue) > 0, 'need at least one value in provided list'
v = exampleValue[0]
prototype = createQgsField(name, v)
subType = prototype.type()
@@ -133,6 +162,7 @@ def createQgsField(name : str, exampleValue, comment:str=None):
raise NotImplemented()
+
def setQgsFieldValue(feature:QgsFeature, field, value):
"""
Wrties the Python value v into a QgsFeature field, taking care of required conversions
@@ -437,6 +467,9 @@ def px2geo(px, gt, pxCenter=True):
return QgsPointXY(gx, gy)
+
+
+
class SpatialExtent(QgsRectangle):
"""
Object to keep QgsRectangle and QgsCoordinateReferenceSystem together
diff --git a/timeseriesviewer/virtualrasters.py b/timeseriesviewer/virtualrasters.py
index 5de55b0a7102dcab90a000a9638f799748dd1e9f..8006d801753c566a257a08e68f24b59b7f3bfc4c 100644
--- a/timeseriesviewer/virtualrasters.py
+++ b/timeseriesviewer/virtualrasters.py
@@ -2,8 +2,8 @@
# noinspection PyPep8Naming
"""
/***************************************************************************
- EO Time Series Viewer
- -------------------
+ Virtual Raster Builder
+ ----------------------
begin : 2015-08-20
git sha : $Format:%H$
copyright : (C) 2017 by HU-Berlin
@@ -19,7 +19,8 @@
* *
***************************************************************************/
"""
-import os, sys, re, pickle, tempfile, unicodedata
+import os, sys, re, pickle, tempfile, uuid
+from xml.etree import ElementTree
from collections import OrderedDict
import tempfile
from osgeo import gdal, osr, ogr, gdalconst as gc
@@ -56,26 +57,56 @@ LUT_GDT_NAME = {gdal.GDT_Byte:'Byte',
gdal.GDT_CFloat64:'Float64'}
+GRA_tooltips = {'NearestNeighbour':'nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).',
+ 'Bilinear':'bilinear resampling.',
+ 'Lanczos':'lanczos windowed sinc resampling.',
+ 'Average':'average resampling, computes the average of all non-NODATA contributing pixels.',
+ 'Cubic':'cubic resampling.',
+ 'CubicSpline':'cubic spline resampling.',
+ 'Mode':'mode resampling, selects the value which appears most often of all the sampled points',
+ 'Max':'maximum resampling, selects the maximum value from all non-NODATA contributing pixels',
+ 'Min':'minimum resampling, selects the minimum value from all non-NODATA contributing pixels.',
+ 'Med':'median resampling, selects the median value of all non-NODATA contributing pixels.',
+ 'Q1':'first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels. ',
+ 'Q3':'third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels'
+ }
RESAMPLE_ALGS = OptionListModel()
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_NearestNeighbour, 'nearest',
- tooltip='nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_Bilinear, 'bilinear', tooltip='bilinear resampling.'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_Lanczos, 'lanczos', tooltip='lanczos windowed sinc resampling.'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_Average, 'average', tooltip='average resampling, computes the average of all non-NODATA contributing pixels.'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_Cubic, 'cubic', tooltip='cubic resampling.'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_CubicSpline, 'cubic_spline', tooltip='cubic spline resampling.'))
-RESAMPLE_ALGS.addOption(Option(gdal.GRA_Mode, 'mode', tooltip='mode resampling, selects the value which appears most often of all the sampled points'))
-
-if int(gdal.VersionInfo()) >= 2020300:
- #respect that older GDAL version do not have python binding to GRA_Max, GRA_Min, GRA_Med, GRA_Q1 and GRA_Q3
- #see https://trac.osgeo.org/gdal/wiki/Release/2.2.3-News
- #https://trac.osgeo.org/gdal/ticket/7153
- RESAMPLE_ALGS.addOption(Option(gdal.GRA_Max, 'max', tooltip='maximum resampling, selects the maximum value from all non-NODATA contributing pixels'))
- RESAMPLE_ALGS.addOption(Option(gdal.GRA_Min, 'min', tooltip='minimum resampling, selects the minimum value from all non-NODATA contributing pixels.'))
- RESAMPLE_ALGS.addOption(Option(gdal.GRA_Med, 'med', tooltip='median resampling, selects the median value of all non-NODATA contributing pixels.'))
- RESAMPLE_ALGS.addOption(Option(gdal.GRA_Q1, 'Q1', tooltip='first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels. '))
- RESAMPLE_ALGS.addOption(Option(gdal.GRA_Q3, 'Q2', tooltip='third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels'))
+for GRAkey in [k for k in list(gdal.__dict__.keys()) if k.startswith('GRA_')]:
+
+ GRA = gdal.__dict__[GRAkey]
+ GRA_Name = GRAkey[4:]
+
+ option = Option(GRA, GRA_Name, tooltip=GRA_tooltips.get(GRA_Name))
+ RESAMPLE_ALGS.addOption(option)
+
+
+# thanks to https://gis.stackexchange.com/questions/75533/how-to-apply-band-settings-using-gdal-python-bindings
+def read_vsimem(fn):
+ """
+ Reads VSIMEM path as string
+ :param fn: vsimem path (str)
+ :return: result of gdal.VSIFReadL(1, vsileng, vsifile)
+ """
+ vsifile = gdal.VSIFOpenL(fn,'r')
+ gdal.VSIFSeekL(vsifile, 0, 2)
+ vsileng = gdal.VSIFTellL(vsifile)
+ gdal.VSIFSeekL(vsifile, 0, 0)
+ return gdal.VSIFReadL(1, vsileng, vsifile)
+
+def write_vsimem(fn:str,data:str):
+ """
+ Writes data to vsimem path
+ :param fn: vsimem path (str)
+ :param data: string to write
+ :return: result of gdal.VSIFCloseL(vsifile)
+ """
+ '''Write GDAL vsimem files'''
+ vsifile = gdal.VSIFOpenL(fn,'w')
+ size = len(data)
+ gdal.VSIFWriteL(data, 1, size, vsifile)
+ return gdal.VSIFCloseL(vsifile)
+
def px2geo(px, gt):
#see http://www.gdal.org/gdal_datamodel.html
@@ -113,7 +144,13 @@ def describeRawFile(pathRaw, pathVrt, xsize, ysize,
assert interleave in ['bsq','bil','bip']
assert byteOrder in ['LSB', 'MSB']
- vrt = ['<VRTDataset rasterXSize="{xsize}" rasterYSize="{ysize}">'.format(xsize=xsize,ysize=ysize)]
+
+ drvVRT = gdal.GetDriverByName('VRT')
+ assert isinstance(drvVRT, gdal.Driver)
+ dsVRT = drvVRT.Create(pathVrt, xsize, ysize, bands=0, eType=eType)
+ assert isinstance(dsVRT, gdal.Dataset)
+
+ #vrt = ['<VRTDataset rasterXSize="{xsize}" rasterYSize="{ysize}">'.format(xsize=xsize,ysize=ysize)]
vrtDir = os.path.dirname(pathVrt)
if pathRaw.startswith(vrtDir):
@@ -124,8 +161,6 @@ def describeRawFile(pathRaw, pathVrt, xsize, ysize,
srcFilename = pathRaw
for b in range(bands):
- vrt.append(' <VRTRasterBand dataType="{dataType}" band="{band}" subClass="VRTRawRasterBand">'.format(
- dataType=LUT_GDT_NAME[eType], band=b+1))
if interleave == 'bsq':
imageOffset = headerOffset
pixelOffset = LUT_GDT_SIZE[eType]
@@ -136,67 +171,47 @@ def describeRawFile(pathRaw, pathVrt, xsize, ysize,
lineOffset = xsize * bands
else:
raise Exception('Interleave {} is not supported'.format(interleave))
- vrt.append(""" <SourceFilename relativetoVRT="{relativeToVRT}">{srcFilename}</SourceFilename>
- <ImageOffset>{imageOffset}</ImageOffset>
- <PixelOffset>{pixelOffset}</PixelOffset>
- <LineOffset>{lineOffset}</LineOffset>
- <ByteOrder>{byteOrder}</ByteOrder>""".format(relativeToVRT=relativeToVRT,
- srcFilename=srcFilename,
- imageOffset=imageOffset,
- pixelOffset=pixelOffset,
- lineOffset=lineOffset,
- byteOrder=byteOrder))
-
- vrt.append(' </VRTRasterBand>')
- vrt.append('</VRTDataset>')
- vrt = '\n'.join(vrt)
-
- drv = gdal.GetDriverByName('VRT')
- assert isinstance(drv, gdal.Driver)
-
- dsVRT = drv.Create(pathVrt, xsize, ysize, bands , eType)
- assert isinstance(dsVRT, gdal.Dataset)
- for b in range(bands):
- band = dsVRT.GetRasterBand(b+1)
- assert isinstance(band, gdal.Band)
- if interleave == 'bsq':
- imageOffset = headerOffset
- pixelOffset = LUT_GDT_SIZE[eType]
- lineOffset = pixelOffset * xsize
- elif interleave == 'bip':
- imageOffset = headerOffset + b * LUT_GDT_SIZE[eType]
- pixelOffset = bands * LUT_GDT_SIZE[eType]
- lineOffset = xsize * bands
- else:
- raise Exception('Interleave {} is not supported'.format(interleave))
+ options = ['subClass=VRTRawRasterBand']
+ options.append('SourceFilename={}'.format(srcFilename))
+ options.append('dataType={}'.format(LUT_GDT_NAME[eType]))
+ options.append('ImageOffset={}'.format(imageOffset))
+ options.append('PixelOffset={}'.format(pixelOffset))
+ options.append('LineOffset={}'.format(lineOffset))
+ options.append('ByteOrder={}'.format(byteOrder))
xml = """<SourceFilename relativetoVRT="{relativeToVRT}">{srcFilename}</SourceFilename>
<ImageOffset>{imageOffset}</ImageOffset>
<PixelOffset>{pixelOffset}</PixelOffset>
<LineOffset>{lineOffset}</LineOffset>
<ByteOrder>{byteOrder}</ByteOrder>""".format(relativeToVRT=relativeToVRT,
- srcFilename=srcFilename,
- imageOffset=imageOffset,
- pixelOffset=pixelOffset,
- lineOffset=lineOffset,
- byteOrder=byteOrder)
-
- band.SetMetadataItem('source_0', xml, "new_vrt_sources")
+ srcFilename=srcFilename,
+ imageOffset=imageOffset,
+ pixelOffset=pixelOffset,
+ lineOffset=lineOffset,
+ byteOrder=byteOrder)
+
+ #md = {}
+ #md['source_0'] = xml
+ #vrtBand = dsVRT.GetRasterBand(b + 1)
+ assert dsVRT.AddBand(eType, options=options) == 0
+
+ vrtBand = dsVRT.GetRasterBand(b+1)
+ assert isinstance(vrtBand, gdal.Band)
+ #vrtBand.SetMetadata(md, 'vrt_sources')
+ #vrt.append(' <VRTRasterBand dataType="{dataType}" band="{band}" subClass="VRTRawRasterBand">'.format(dataType=LUT_GDT_NAME[eType], band=b+1))
dsVRT.FlushCache()
-
- #file = open(pathVrt, 'w')
- #file.write(vrt)
- #file.flush()
- #file.close()
- #ds = gdal.Open(pathVrt)
- #return ds
return dsVRT
class VRTRasterInputSourceBand(object):
@staticmethod
def fromGDALDataSet(pathOrDataSet):
+ """
+ Returns the VRTRasterInputSourceBands from a raster data source
+ :param pathOrDataSet: str | gdal.Dataset
+ :return: [list-of-VRTRasterInputSourceBand]
+ """
srcBands = []
@@ -211,15 +226,14 @@ class VRTRasterInputSourceBand(object):
- def __init__(self, path, bandIndex, bandName=''):
- self.mPath = os.path.normpath(path)
+ def __init__(self, path:str, bandIndex:int, bandName:str=''):
+ self.mPath = path
self.mBandIndex = bandIndex
self.mBandName = bandName
self.mNoData = None
self.mVirtualBand = None
-
def isEqual(self, other):
if isinstance(other, VRTRasterInputSourceBand):
return self.mPath == other.mPath and self.mBandIndex == other.mBandIndex
@@ -248,16 +262,16 @@ class VRTRasterBand(QObject):
sigSourceRemoved = pyqtSignal(int, VRTRasterInputSourceBand)
def __init__(self, name='', parent=None):
super(VRTRasterBand, self).__init__(parent)
- self.sources = []
+ self.mSources = []
self.mName = ''
self.setName(name)
self.mVRT = None
+ def __len__(self):
+ return len(self.mSources)
def setName(self, name):
-
- if isinstance(name, str):
- name = unicode(name)
+ assert isinstance(name, str)
oldName = self.mName
self.mName = name
if oldName != self.mName:
@@ -270,16 +284,17 @@ class VRTRasterBand(QObject):
def addSource(self, virtualBandInputSource):
assert isinstance(virtualBandInputSource, VRTRasterInputSourceBand)
- self.insertSource(len(self.sources), virtualBandInputSource)
+ self.insertSource(len(self.mSources), virtualBandInputSource)
def insertSource(self, index, virtualBandInputSource):
assert isinstance(virtualBandInputSource, VRTRasterInputSourceBand)
virtualBandInputSource.mVirtualBand = self
- if index <= len(self.sources):
- self.sources.insert(index, virtualBandInputSource)
+ if index <= len(self.mSources):
+ self.mSources.insert(index, virtualBandInputSource)
self.sigSourceInserted.emit(index, virtualBandInputSource)
else:
- print('DEBUG: index <= len(self.sources)')
+ pass
+ #print('DEBUG: index <= len(self.sources)')
def bandIndex(self):
if isinstance(self.mVRT, VRTRaster):
return self.mVRT.mBands.index(self)
@@ -294,10 +309,10 @@ class VRTRasterBand(QObject):
:return: The VRTRasterInputSourceBand that was removed
"""
if not isinstance(vrtRasterInputSourceBand, VRTRasterInputSourceBand):
- vrtRasterInputSourceBand = self.sources[vrtRasterInputSourceBand]
- if vrtRasterInputSourceBand in self.sources:
- i = self.sources.index(vrtRasterInputSourceBand)
- self.sources.remove(vrtRasterInputSourceBand)
+ vrtRasterInputSourceBand = self.mSources[vrtRasterInputSourceBand]
+ if vrtRasterInputSourceBand in self.mSources:
+ i = self.mSources.index(vrtRasterInputSourceBand)
+ self.mSources.remove(vrtRasterInputSourceBand)
self.sigSourceRemoved.emit(i, vrtRasterInputSourceBand)
@@ -305,12 +320,12 @@ class VRTRasterBand(QObject):
"""
:return: list of file-paths to all source files
"""
- files = set([inputSource.mPath for inputSource in self.sources])
+ files = set([inputSource.mPath for inputSource in self.mSources])
return sorted(list(files))
def __repr__(self):
infos = ['VirtualBand name="{}"'.format(self.mName)]
- for i, info in enumerate(self.sources):
+ for i, info in enumerate(self.mSources):
assert isinstance(info, VRTRasterInputSourceBand)
infos.append('\t{} SourceFileName {} SourceBand {}'.format(i + 1, info.mPath, info.mBandIndex))
return '\n'.join(infos)
@@ -449,7 +464,8 @@ class VRTRaster(QObject):
if crs != self.mCrs:
extent = self.extent()
if isinstance(extent, QgsRectangle):
- trans = QgsCoordinateTransform(self.mCrs, crs)
+ trans = QgsCoordinateTransform()
+ trans.setDestinationCrs(self.mCrs, crs)
extent = trans.transform(extent)
self.setExtent(extent)
self.mCrs = crs
@@ -527,7 +543,7 @@ class VRTRaster(QObject):
assert path in self.sourceRaster()
for vBand in self.mBands:
assert isinstance(vBand, VRTRasterBand)
- if path in vBand.sources():
+ if path in vBand.mSources():
vBand.removeSource(path)
def removeVirtualBand(self, bandOrIndex):
@@ -621,14 +637,14 @@ class VRTRaster(QObject):
ds = gdal.Open(pathVRT)
assert isinstance(ds, gdal.Dataset)
assert ds.GetDriver().GetDescription() == 'VRT'
- from xml.etree import ElementTree
+
for b in range(ds.RasterCount):
srcBand = ds.GetRasterBand(b+1)
vrtBand = VRTRasterBand(name=srcBand.GetDescription().decode('utf-8'))
for key, xml in srcBand.GetMetadata(str('vrt_sources')).items():
tree = ElementTree.fromstring(xml)
- srcPath = os.path.normpath(tree.find('SourceFilename').text)
+ srcPath = tree.find('SourceFilename').text
srcBandIndex = int(tree.find('SourceBand').text)
vrtBand.addSource(VRTRasterInputSourceBand(srcPath, srcBandIndex))
@@ -638,22 +654,32 @@ class VRTRaster(QObject):
- def saveVRT(self, pathVRT):
-
- if len(self.mBands) == 0:
- print('No VRT Inputs defined.')
- return None
+ def saveVRT(self, pathVRT, warpedImageFolder = '.warpedimage'):
+ """
+ Save the VRT to path.
+ If source images need to be warped to the final CRS warped VRT image will be created in a folder <directory>/<basename>+<warpedImageFolder>/
+ :param pathVRT: str, path of final VRT.
+ :param warpedImageFolder: basename of folder that is created
+ :return:
+ """
+ """
+ :param pathVRT:
+ :return:
+ """
+ assert len(self) >= 1, 'VRT needs to define at least 1 band'
assert os.path.splitext(pathVRT)[-1].lower() == '.vrt'
- dirVrt = os.path.dirname(pathVRT)
- dirWarpedVRT = os.path.join(dirVrt, 'WarpedVRTs')
- if not os.path.isdir(dirVrt):
- os.mkdir(dirVrt)
-
srcLookup = dict()
srcNodata = None
+ inMemory = pathVRT.startswith('/vsimem/')
+
+ if inMemory:
+ dirWarped = '/vsimem/'
+ else:
+ dirWarped = os.path.join(os.path.splitext(pathVRT)[0] + '.WarpedImages')
+ drvVRT = gdal.GetDriverByName('VRT')
for i, pathSrc in enumerate(self.sourceRaster()):
dsSrc = gdal.Open(pathSrc)
assert isinstance(dsSrc, gdal.Dataset)
@@ -667,68 +693,114 @@ class VRTRaster(QObject):
if crs == self.mCrs:
srcLookup[pathSrc] = pathSrc
else:
+ #do a CRS transformation using VRTs
+
+ warpedFileName = 'warped.{}.vrt'.format(os.path.basename(pathSrc))
+ if inMemory:
+ warpedFileName = dirWarped + warpedFileName
+ else:
+ os.makedirs(dirWarped, exist_ok=True)
+ warpedFileName = os.path.join(dirWarped, warpedFileName)
- if not os.path.isdir(dirWarpedVRT):
- os.mkdir(dirWarpedVRT)
- pathVRT2 = os.path.join(dirWarpedVRT, 'warped.{}.vrt'.format(os.path.basename(pathSrc)))
wops = gdal.WarpOptions(format='VRT',
dstSRS=self.mCrs.toWkt())
- tmp = gdal.Warp(pathVRT2, dsSrc, options=wops)
+ tmp = gdal.Warp(warpedFileName, dsSrc, options=wops)
assert isinstance(tmp, gdal.Dataset)
- tmp = None
- srcLookup[pathSrc] = pathVRT2
+ vrtXML = read_vsimem(warpedFileName)
+ xml = ElementTree.fromstring(vrtXML)
+ #print(vrtXML.decode('utf-8'))
- srcFiles = [srcLookup[src] for src in self.sourceRaster()]
+ if False:
+ dsTmp = gdal.Open(warpedFileName)
+ assert isinstance(dsTmp, gdal.Dataset)
+ drvVRT.Delete(warpedFileName)
+ dsTmp = gdal.Open(warpedFileName)
+ assert not isinstance(dsTmp, gdal.Dataset)
- #1. build a temporary VRT that describes the spatial shifts of all input sources
+ srcLookup[pathSrc] = warpedFileName
+ srcFiles = [srcLookup[src] for src in self.sourceRaster()]
- kwds = {}
+ #these need to be set
+ ns = nl = gt = crs = eType = None
res = self.resolution()
+ extent = self.extent()
- if res is None:
- res = 'average'
- if isinstance(res, QSizeF):
- kwds['resolution'] = 'user'
- kwds['xRes'] = res.width()
- kwds['yRes'] = res.height()
- else:
- assert res in ['highest','lowest','average']
- kwds['resolution'] = res
+ srs = None
+ if isinstance(self.crs(), QgsCoordinateReferenceSystem):
+ srs = self.crs().toWkt()
+
+ if len(srcFiles) > 0:
+ # 1. build a temporary VRT that describes the spatial shifts of all input sources
+ kwds = {}
+ if res is None:
+ res = 'average'
+ if isinstance(res, QSizeF):
+ kwds['resolution'] = 'user'
+ kwds['xRes'] = res.width()
+ kwds['yRes'] = res.height()
+ else:
+ assert res in ['highest','lowest','average']
+ kwds['resolution'] = res
- extent = self.extent()
- if isinstance(extent, QgsRectangle):
- kwds['outputBounds'] = (extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum())
+ if isinstance(extent, QgsRectangle):
+ kwds['outputBounds'] = (extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum())
- vro = gdal.BuildVRTOptions(separate=True, **kwds)
+ if srs is not None:
+ kwds['outputSRS'] = srs
- gdal.BuildVRT(pathVRT, srcFiles, options=vro)
- dsVRTDst = gdal.Open(pathVRT)
- assert isinstance(dsVRTDst, gdal.Dataset)
- assert len(srcLookup) == dsVRTDst.RasterCount
- ns, nl = dsVRTDst.RasterXSize, dsVRTDst.RasterYSize
- gt = dsVRTDst.GetGeoTransform()
- crs = dsVRTDst.GetProjectionRef()
- eType = dsVRTDst.GetRasterBand(1).DataType
- SOURCE_TEMPLATES = dict()
- for i, srcFile in enumerate(srcFiles):
- vrt_sources = dsVRTDst.GetRasterBand(i+1).GetMetadata(str('vrt_sources'))
- assert len(vrt_sources) == 1
- srcXML = vrt_sources['source_0']
- assert os.path.basename(srcFile)+'</SourceFilename>' in srcXML
- assert '<SourceBand>1</SourceBand>' in srcXML
- SOURCE_TEMPLATES[srcFile] = srcXML
- dsVRTDst = None
- #remove the temporary VRT, we don't need it any more
- os.remove(pathVRT)
+
+
+ pathInMEMVRT = '/vsimem/{}.vrt'.format(uuid.uuid4())
+ vro = gdal.BuildVRTOptions(separate=True, **kwds)
+ dsVRTDst = gdal.BuildVRT(pathInMEMVRT, srcFiles, options=vro)
+
+ assert isinstance(dsVRTDst, gdal.Dataset)
+
+ ns, nl = dsVRTDst.RasterXSize, dsVRTDst.RasterYSize
+ gt = dsVRTDst.GetGeoTransform()
+ crs = dsVRTDst.GetProjectionRef()
+ eType = dsVRTDst.GetRasterBand(1).DataType
+ SOURCE_TEMPLATES = dict()
+ for i, srcFile in enumerate(srcFiles):
+ vrt_sources = dsVRTDst.GetRasterBand(i+1).GetMetadata(str('vrt_sources'))
+ assert len(vrt_sources) == 1
+ srcXML = vrt_sources['source_0']
+ assert os.path.basename(srcFile)+'</SourceFilename>' in srcXML
+ assert '<SourceBand>1</SourceBand>' in srcXML
+ SOURCE_TEMPLATES[srcFile] = srcXML
+
+ drvVRT.Delete(pathInMEMVRT)
+
+ else:
+ # special case: no source files defined
+ ns = nl = 1 #this is the minimum size
+ if isinstance(extent, QgsRectangle):
+ x0 = extent.xMinimum()
+ y1 = extent.yMaximum()
+ else:
+ x0 = 0
+ y1 = 0
+
+ if isinstance(res, QSizeF):
+ resx = res.width()
+ resy = res.height()
+ else:
+ resx = 1
+ resy = 1
+
+ gt = (x0, resx, 0, y1, 0, -resy)
+ eType = gdal.GDT_Float32
#2. build final VRT from scratch
- drvVRT = gdal.GetDriverByName(str('VRT'))
+ drvVRT = gdal.GetDriverByName('VRT')
assert isinstance(drvVRT, gdal.Driver)
dsVRTDst = drvVRT.Create(pathVRT, ns, nl,0, eType=eType)
#2.1. set general properties
assert isinstance(dsVRTDst, gdal.Dataset)
- dsVRTDst.SetProjection(crs)
+
+ if srs is not None:
+ dsVRTDst.SetProjection(srs)
dsVRTDst.SetGeoTransform(gt)
#2.2. add virtual bands
@@ -740,23 +812,21 @@ class VRTRaster(QObject):
vrtBandDst.SetDescription(vBand.name())
md = {}
#add all input sources for this virtual band
- for iSrc, sourceInfo in enumerate(vBand.sources):
+ for iSrc, sourceInfo in enumerate(vBand.mSources):
assert isinstance(sourceInfo, VRTRasterInputSourceBand)
bandIndex = sourceInfo.mBandIndex
xml = SOURCE_TEMPLATES[srcLookup[sourceInfo.mPath]]
- xml = re.sub('<SourceBand>1</SourceBand>','<SourceBand>{}</SourceBand>'.format(bandIndex+1), xml)
+ xml = re.sub('<SourceBand>1</SourceBand>', '<SourceBand>{}</SourceBand>'.format(bandIndex+1), xml)
md['source_{}'.format(iSrc)] = xml
- vrtBandDst.SetMetadata(md,str('vrt_sources'))
- if False:
- vrtBandDst.ComputeBandStats(1)
+ vrtBandDst.SetMetadata(md,'vrt_sources')
dsVRTDst = None
#check if we get what we like to get
dsCheck = gdal.Open(pathVRT)
+ assert isinstance(dsCheck, gdal.Dataset)
- s = ""
return dsCheck
def __repr__(self):