diff --git a/CHANGELOG b/CHANGELOG index fb1f7beb97cf5fb4445f889a860f9b9f1acd9523..dcb5b8d268b4e671d74aca8de83a964019a8b8a1 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,6 +1,10 @@ ============== Changelog ============== +2019-08-06 (version 1.7): + * increased contrast for default map view text + * improved detect of wavelength information, e.g. from Pleiades, Sentinel-2 and RapidEye data + 2019-07-16 (version 1.6): * re-design of map visualization: faster and more compact, the number of maps is fixed to n dates x m map views * date, sensor or map view information can be plotted within each map and become available in screenshots diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml index a36fc9627673d99b9d7214df3462e685725e82d4..abb6ec1128a99211357dad2070e4cf3d90fee182 100644 --- a/bitbucket-pipelines.yml +++ b/bitbucket-pipelines.yml @@ -17,7 +17,7 @@ pipelines: - python -m pip install -r requirements.txt - python -m pip install nose2 - apt-get update - - apt-get install xvfb + - apt-get -y install xvfb - Xvfb :1 -screen 0 1024x768x16 &> xvfb.log & - ps aux | grep X - DISPLAY=:1.0 diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index ce802c2b93a407f0278d373a5ecbc9bb51a88bfc..36f46eee72cef3cd8a8d3b191b6b1d11903db917 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -1,6 +1,10 @@ ============== Changelog ============== +2019-08-06 (version 1.7): + * increased contrast for default map view text + * improved detect of wavelength information, e.g. from Pleiades, Sentinel-2 and RapidEye data + 2019-07-16 (version 1.6): * re-design of map visualization: faster and more compact, the number of maps is fixed to n dates x m map views * date, sensor or map view information can be plotted within each map and become available in screenshots diff --git a/eotimeseriesviewer/__init__.py b/eotimeseriesviewer/__init__.py index 3f15799e8efc3750d1b4f8142f4fe28d61951569..36175536466d547efbf9070f4cb550c3acd31d41 100644 --- a/eotimeseriesviewer/__init__.py +++ b/eotimeseriesviewer/__init__.py @@ -21,7 +21,7 @@ # noinspection PyPep8Naming -__version__ = '1.6' # sub-subversion number is added automatically +__version__ = '1.7' # sub-subversion number is added automatically LICENSE = 'GNU GPL-3' TITLE = 'EO Time Series Viewer' DESCRIPTION = 'Visualization of multi-sensor Earth observation time series data.' diff --git a/eotimeseriesviewer/main.py b/eotimeseriesviewer/main.py index 7ae8cd3080d9f121e8b4730e0dd4129fb1e8f45d..9dcde82d130cb6591348af70d933dd5526683026 100644 --- a/eotimeseriesviewer/main.py +++ b/eotimeseriesviewer/main.py @@ -754,10 +754,10 @@ class TimeSeriesViewer(QgisInterface, QObject): iface = qgis.utils.iface assert isinstance(iface, QgisInterface) - self.ui.actionImportExtent.triggered.connect(lambda: self.spatialTemporalVis.setSpatialExtent(SpatialExtent.fromMapCanvas(iface.mapCanvas()))) - self.ui.actionExportExtent.triggered.connect(lambda: iface.mapCanvas().setExtent(self.spatialTemporalVis.spatialExtent().toCrs(iface.mapCanvas().mapSettings().destinationCrs()))) - self.ui.actionExportCenter.triggered.connect(lambda: iface.mapCanvas().setCenter(self.spatialTemporalVis.spatialExtent().spatialCenter())) - self.ui.actionImportCenter.triggered.connect(lambda: self.spatialTemporalVis.setSpatialCenter(SpatialPoint.fromMapCanvasCenter(iface.mapCanvas()))) + self.ui.actionImportExtent.triggered.connect(lambda: self.setSpatialExtent(SpatialExtent.fromMapCanvas(iface.mapCanvas()))) + self.ui.actionExportExtent.triggered.connect(lambda: iface.mapCanvas().setExtent(self.spatialExtent().toCrs(iface.mapCanvas().mapSettings().destinationCrs()))) + self.ui.actionExportCenter.triggered.connect(lambda: iface.mapCanvas().setCenter(self.spatialCenter().toCrs(iface.mapCanvas().mapSettings().destinationCrs()))) + self.ui.actionImportCenter.triggered.connect(lambda: self.setSpatialCenter(SpatialPoint.fromMapCanvasCenter(iface.mapCanvas()))) def onSyncRequest(qgisChanged:bool): if self.ui.optionSyncMapCenter.isChecked(): @@ -865,7 +865,18 @@ class TimeSeriesViewer(QgisInterface, QObject): """ self.mapWidget().setSpatialCenter(spatialPoint) + def spatialExtent(self)->SpatialExtent: + """ + Returns the map extent + :return: SpatialExtent + """ + return self.mapWidget().spatialExtent() + def spatialCenter(self)->SpatialPoint: + """ + Returns the map center + :return: SpatialPoint + """ return self.mapWidget().spatialCenter() def setCurrentLocation(self, spatialPoint:SpatialPoint, mapCanvas:QgsMapCanvas=None): diff --git a/eotimeseriesviewer/mapcanvas.py b/eotimeseriesviewer/mapcanvas.py index e6fe7b8c87412e5412b0906c0058489db721b47a..bc856e34775886d9f64f1a56fa1c0ae424e921be 100644 --- a/eotimeseriesviewer/mapcanvas.py +++ b/eotimeseriesviewer/mapcanvas.py @@ -998,17 +998,29 @@ class MapCanvas(QgsMapCanvas): m = menu.addMenu('Copy...') action = m.addAction('Date') action.triggered.connect(lambda: QApplication.clipboard().setText(str(tsd.date()))) - action.setToolTip('Sends "" to the clipboard.'.format(str(tsd.date()))) + action.setToolTip('Sends "{}" to the clipboard.'.format(str(tsd.date()))) action = m.addAction('Sensor') action.triggered.connect(lambda: QApplication.clipboard().setText(tsd.sensor().name())) - action.setToolTip('Sends "" to the clipboard.'.format(tsd.sensor().name())) + action.setToolTip('Sends "{}" to the clipboard.'.format(tsd.sensor().name())) action = m.addAction('Path') - action.triggered.connect(lambda: QApplication.clipboard().setText('\n'.join(tsd.sourceUris()))) - action.setToolTip('Sends the {} source URI(s) to the clipboard.'.format(len(tsd))) + + paths = [QDir.toNativeSeparators(p) for p in tsd.sourceUris()] + + action.triggered.connect(lambda _, paths=paths: QApplication.clipboard().setText('\n'.join(paths))) + action.setToolTip('Sends {} source URI(s) to the clipboard.'.format(len(tsd))) action = m.addAction('Map') action.triggered.connect(lambda: QApplication.clipboard().setPixmap(self.pixmap())) action.setToolTip('Copies this map into the clipboard.') + from .utils import findParent + from .mapvisualization import MapWidget + mw = findParent(self, MapWidget) + if isinstance(mw, MapWidget): + action = m.addAction('All Maps') + action.triggered.connect(lambda: QApplication.clipboard().setPixmap(mw.grab())) + action.setToolTip('Copies all maps into the clipboard.') + + m = menu.addMenu('Map Coordinates...') ext = self.spatialExtent() diff --git a/eotimeseriesviewer/mapvisualization.py b/eotimeseriesviewer/mapvisualization.py index 61fa1c80e0789885e840a1b685c576d7c7aa7487..dba3a04d9f27b5bd734cf9ae6a429b5592a0c048 100644 --- a/eotimeseriesviewer/mapvisualization.py +++ b/eotimeseriesviewer/mapvisualization.py @@ -626,6 +626,7 @@ class MapView(QFrame, loadUIFormClass(jp(DIR_UI, 'mapview.ui'))): """ assert isinstance(sensor, SensorInstrument) if sensor not in self.sensors(): + sensor.sigNameChanged.connect(self.sigCanvasAppearanceChanged) dummyLayer = sensor.proxyLayer() assert isinstance(dummyLayer.renderer(), QgsRasterRenderer) dummyLayer.rendererChanged.connect(lambda sensor=sensor: self.onSensorRendererChanged(sensor)) @@ -1054,6 +1055,11 @@ class MapWidget(QFrame, loadUIFormClass(jp(DIR_UI, 'mapwidget.ui'))): n = len(self.timeSeries()) self.mTimeSlider.setRange(0, n) self.mTimeSlider.setEnabled(n > 0) + if n > 10: + pageStep = int(n/100)*10 + else: + pageStep = 5 + self.timeSlider().setPageStep(pageStep) if n > 0: tsd = self.currentDate() @@ -1657,7 +1663,6 @@ class MapViewDock(QgsDockWidget, loadUI('mapviewdock.ui')): self.actionAddMapView.triggered.connect(self.createMapView) self.actionRemoveMapView.triggered.connect(lambda: self.removeMapView(self.currentMapView()) if self.currentMapView() else None) - self.toolBox.currentChanged.connect(self.onToolboxIndexChanged) self.spinBoxMapSizeX.valueChanged.connect(lambda: self.onMapSizeChanged('X')) @@ -1672,8 +1677,6 @@ class MapViewDock(QgsDockWidget, loadUI('mapviewdock.ui')): self.sigMapSizeChanged.emit(QSize(self.spinBoxMapSizeX.value(), self.spinBoxMapSizeY.value())) self.sigMapsPerMapViewChanged.emit(self.mapsPerMapView()) - - def setMapWidget(self, mw)->MapWidget: """ Connects this MapViewDock with a MapWidget diff --git a/eotimeseriesviewer/pixelloader.py b/eotimeseriesviewer/pixelloader.py index e80105ca476917d27255226d07d3fd4163dbf3b5..48fea99e181a70ac1dc7335d17c56ae82613c44c 100644 --- a/eotimeseriesviewer/pixelloader.py +++ b/eotimeseriesviewer/pixelloader.py @@ -38,6 +38,9 @@ if DEBUG: logger = multiprocessing.log_to_stderr() logger.setLevel(multiprocessing.SUBDEBUG) +class TaskMock(QgsTask): + def __init__(self): + super(TaskMock, self).__init__() def dprint(msg): if DEBUG: @@ -186,171 +189,173 @@ def transformPoint2Px(trans, pt, gt): def doLoaderTask(taskWrapper:QgsTask, dump): - # assert isinstance(taskWrapper, QgsTask) - if isinstance(dump, PixelLoaderTask): - task = dump - else: - task = PixelLoaderTask.fromDump(dump) - assert isinstance(task, PixelLoaderTask) - result = task - ds = gdal.Open(task.sourcePath, gdal.GA_ReadOnly) - nb, ns, nl = ds.RasterCount, ds.RasterXSize, ds.RasterYSize - - bandIndices = list(range(nb)) if task.bandIndices is None else list(task.bandIndices) - # ensure to load valid indices only - bandIndices = [i for i in bandIndices if i >= 0 and i < nb] - - task.bandIndices = bandIndices - - gt = ds.GetGeoTransform() - result.resGeoTransformation = gt - result.resCrsWkt = ds.GetProjection() - if result.resCrsWkt == '': - result.resCrsWkt = QgsRasterLayer(ds.GetDescription()).crs().toWkt() - crsSrc = osr.SpatialReference(result.resCrsWkt) - - # convert Geometries into pixel indices to be extracted - PX_SUBSETS = [] - - for geom in task.geometries: - crsRequest = osr.SpatialReference() - - if geom.crs().isValid(): - crsRequest.ImportFromWkt(geom.crs().toWkt()) - else: - crsRequest.ImportFromWkt(crsSrc.ExportToWkt()) - trans = osr.CoordinateTransformation(crsRequest, crsSrc) - - if isinstance(geom, QgsPointXY): - ptUL = ptLR = QgsPointXY(geom) - elif isinstance(geom, QgsRectangle): - TYPE = 'RECTANGLE' - ptUL = QgsPointXY(geom.xMinimum(), geom.yMaximum()) - ptLR = QgsPointXY(geom.xMaximum(), geom.yMinimum()) - else: - raise NotImplementedError('Unsupported geometry {} {}'.format(type(geom), str(geom))) + assert isinstance(taskWrapper, QgsTask) - pxUL = transformPoint2Px(trans, ptUL, gt) - pxLR = transformPoint2Px(trans, ptLR, gt) + tasks = pickle.loads(dump) + results = [] + for task in tasks: + assert isinstance(task, PixelLoaderTask) - bUL = isOutOfImage(ds, pxUL) - bLR = isOutOfImage(ds, pxLR) + result = task + ds = gdal.Open(task.sourcePath, gdal.GA_ReadOnly) + nb, ns, nl = ds.RasterCount, ds.RasterXSize, ds.RasterYSize - if all([bUL, bLR]): - PX_SUBSETS.append(INFO_OUT_OF_IMAGE) - continue + bandIndices = list(range(nb)) if task.bandIndices is None else list(task.bandIndices) + # ensure to load valid indices only + bandIndices = [i for i in bandIndices if i >= 0 and i < nb] + task.bandIndices = bandIndices - def shiftIntoImageBounds(pt, xMax, yMax): - assert isinstance(pt, QPoint) - if pt.x() < 0: - pt.setX(0) - elif pt.x() > xMax: - pt.setX(xMax) - if pt.y() < 0: - pt.setY(0) - elif pt.y() > yMax: - pt.setY(yMax) + gt = ds.GetGeoTransform() + result.resGeoTransformation = gt + result.resCrsWkt = ds.GetProjection() + if result.resCrsWkt == '': + result.resCrsWkt = QgsRasterLayer(ds.GetDescription()).crs().toWkt() + crsSrc = osr.SpatialReference(result.resCrsWkt) + # convert Geometries into pixel indices to be extracted + PX_SUBSETS = [] - shiftIntoImageBounds(pxUL, ds.RasterXSize, ds.RasterYSize) - shiftIntoImageBounds(pxLR, ds.RasterXSize, ds.RasterYSize) + for geom in task.geometries: + crsRequest = osr.SpatialReference() - if pxUL == pxLR: - size_x = size_y = 1 - else: - size_x = abs(pxUL.x() - pxLR.x()) - size_y = abs(pxUL.y() - pxLR.y()) - - if size_x < 1: size_x = 1 - if size_y < 1: size_y = 1 + if geom.crs().isValid(): + crsRequest.ImportFromWkt(geom.crs().toWkt()) + else: + crsRequest.ImportFromWkt(crsSrc.ExportToWkt()) + trans = osr.CoordinateTransformation(crsRequest, crsSrc) + + if isinstance(geom, QgsPointXY): + ptUL = ptLR = QgsPointXY(geom) + elif isinstance(geom, QgsRectangle): + TYPE = 'RECTANGLE' + ptUL = QgsPointXY(geom.xMinimum(), geom.yMaximum()) + ptLR = QgsPointXY(geom.xMaximum(), geom.yMinimum()) + else: + raise NotImplementedError('Unsupported geometry {} {}'.format(type(geom), str(geom))) - PX_SUBSETS.append((pxUL, pxUL, size_x, size_y)) + pxUL = transformPoint2Px(trans, ptUL, gt) + pxLR = transformPoint2Px(trans, ptLR, gt) - PROFILE_DATA = [] + bUL = isOutOfImage(ds, pxUL) + bLR = isOutOfImage(ds, pxLR) - if bandIndices == range(ds.RasterCount): - # we have to extract all bands - # in this case we use gdal.Dataset.ReadAsArray() - noData = ds.GetRasterBand(1).GetNoDataValue() - for px in PX_SUBSETS: - if px == INFO_OUT_OF_IMAGE: - PROFILE_DATA.append(INFO_OUT_OF_IMAGE) + if all([bUL, bLR]): + PX_SUBSETS.append(INFO_OUT_OF_IMAGE) continue - pxUL, pxUL, size_x, size_y = px - - bandData = ds.ReadAsArray(pxUL.x(), pxUL.y(), size_x, size_y).reshape((nb, size_x*size_y)) - if noData: - isValid = np.ones(bandData.shape[1], dtype=np.bool) - for b in range(bandData.shape[0]): - isValid *= bandData[b, :] != ds.GetRasterBand(b+1).GetNoDataValue() - bandData = bandData[:, np.where(isValid)[0]] - PROFILE_DATA.append(bandData) - else: - # access band values band-by-band - # in this case we use gdal.Band.ReadAsArray() - # and need to iterate over the requested band indices - - # save the returned band values for each geometry in a separate list - # empty list == invalid geometry - for i in range(len(PX_SUBSETS)): - if PX_SUBSETS[i] == INFO_OUT_OF_IMAGE: - PROFILE_DATA.append(INFO_OUT_OF_IMAGE) + + def shiftIntoImageBounds(pt, xMax, yMax): + assert isinstance(pt, QPoint) + if pt.x() < 0: + pt.setX(0) + elif pt.x() > xMax: + pt.setX(xMax) + if pt.y() < 0: + pt.setY(0) + elif pt.y() > yMax: + pt.setY(yMax) + + + shiftIntoImageBounds(pxUL, ds.RasterXSize, ds.RasterYSize) + shiftIntoImageBounds(pxLR, ds.RasterXSize, ds.RasterYSize) + + if pxUL == pxLR: + size_x = size_y = 1 else: - PROFILE_DATA.append([]) + size_x = abs(pxUL.x() - pxLR.x()) + size_y = abs(pxUL.y() - pxLR.y()) + if size_x < 1: size_x = 1 + if size_y < 1: size_y = 1 - for bandIndex in bandIndices: - band = ds.GetRasterBand(bandIndex+1) - noData = band.GetNoDataValue() - assert isinstance(band, gdal.Band) + PX_SUBSETS.append((pxUL, pxUL, size_x, size_y)) - for i, px in enumerate(PX_SUBSETS): + PROFILE_DATA = [] + + if bandIndices == range(ds.RasterCount): + # we have to extract all bands + # in this case we use gdal.Dataset.ReadAsArray() + noData = ds.GetRasterBand(1).GetNoDataValue() + for px in PX_SUBSETS: if px == INFO_OUT_OF_IMAGE: + PROFILE_DATA.append(INFO_OUT_OF_IMAGE) continue + pxUL, pxUL, size_x, size_y = px - bandData = band.ReadAsArray(pxUL.x(), pxUL.y(), size_x, size_y).flatten() + + bandData = ds.ReadAsArray(pxUL.x(), pxUL.y(), size_x, size_y).reshape((nb, size_x*size_y)) if noData: - bandData = bandData[np.where(bandData != noData)[0]] - PROFILE_DATA[i].append(bandData) - - for i in range(len(PX_SUBSETS)): - pd = PROFILE_DATA[i] - if isinstance(pd, list): - if len(pd) == 0: - PROFILE_DATA[i] = INFO_OUT_OF_IMAGE + isValid = np.ones(bandData.shape[1], dtype=np.bool) + for b in range(bandData.shape[0]): + isValid *= bandData[b, :] != ds.GetRasterBand(b+1).GetNoDataValue() + bandData = bandData[:, np.where(isValid)[0]] + PROFILE_DATA.append(bandData) + else: + # access band values band-by-band + # in this case we use gdal.Band.ReadAsArray() + # and need to iterate over the requested band indices + + # save the returned band values for each geometry in a separate list + # empty list == invalid geometry + for i in range(len(PX_SUBSETS)): + if PX_SUBSETS[i] == INFO_OUT_OF_IMAGE: + PROFILE_DATA.append(INFO_OUT_OF_IMAGE) else: - #PROFILE_DATA[i] = np.dstack(pd).transpose(2,0,1) - PROFILE_DATA[i] = np.vstack(pd) + PROFILE_DATA.append([]) + + + for bandIndex in bandIndices: + band = ds.GetRasterBand(bandIndex+1) + noData = band.GetNoDataValue() + assert isinstance(band, gdal.Band) + + for i, px in enumerate(PX_SUBSETS): + if px == INFO_OUT_OF_IMAGE: + continue + pxUL, pxUL, size_x, size_y = px + bandData = band.ReadAsArray(pxUL.x(), pxUL.y(), size_x, size_y).flatten() + if noData: + bandData = bandData[np.where(bandData != noData)[0]] + PROFILE_DATA[i].append(bandData) + + for i in range(len(PX_SUBSETS)): + pd = PROFILE_DATA[i] + if isinstance(pd, list): + if len(pd) == 0: + PROFILE_DATA[i] = INFO_OUT_OF_IMAGE + else: + #PROFILE_DATA[i] = np.dstack(pd).transpose(2,0,1) + PROFILE_DATA[i] = np.vstack(pd) - # finally, ensure that there is on 2D array only - for i in range(len(PROFILE_DATA)): - d = PROFILE_DATA[i] - if isinstance(d, np.ndarray): - assert d.ndim == 2 - b, yx = d.shape - assert b == len(bandIndices) + # finally, ensure that there is on 2D array only + for i in range(len(PROFILE_DATA)): + d = PROFILE_DATA[i] + if isinstance(d, np.ndarray): + assert d.ndim == 2 + b, yx = d.shape + assert b == len(bandIndices) - _, _, size_x, size_y = PX_SUBSETS[i] + _, _, size_x, size_y = PX_SUBSETS[i] - if yx > 0: - d = d.reshape((b, yx)) - vMean = d.mean(axis=1) - vStd = d.std(axis=1) + if yx > 0: + d = d.reshape((b, yx)) + vMean = d.mean(axis=1) + vStd = d.std(axis=1) - assert len(vMean) == len(bandIndices) - assert len(vStd) == len(bandIndices) - PROFILE_DATA[i] = (vMean, vStd) - else: - PROFILE_DATA[i] = INFO_NO_DATA + assert len(vMean) == len(bandIndices) + assert len(vStd) == len(bandIndices) + PROFILE_DATA[i] = (vMean, vStd) + else: + PROFILE_DATA[i] = INFO_NO_DATA - s = "" - task.resProfiles = PROFILE_DATA - task.mIsDone = True - return task.toDump() + s = "" + task.resProfiles = PROFILE_DATA + task.mIsDone = True + results.append(task) + return pickle.dumps(results) @@ -416,31 +421,38 @@ class PixelLoader(QObject): assert isinstance(tasks, list) - tm = self.taskManager() + dump = pickle.dumps(tasks) - #self.sigLoadingStarted.emit() - #todo: create chuncks - import uuid - for plt in tasks: - assert isinstance(plt, PixelLoaderTask) - taskName = 'pltTask.{}'.format(uuid.uuid4()) - plt.setId(taskName) - dump = plt.toDump() - qgsTask = QgsTask.fromFunction(taskName, doLoaderTask, dump, on_finished=self.onLoadingFinished) + if False: + qgsTask = TaskMock() + resultDump = doLoaderTask(qgsTask, dump) + self.onLoadingFinished(TaskMock(), resultDump) + + else: + qgsTask = QgsTask.fromFunction('Load band values', doLoaderTask, dump, on_finished=self.onLoadingFinished) + tid = id(qgsTask) + qgsTask.taskCompleted.connect(lambda *args, tid=tid: self.onRemoveTask(tid)) + qgsTask.taskTerminated.connect(lambda *args, tid=tid: self.onRemoveTask(tid)) + # todo: progress changed? + self.mTasks[tid] = qgsTask + tm = self.taskManager() tm.addTask(qgsTask) - self.mTasks[taskName] = qgsTask + def onRemoveTask(self, tid): + if tid in self.mTasks.keys(): + del self.mTasks[tid] def onLoadingFinished(self, *args, **kwds): error = args[0] if error is None: dump = args[1] - plt = PixelLoaderTask.fromDump(dump) - if isinstance(plt, PixelLoaderTask): - self.mTasks.pop(plt.id()) + tasks = pickle.loads(dump) + for plt in tasks: + assert isinstance(plt, PixelLoaderTask) + #self.mTasks.pop(plt.id()) self.sigPixelLoaded.emit(plt) diff --git a/eotimeseriesviewer/settings.py b/eotimeseriesviewer/settings.py index 881762d22edd16dda96c62dd13a94191dfa1370c..81b5bf9e17fe4b89d785d0713e3f9565a7f9855e 100644 --- a/eotimeseriesviewer/settings.py +++ b/eotimeseriesviewer/settings.py @@ -48,14 +48,14 @@ def defaultValues() -> dict: d[Keys.MapBackgroundColor] = QColor('black') textFormat = QgsTextFormat() - textFormat.setColor(QColor('yellow')) + textFormat.setColor(QColor('black')) textFormat.setSizeUnit(QgsUnitTypes.RenderPoints) textFormat.setFont(QFont('Helvetica')) textFormat.setSize(11) buffer = QgsTextBufferSettings() - buffer.setColor(QColor('black')) - buffer.setSize(1) + buffer.setColor(QColor('white')) + buffer.setSize(5) buffer.setSizeUnit(QgsUnitTypes.RenderPixels) buffer.setEnabled(True) textFormat.setBuffer(buffer) diff --git a/eotimeseriesviewer/timeseries.py b/eotimeseriesviewer/timeseries.py index 6214e08222961b05be2b13c07d1eef67b9af7f33..f318a53b706c7562d7cf58e405ed0aaf597b8f5a 100644 --- a/eotimeseriesviewer/timeseries.py +++ b/eotimeseriesviewer/timeseries.py @@ -21,8 +21,9 @@ # noinspection PyPep8Naming import sys, re, collections, traceback, time, json, urllib, types, enum, typing, pickle, json, uuid +from xml.etree import ElementTree - +from qgis.PyQt.QtXml import QDomDocument import bisect from qgis import * @@ -33,6 +34,17 @@ from qgis.PyQt.QtWidgets import * from qgis.PyQt.QtCore import * +LUT_WAVELENGTH_UNITS = {} +for siUnit in [r'nm', r'μm', r'mm', r'cm', r'dm']: + LUT_WAVELENGTH_UNITS[siUnit] = siUnit +LUT_WAVELENGTH_UNITS[r'nanometers'] = r'nm' +LUT_WAVELENGTH_UNITS[r'micrometers'] = r'μm' +LUT_WAVELENGTH_UNITS[r'um'] = r'μm' +LUT_WAVELENGTH_UNITS[r'millimeters'] = r'mm' +LUT_WAVELENGTH_UNITS[r'centimeters'] = r'cm' +LUT_WAVELENGTH_UNITS[r'decimeters'] = r'dm' + + from osgeo import gdal from eotimeseriesviewer.dateparser import DOYfromDatetime64 from eotimeseriesviewer.utils import SpatialExtent, loadUI, px2geo @@ -879,15 +891,43 @@ class TimeSeriesTreeView(QTreeView): """ idx = self.indexAt(event.pos()) - tsd = self.model().data(idx, role=Qt.UserRole) + node = self.model().data(idx, role=Qt.UserRole) menu = QMenu(self) - if isinstance(tsd, TimeSeriesDate): - a = menu.addAction('Show map for {}'.format(tsd.date())) + + def setUri(paths): + urls = [] + paths2 = [] + for p in paths: + if os.path.isfile(p): + url = QUrl.fromLocalFile(p) + paths2.append(QDir.toNativeSeparators(p)) + else: + url = QUrl(p) + paths2.append(p) + urls.append(url) + md = QMimeData() + md.setText('\n'.join(paths2)) + md.setUrls(urls) + + QApplication.clipboard().setMimeData(md) + + + if isinstance(node, TimeSeriesDate): + a = menu.addAction('Show map for {}'.format(node.date())) a.setToolTip('Shows the map related to this time series date.') - a.triggered.connect(lambda _, tsd=tsd: self.sigMoveToDateRequest.emit(tsd)) + a.triggered.connect(lambda _, tsd=node: self.sigMoveToDateRequest.emit(tsd)) menu.addSeparator() + a = menu.addAction('Copy path(s)') + a.triggered.connect(lambda _, paths=node.sourceUris(): setUri(paths)) + a.setToolTip('Copy path to cliboard') + + if isinstance(node, TimeSeriesSource): + a = menu.addAction('Copy path') + a.triggered.connect(lambda _, paths=[node.uri()]: setUri(paths)) + a.setToolTip('Copy path to cliboard') + a = menu.addAction('Copy value(s)') a.triggered.connect(lambda: self.onCopyValues()) a = menu.addAction('Hide date(s)') @@ -1869,69 +1909,180 @@ def getSpatialPropertiesFromDataset(ds): return nb, nl, ns, crs, px_x, px_y +def extractWavelengthsFromGDALMetaData(ds:gdal.Dataset)->(list, str): + """ + Reads the wavelength info from standard metadata strings + :param ds: gdal.Dataset + :return: (list, str) + """ + + regWLkey = re.compile('^(center )?wavelength[_ ]*$', re.I) + regWLUkey = re.compile('^wavelength[_ ]*units?$', re.I) + regNumeric = re.compile(r"([-+]?\d*\.\d+|[-+]?\d+)", re.I) + def findKey(d:dict, regex)->str: + for key in d.keys(): + if regex.search(key): + return key + # 1. try band level + wlu = [] + wl = [] + for b in range(ds.RasterCount): + band = ds.GetRasterBand(b + 1) + assert isinstance(band, gdal.Band) + md = band.GetMetadata_Dict() + keyWLU = findKey(md, regWLUkey) + keyWL = findKey(md, regWLkey) + if isinstance(keyWL, str) and isinstance(keyWLU, str): + wl.append(float(md[keyWL])) + wlu.append(LUT_WAVELENGTH_UNITS[md[keyWLU].lower()]) + if len(wlu) == len(wl) and len(wl) == ds.RasterCount: + return wl, wlu[0] + # 2. try data set level + for domain in ds.GetMetadataDomainList(): + md = ds.GetMetadata_Dict(domain) + keyWLU = findKey(md, regWLUkey) + keyWL = findKey(md, regWLkey) + if isinstance(keyWL, str) and isinstance(keyWLU, str): -def extractWavelengths(ds): - wl = None - wlu = None + wlu = LUT_WAVELENGTH_UNITS[md[keyWLU].lower()] + matches = regNumeric.findall(md[keyWL]) + wl = [float(n) for n in matches] - # see http://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html for supported wavelength units - regWLkey = re.compile('.*wavelength[_ ]*$', re.I) - regWLUkey = re.compile('.*wavelength[_ ]*units?$', re.I) - regNumeric = re.compile(r"([-+]?\d*\.\d+|[-+]?\d+)", re.I) - regWLU = re.compile('((micro|nano|centi)meters)|(um|nm|mm|cm|m|GHz|MHz)', re.I) + + + if len(wl) == ds.RasterCount: + return wl, wlu + + return None, None + + + +def extractWavelengthsFromRapidEyeXML(ds:gdal.Dataset, dom:QDomDocument)->(list, str): + nodes = dom.elementsByTagName('re:bandSpecificMetadata') + # see http://schemas.rapideye.de/products/re/4.0/RapidEye_ProductMetadata_GeocorrectedLevel.xsd + # wavelength and units not given in the XML + # -> use values from https://www.satimagingcorp.com/satellite-sensors/other-satellite-sensors/rapideye/ + if nodes.count() == ds.RasterCount and ds.RasterCount == 5: + wlu = r'nm' + wl = [0.5 * (440 + 510), + 0.5 * (520 + 590), + 0.5 * (630 + 685), + 0.5 * (760 + 850), + 0.5 * (760 - 850) + ] + return wl, wlu + return None, None + + +def extractWavelengthsFromDIMAPXML(ds:gdal.Dataset, dom:QDomDocument)->(list, str): + """ + :param dom: QDomDocument | gdal.Dataset + :return: (list of wavelengths, str wavelength unit) + """ + # DIMAP XML metadata? + assert isinstance(dom, QDomDocument) + nodes = dom.elementsByTagName('Band_Spectral_Range') + if nodes.count() > 0: + candidates = [] + for element in [nodes.item(i).toElement() for i in range(nodes.count())]: + _band = element.firstChildElement('BAND_ID').text() + _wlu = element.firstChildElement('MEASURE_UNIT').text() + wlMin = float(element.firstChildElement('MIN').text()) + wlMax = float(element.firstChildElement('MAX').text()) + _wl = 0.5 * wlMin + wlMax + candidates.append((_band, _wl, _wlu)) + + if len(candidates) == ds.RasterCount: + candidates = sorted(candidates, key=lambda t: t[0]) + + wlu = candidates[0][2] + wlu = LUT_WAVELENGTH_UNITS[wlu] + wl = [c[1] for c in candidates] + return wl, wlu + return None, None + +def extractWavelengths(ds): + """ + Returns the wavelength and wavelength units + :param ds: gdal.Dataset + :return: (float [list-of-wavelengths], str with wavelength unit) + """ if isinstance(ds, QgsRasterLayer): - lyr = ds - md = [l.split('=') for l in str(lyr.metadata()).splitlines() if 'wavelength' in l.lower()] - #see http://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html for supported wavelength units - for kv in md: - key, value = kv - key = key.lower() - if key == 'center wavelength': - tmp = re.findall(r'\d*\.\d+|\d+', value) #find floats - if len(tmp) == 0: - tmp = re.findall(r'\d+', value) #find integers - if len(tmp) == lyr.bandCount(): - wl = [float(w) for w in tmp] - - if key == 'wavelength units': - match = regWLU.search(value) - if match: - wlu = match.group() - - names = ['nanometers','micrometers','millimeters','centimeters','decimenters'] - si = ['nm','um','mm','cm','dm'] - if wlu in names: - wlu = si[names.index(wlu)] + + if ds.dataProvider().name() == 'gdal': + uri = ds.source() + return extractWavelengths(gdal.Open(uri)) + else: + + md = [l.split('=') for l in str(ds.metadata()).splitlines() if 'wavelength' in l.lower()] + + wl = wlu = None + for kv in md: + key, value = kv + key = key.lower() + value = value.strip() + + if key == 'wavelength': + tmp = re.findall(r'\d*\.\d+|\d+', value) #find floats + if len(tmp) == 0: + tmp = re.findall(r'\d+', value) #find integers + if len(tmp) == ds.bandCount(): + wl = [float(w) for w in tmp] + + if key == 'wavelength units': + wlu = value + if wlu in LUT_WAVELENGTH_UNITS.keys(): + wlu = LUT_WAVELENGTH_UNITS[wlu] + + if isinstance(wl, list) and isinstance(wlu, str): + return wl, wlu + elif isinstance(ds, gdal.Dataset): - for domain in ds.GetMetadataDomainList(): - md = ds.GetMetadata_Dict(domain) - for key, value in md.items(): - if wl is None and regWLkey.search(key): - numbers = regNumeric.findall(value) - if len(numbers) == ds.RasterCount: - wl = [float(n) for n in numbers] - - if wlu is None and regWLUkey.search(key): - match = regWLU.search(value) - if match: - wlu = match.group().lower() - names = ['nanometers', 'micrometers', 'millimeters', 'centimeters', 'decimeters'] - si = ['nm', 'um', 'mm', 'cm', 'dm'] - if wlu in names: - wlu = si[names.index(wlu)] - - return wl, wlu + def testWavelLengthInfo(wl, wlu)->bool: + return isinstance(wl, list) and len(wl) == ds.RasterCount and isinstance(wlu, str) and wlu in LUT_WAVELENGTH_UNITS.keys() + + # try band-specific metadata + wl, wlu = extractWavelengthsFromGDALMetaData(ds) + if testWavelLengthInfo(wl, wlu): + return wl, wlu + + # try internal locations with XML info + # SPOT DIMAP + if 'xml:dimap' in ds.GetMetadataDomainList(): + md = ds.GetMetadata_Dict('xml:dimap') + for key in md.keys(): + dom = QDomDocument() + dom.setContent(key + '=' + md[key]) + wl, wlu = extractWavelengthsFromDIMAPXML(ds, dom) + if testWavelLengthInfo(wl, wlu): + return wl, wlu + + # try separate XML files + xmlReaders = [extractWavelengthsFromDIMAPXML, extractWavelengthsFromRapidEyeXML] + for path in ds.GetFileList(): + if re.search(r'\.xml$', path, re.I) and not re.search(r'\.aux.xml$', path, re.I): + dom = QDomDocument() + with open(path, encoding='utf-8') as f: + dom.setContent(f.read()) + + if dom.hasChildNodes(): + for xmlReader in xmlReaders: + wl, wlu = xmlReader(ds, dom) + if testWavelLengthInfo(wl, wlu): + return wl, wlu + + return None, None diff --git a/make/createtestdata.py b/make/createtestdata.py index a6a5f681897099ef8b005ac8c4b1ab8a507b16ca..e398f204d69a6b9b890722ef3dd1e6901b57e63a 100644 --- a/make/createtestdata.py +++ b/make/createtestdata.py @@ -136,7 +136,7 @@ def groupLandsat(dirIn, dirOut, pattern='L*_sr_band*.img'): raise NotImplementedError() #https://www.harrisgeospatial.com/docs/ENVIHeaderFiles.html - dsVRT.SetMetadataItem('wavelength units','Micrometers', 'ENVI') + dsVRT.SetMetadataItem('wavelength units', 'Micrometers', 'ENVI') dsVRT.SetMetadataItem('wavelength', '{{{}}}'.format(','.join([str(w) for w in cwl])), 'ENVI') dsVRT.SetMetadataItem('sensor type', 'Landsat-8 OLI', 'ENVI') from eotimeseriesviewer.dateparser import datetime64FromYYYYDOY diff --git a/tests/test_pixelloader.py b/tests/test_pixelloader.py index 8ec7724262834483095a4fe05b25a02b7e456400..a3fa6583ad3939c3507ae83e7ed0ccaca8e86fe1 100644 --- a/tests/test_pixelloader.py +++ b/tests/test_pixelloader.py @@ -175,36 +175,38 @@ class PixelLoaderTest(unittest.TestCase): #test a valid pixels plt = PixelLoaderTask(source, [ptValid1, ptValid2]) + dump = pickle.dumps([plt]) task = QgsTask.fromFunction('', doLoaderTask, plt) try: - result = doLoaderTask(task, plt.toDump()) - result = PixelLoaderTask.fromDump(result) + results = doLoaderTask(task, dump) + results = pickle.loads(results) except Exception as ex: self.fail('Failed to return the pixels for two geometries') + for result in results: - self.assertIsInstance(result, PixelLoaderTask) - self.assertTrue(result.success()) - self.assertEqual(result.sourcePath, source) - self.assertSequenceEqual(result.bandIndices, [0,1,2,3,4,5]) - self.assertIs(result.exception, None) + self.assertIsInstance(result, PixelLoaderTask) + self.assertTrue(result.success()) + self.assertEqual(result.sourcePath, source) + self.assertSequenceEqual(result.bandIndices, [0,1,2,3,4,5]) + self.assertIs(result.exception, None) - self.assertEqual(result.resCrsWkt, wkt) - self.assertEqual(result.resGeoTransformation, gt) - self.assertEqual(result.resNoDataValues, None) - self.assertIsInstance(result.resProfiles, list) - self.assertEqual(len(result.resProfiles), 2, msg='did not return results for two geometries.') + self.assertEqual(result.resCrsWkt, wkt) + self.assertEqual(result.resGeoTransformation, gt) + self.assertEqual(result.resNoDataValues, None) + self.assertIsInstance(result.resProfiles, list) + self.assertEqual(len(result.resProfiles), 2, msg='did not return results for two geometries.') - for i in range(len(result.resProfiles)): - dn, stdDev = result.resProfiles[i] - self.assertIsInstance(dn, np.ndarray) - self.assertIsInstance(stdDev, np.ndarray) + for i in range(len(result.resProfiles)): + dn, stdDev = result.resProfiles[i] + self.assertIsInstance(dn, np.ndarray) + self.assertIsInstance(stdDev, np.ndarray) - self.assertEqual(len(dn), len(stdDev)) - self.assertEqual(len(dn), 6) - self.assertEqual(stdDev.max(), 0) #std deviation of a single pixel is always zero - self.assertEqual(stdDev.min(), 0) + self.assertEqual(len(dn), len(stdDev)) + self.assertEqual(len(dn), 6) + self.assertEqual(stdDev.max(), 0) #std deviation of a single pixel is always zero + self.assertEqual(stdDev.min(), 0) diff --git a/tests/test_timeseries.py b/tests/test_timeseries.py index 760fbf98635a5e0e48e0747e597d2307fc023ed9..f40590f37cd2bf6ba63ba6877b6709a9d758e436 100644 --- a/tests/test_timeseries.py +++ b/tests/test_timeseries.py @@ -1,4 +1,4 @@ -# coding=utf-8 +# -*- coding: utf-8 -*- """Tests QGIS plugin init.""" import os @@ -13,7 +13,7 @@ from eotimeseriesviewer.tests import initQgisApplication QAPP = initQgisApplication() -SHOW_GUI = True and os.environ.get('CI') is None +SHOW_GUI = False and os.environ.get('CI') is None class TestInit(unittest.TestCase): @@ -87,6 +87,7 @@ class TestInit(unittest.TestCase): c2 = sensorIDtoProperties(sid) self.assertListEqual(list(conf), list(c2)) + s = "" def test_TimeSeriesDate(self): @@ -318,6 +319,82 @@ class TestInit(unittest.TestCase): self.assertIsInstance(extent, SpatialExtent) + def test_pleiades(self): + + paths = [r'Y:\Pleiades\GFIO_Gp13_Novo_SO16018091-4-01_DS_PHR1A_201703031416139_FR1_PX_W056S07_0906_01636\TPP1600581943\IMG_PHR1A_PMS_001\DIM_PHR1A_PMS_201703031416139_ORT_2224693101-001.XML' + ,r'Y:\Pleiades\GFIO_Gp13_Novo_SO16018091-4-01_DS_PHR1A_201703031416139_FR1_PX_W056S07_0906_01636\TPP1600581943\IMG_PHR1A_PMS_001\IMG_PHR1A_PMS_201703031416139_ORT_2224693101-001_R1C1.JP2' + ] + for p in paths: + if not os.path.isfile(p): + continue + + ds = gdal.Open(p) + self.assertIsInstance(ds, gdal.Dataset) + band = ds.GetRasterBand(1) + self.assertIsInstance(band, gdal.Band) + + + tss = TimeSeriesSource(ds) + self.assertIsInstance(tss, TimeSeriesSource) + self.assertEqual(tss.mWLU, r'μm') + self.assertListEqual(tss.mWL, [0.775, 0.867, 1.017, 1.315]) + + s = "" + + def test_rapideye(self): + from example.Images import re_2014_06_25 + paths = [r'Y:\RapidEye\3A\2135821_2014-06-25_RE2_3A_328202\2135821_2014-06-25_RE2_3A_328202.tif'] + + for p in paths: + if not os.path.isfile(p): + continue + + ds = gdal.Open(p) + self.assertIsInstance(ds, gdal.Dataset) + band = ds.GetRasterBand(1) + self.assertIsInstance(band, gdal.Band) + + + tss = TimeSeriesSource(ds) + self.assertIsInstance(tss, TimeSeriesSource) + + # see https://www.satimagingcorp.com/satellite-sensors/other-satellite-sensors/rapideye/ + wlu = r'nm' + wl = [0.5 * (440 + 510), + 0.5 * (520 + 590), + 0.5 * (630 + 685), + 0.5 * (760 + 850), + 0.5 * (760 - 850) + ] + self.assertEqual(tss.mWLU, wlu) + self.assertListEqual(tss.mWL, wl) + + def test_sentinel2(self): + + p = r'Q:\Processing_BJ\01_Data\Sentinel\T21LXL\S2A_MSIL1C_20161221T141042_N0204_R110_T21LXL_20161221T141040.SAFE\MTD_MSIL1C.xml' + + if not os.path.isfile(p): + return + + dsC = gdal.Open(p) + self.assertIsInstance(dsC, gdal.Dataset) + for item in dsC.GetSubDatasets(): + path = item[0] + ds = gdal.Open(path) + gt = ds.GetGeoTransform() + self.assertIsInstance(ds, gdal.Dataset) + + band = ds.GetRasterBand(1) + self.assertIsInstance(band, gdal.Band) + + wlu = ds.GetRasterBand(1).GetMetadata_Dict()['WAVELENGTH_UNIT'] + wl = [float(ds.GetRasterBand(b+1).GetMetadata_Dict()['WAVELENGTH']) for b in range(ds.RasterCount)] + + tss = TimeSeriesSource(ds) + self.assertIsInstance(tss, TimeSeriesSource) + + self.assertEqual(tss.mWLU, wlu) + self.assertEqual(tss.mWL, wl) def test_sensors(self):