Newer
Older
n = len(sources)
for i, source in enumerate(sources):
if qgsTask.isCanceled():
return pickle.dumps(results)
s = TimeSeriesSource.create(source)
if isinstance(s, TimeSeriesSource):
results.append(s)
qgsTask.setProgress(i + 1)
return pickle.dumps(results)
class TimeSeries(QAbstractItemModel):
"""
The sorted list of data sources that specify the time series
"""
sigTimeSeriesDatesAdded = pyqtSignal(list)
sigTimeSeriesDatesRemoved = pyqtSignal(list)
sigSensorAdded = pyqtSignal(SensorInstrument)
sigSensorRemoved = pyqtSignal(SensorInstrument)
sigSourcesAdded = pyqtSignal(list)
sigSourcesRemoved = pyqtSignal(list)

Benjamin Jakimow
committed
sigVisibilityChanged = pyqtSignal()
_sep = ';'
def __init__(self, imageFiles=None):
super(TimeSeries, self).__init__()
self.mTSDs = list()
self.mSensors = []
self.mShape = None

Benjamin Jakimow
committed
self.mDateTimePrecision = DateTimePrecision.Original
self.mLoadingProgressDialog = None
self.mVisibleDate = []
self.cnDate = 'Date'
self.cnSensor = 'Sensor'
self.cnNS = 'ns'
self.cnNL = 'nl'
self.cnNB = 'nb'
self.cnCRS = 'CRS'
self.cnImages = 'Source Image(s)'
self.mColumnNames = [self.cnDate, self.cnSensor,
self.cnNS, self.cnNL, self.cnNB,
self.cnCRS, self.cnImages]
self.mRootIndex = QModelIndex()
self.mTasks = dict()
def setCurrentSpatialExtent(self, spatialExtent:SpatialExtent):
"""
Sets the spatial extent currently shown
:param spatialExtent:
"""
if isinstance(spatialExtent, SpatialExtent) and self.mCurrentSpatialExtent != spatialExtent:
self.mCurrentSpatialExtent = spatialExtent
def focusVisibilityToExtent(self, ext:SpatialExtent=None):
"""
Changes TSDs visibility according to its intersection with a SpatialExtent
:param ext: SpatialExtent
"""
if ext is None:
ext = self.currentSpatialExtent()
changed = False
assert isinstance(tsd, TimeSeriesDate)
if b != tsd.isVisible():
changed = True
if changed:
ul = self.index(0, 0)
lr = self.index(self.rowCount()-1, 0)
self.sigVisibilityChanged.emit()
def currentSpatialExtent(self)->SpatialExtent:
"""
Returns the current spatial extent
:return: SpatialExtent
"""
return self.mCurrentSpatialExtent
def setVisibleDates(self, tsds:list):
:param tsds: [list-of-TimeSeriesDate]
self.mVisibleDate.clear()
self.mVisibleDate.extend(tsds)
for tsd in tsds:
assert isinstance(tsd, TimeSeriesDate)
if tsd in self:
idx = self.tsdToIdx(tsd)
# force reset of background color
idx2 = self.index(idx.row(), self.columnCount()-1)
self.dataChanged.emit(idx, idx2, [Qt.BackgroundColorRole])
def sensor(self, sensorID:str)->SensorInstrument:
for sensor in self.mSensors:
assert isinstance(sensor, SensorInstrument)
return sensor
return None
def sensors(self)->list:
"""
Returns the list of sensors derived from the TimeSeries data sources
:return: [list-of-SensorInstruments]
"""
def loadFromFile(self, path, n_max=None, progressDialog:QProgressDialog=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 = []
with open(path, 'r') as f:
lines = f.readlines()
for l in lines:
if re.match('^[ ]*[;#&]', l):
continue
parts = re.split('[\n'+TimeSeries._sep+']', l)
parts = [p for p in parts if p != '']
images.append(parts[0])
if len(parts) > 1:
masks.append(parts[1])
images = images[0:n_max]
if isinstance(progressDialog, QProgressDialog):
progressDialog.setMaximum(len(images))
progressDialog.setMinimum(0)
progressDialog.setValue(0)
progressDialog.setLabelText('Start loading {} images....'.format(len(images)))
self.addSources(images, progressDialog=progressDialog)
"""
Saves the TimeSeries sources into a CSV file
:param path: str, path of CSV file
:return: path of CSV file
"""
lines = []
lines.append('#Time series definition file: {}'.format(np.datetime64('now').astype(str)))
assert isinstance(TSD, TimeSeriesDate)
for pathImg in TSD.sourceUris():
lines.append(pathImg)
messageLog('Time series source images written to {}'.format(path))
"""
Returns the pixel sizes of all SensorInstruments
:return: [list-of-QgsRectangles]
"""

benjamin.jakimow@geo.hu-berlin.de
committed
r = []

benjamin.jakimow@geo.hu-berlin.de
committed
r.append((QgsRectangle(sensor.px_size_x, sensor.px_size_y)))
return r
"""
Returns the maximum SpatialExtent of all images of the TimeSeries
:param crs: QgsCoordinateSystem to express the SpatialExtent coordinates.
:return:
"""
assert isinstance(tsd, TimeSeriesDate)
ext = tsd.spatialExtent()
if isinstance(extent, SpatialExtent):
extent = extent.combineExtentWith(ext)
else:
extent = ext

benjamin.jakimow@geo.hu-berlin.de
committed
def getTSD(self, pathOfInterest):
Returns the TimeSeriesDate related to an image source
:param pathOfInterest: str, image source uri
tsd = self.mLUT_Path2TSD.get(pathOfInterest)
if isinstance(tsd, TimeSeriesDate):
return tsd
else:
for tsd in self.mTSDs:
assert isinstance(tsd, TimeSeriesDate)
if pathOfInterest in tsd.sourceUris():
return tsd

benjamin.jakimow@geo.hu-berlin.de
committed
return None
def tsd(self, date: np.datetime64, sensor)->TimeSeriesDate:
Returns the TimeSeriesDate identified by date and sensorID

Benjamin Jakimow
committed
:param date: numpy.datetime64
:param sensor: SensorInstrument | str with sensor id
:return:
"""
assert isinstance(date, np.datetime64)
if isinstance(sensor, str):
sensor = self.sensor(sensor)
if isinstance(sensor, SensorInstrument):
for tsd in self.mTSDs:
if tsd.date() == date and tsd.sensor() == sensor:
return tsd

Benjamin Jakimow
committed
else:
for tsd in self.mTSDs:
if tsd.date() == date:
return tsd
def insertTSD(self, tsd: TimeSeriesDate)->TimeSeriesDate:
Inserts a TimeSeriesDate
:param tsd: TimeSeriesDate
"""
#insert sorted by time & sensor
assert tsd not in self.mTSDs
assert tsd.sensor() in self.mSensors
tsd.sigRemoveMe.connect(lambda tsd=tsd: self.removeTSDs([tsd]))
tsd.rowsAboutToBeRemoved.connect(lambda p, first, last, tsd=tsd: self.beginRemoveRows(self.tsdToIdx(tsd), first, last))
tsd.rowsRemoved.connect(self.endRemoveRows)
tsd.rowsAboutToBeInserted.connect(lambda p, first, last, tsd=tsd: self.beginInsertRows(self.tsdToIdx(tsd), first, last))
tsd.rowsInserted.connect(self.endInsertRows)
tsd.sigSourcesAdded.connect(self.sigSourcesAdded)
tsd.sigSourcesRemoved.connect(self.sigSourcesRemoved)
row = bisect.bisect(self.mTSDs, tsd)
self.beginInsertRows(self.mRootIndex, row, row)
self.mTSDs.insert(row, tsd)
self.endInsertRows()
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
def showTSDs(self, tsds:list, b:bool=True):
tsds = [t for t in tsds if t in self]
minRow = None
maxRow = None
for t in tsds:
idx = self.tsdToIdx(t)
if minRow is None:
minRow = idx.row()
maxRow = idx.row()
else:
minRow = min(minRow, idx.row())
maxRow = min(maxRow, idx.row())
assert isinstance(t, TimeSeriesDate)
t.setVisibility(b)
if minRow:
ul = self.index(minRow, 0)
lr = self.index(maxRow, 0)
self.dataChanged.emit(ul, lr, [Qt.CheckStateRole])
self.sigVisibilityChanged.emit()
def hideTSDs(self, tsds):
self.showTSDs(tsds, False)
Removes a list of TimeSeriesDate
:param tsds: [list-of-TimeSeriesDate]
toRemove = set()
for t in tsds:
if isinstance(t, TimeSeriesDate):
toRemove.add(t)
if isinstance(t, TimeSeriesSource):
toRemove.add(t.timeSeriesDate())
for tsd in list(sorted(list(toRemove), reverse=True)):
assert isinstance(tsd, TimeSeriesDate)
tsd.sigSourcesRemoved.disconnect()
tsd.sigSourcesAdded.disconnect()
tsd.sigRemoveMe.disconnect()
row = self.mTSDs.index(tsd)
self.beginRemoveRows(self.mRootIndex, row, row)
self.mTSDs.remove(tsd)
tsd.mTimeSeries = None
removed.append(tsd)
self.endRemoveRows()

Benjamin Jakimow
committed
if len(removed) > 0:
pathsToRemove = [path for path, tsd in self.mLUT_Path2TSD.items() if tsd in removed]
for path in pathsToRemove:
self.mLUT_Path2TSD.pop(path)
self.checkSensorList()

Benjamin Jakimow
committed
self.sigTimeSeriesDatesRemoved.emit(removed)
def tsds(self, date:np.datetime64=None, sensor:SensorInstrument=None)->list:
Returns a list of TimeSeriesDate of the TimeSeries. By default all TimeSeriesDate will be returned.
:param date: numpy.datetime64 to return the TimeSeriesDate for
:param sensor: SensorInstrument of interest to return the [list-of-TimeSeriesDate] for.
:return: [list-of-TimeSeriesDate]
tsds = self.mTSDs[:]
if date:
tsds = [tsd for tsd in tsds if tsd.date() == date]
if sensor:

Benjamin Jakimow
committed
tsds = [tsd for tsd in tsds if tsd.sensor() == sensor]
"""
Removes all data sources from the TimeSeries (which will be empty after calling this routine).
"""
self.removeTSDs(self[:])
def addSensor(self, sensor:SensorInstrument):
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
if not sensor in self.mSensors:
self.mSensors.append(sensor)
self.sigSensorAdded.emit(sensor)
return sensor
else:
return None
def checkSensorList(self):
"""
Removes sensors without linked TSD / no data
"""
to_remove = []
for sensor in self.sensors():
tsds = [tsd for tsd in self.mTSDs if tsd.sensor() == sensor]
if len(tsds) == 0:
to_remove.append(sensor)
for sensor in to_remove:
self.removeSensor(sensor)
def removeSensor(self, sensor:SensorInstrument)->SensorInstrument:
"""
Removes a sensor and all linked images
:param sensor: SensorInstrument
:return: SensorInstrument or none, if sensor was not defined in the TimeSeries
"""
assert isinstance(sensor, SensorInstrument)
if sensor in self.mSensors:
tsds = [tsd for tsd in self.mTSDs if tsd.sensor() == sensor]
self.removeTSDs(tsds)
self.mSensors.remove(sensor)
self.sigSensorRemoved.emit(sensor)
return sensor
return None
def addSources(self, sources:list, nWorkers:int = 1, progressDialog:QProgressDialog=None, runAsync=True):
"""
Adds source images to the TimeSeries
:param sources: list of source images, e.g. a list of file paths
:param nWorkers: not used yet
:param progressDialog: QProgressDialog
:param runAsync: bool
"""
tm = QgsApplication.taskManager()
assert isinstance(tm, QgsTaskManager)
assert isinstance(nWorkers, int) and nWorkers >= 1
self.mLoadingProgressDialog = progressDialog
n = len(sources)
taskDescription = 'Load {} images'.format(n)
dump = pickle.dumps(sources)
if runAsync:
qgsTask = QgsTask.fromFunction(taskDescription, doLoadTimeSeriesSourcesTask, dump,
on_finished = self.onAddSourcesAsyncFinished)
else:
from .utils import TaskMock
qgsTask = TaskMock()
tid = id(qgsTask)
qgsTask.taskCompleted.connect(lambda *args, tid=tid: self.onRemoveTask(tid))
qgsTask.taskTerminated.connect(lambda *args, tid=tid: self.onRemoveTask(tid))
if isinstance(progressDialog, QProgressDialog):
qgsTask.progressChanged.connect(progressDialog.setValue)
progressDialog.setLabelText('Read images')
progressDialog.setRange(0, n)
progressDialog.setValue(0)
self.mTasks[tid] = qgsTask
tm.addTask(qgsTask)
else:
resultDump = doLoadTimeSeriesSourcesTask(qgsTask, dump)
self.onAddSourcesAsyncFinished(None, resultDump)
def onRemoveTask(self, key):
self.mTasks.pop(key)
def onAddSourcesAsyncFinished(self, *args):
hasProgressDialog = isinstance(self.mLoadingProgressDialog, QProgressDialog)
if error is None:
try:
addedDates = []
dump = args[1]
sources = pickle.loads(dump)
if hasProgressDialog:
self.mLoadingProgressDialog.setLabelText('Add Images')
self.mLoadingProgressDialog.setRange(0, len(sources))
self.mLoadingProgressDialog.setValue(0)
for i, source in enumerate(sources):
if isinstance(newTSD, TimeSeriesDate):
if hasProgressDialog:
self.mLoadingProgressDialog.setValue(i+1)
if len(addedDates) > 0:
self.sigTimeSeriesDatesAdded.emit(addedDates)
except Exception as ex:

Benjamin Jakimow
committed
import traceback
traceback.print_exc()
if isinstance(self.mLoadingProgressDialog, QProgressDialog):
self.mLoadingProgressDialog.hide()
self.mLoadingProgressDialog = None
def _addSource(self, source:TimeSeriesSource)->TimeSeriesDate:
"""
:param source:
:return: TimeSeriesDate (if new created)
"""
if isinstance(source, TimeSeriesSource):
tss = source
else:
tss = TimeSeriesSource.create(source)
assert isinstance(tss, TimeSeriesSource)
newTSD = None
tsdDate = self.date2date(tss.date())
tssDate = tss.date()
sid = tss.sid()
sensor = self.sensor(sid)
# if necessary, add a new sensor instance
if not isinstance(sensor, SensorInstrument):
sensor = self.addSensor(SensorInstrument(sid))
assert isinstance(sensor, SensorInstrument)
# if necessary, add a new TimeSeriesDate instance
if not isinstance(tsd, TimeSeriesDate):
tsd = self.insertTSD(TimeSeriesDate(self, tsdDate, sensor))
newTSD = tsd
# addedDates.append(tsd)
assert isinstance(tsd, TimeSeriesDate)
# add the source
tsd.addSource(tss)
return newTSD

Benjamin Jakimow
committed
def setDateTimePrecision(self, mode:DateTimePrecision):
"""
Sets the precision with which the parsed DateTime information will be handled.
:param mode: TimeSeriesViewer:DateTimePrecision
:return:
"""
self.mDateTimePrecision = mode
#do we like to update existing sources?
def date2date(self, date:np.datetime64)->np.datetime64:
"""

Benjamin Jakimow
committed
Converts a date of arbitrary precision into the date with precision according to the EOTSV settions.
:param date: numpy.datetime64
:return: numpy.datetime64
"""

Benjamin Jakimow
committed
assert isinstance(date, np.datetime64)
if self.mDateTimePrecision == DateTimePrecision.Original:
return date
else:
date = np.datetime64(date, self.mDateTimePrecision.value)
return date
def sources(self) -> list:
"""
Returns the input sources
:return: iterator over [list-of-TimeSeriesSources]
"""
for tsd in self:
for source in tsd:
yield source

Benjamin Jakimow
committed
def sourceUris(self)->list:
"""
Returns the uris of all sources
:return: [list-of-str]
"""
uris = []
for tsd in self:
assert isinstance(tsd, TimeSeriesDate)
def __len__(self):
def __iter__(self):
def __getitem__(self, slice):
def __delitem__(self, slice):
def __contains__(self, item):
def __repr__(self):
info = []
info.append('TimeSeries:')
l = len(self)
info.append(' Scenes: {}'.format(l))
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
def headerData(self, section, orientation, role):
assert isinstance(section, int)
if orientation == Qt.Horizontal and role == Qt.DisplayRole:
if len(self.mColumnNames) > section:
return self.mColumnNames[section]
else:
return ''
else:
return None
def parent(self, index: QModelIndex) -> QModelIndex:
"""
Returns the parent index of a QModelIndex `index`
:param index: QModelIndex
:return: QModelIndex
"""
if not index.isValid():
return QModelIndex()
node = index.internalPointer()
tsd = None
tss = None
if isinstance(node, TimeSeriesDate):
return self.mRootIndex
elif isinstance(node, TimeSeriesSource):
tss = node
return self.createIndex(self.mTSDs.index(tsd), 0, tsd)
def rowCount(self, index: QModelIndex=None) -> int:
"""
Return the row-count, i.e. number of child node for a TreeNode as index `index`.
:param index: QModelIndex
:return: int
"""
if index is None:
index = QModelIndex()
if not index.isValid():
return len(self)
node = index.internalPointer()
if isinstance(node, TimeSeriesDate):
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
return len(node)
if isinstance(node, TimeSeriesSource):
return 0
def columnNames(self) -> list:
"""
Returns the column names
:return: [list-of-string]
"""
return self.mColumnNames[:]
def columnCount(self, index:QModelIndex = None) -> int:
"""
Returns the number of columns
:param index: QModelIndex
:return:
"""
return len(self.mColumnNames)
def connectTreeView(self, treeView):
self.mTreeView = treeView
def index(self, row: int, column: int, parent: QModelIndex = None) -> QModelIndex:
"""
Returns the QModelIndex
:param row: int
:param column: int
:param parent: QModelIndex
:return: QModelIndex
"""
if parent is None:
parent = self.mRootIndex
else:
assert isinstance(parent, QModelIndex)
if row < 0 or row >= len(self):
return QModelIndex()
if column < 0 or column >= len(self.mColumnNames):
return QModelIndex()
if parent == self.mRootIndex:
# TSD node
if row < 0 or row >= len(self):
return QModelIndex()
return self.createIndex(row, column, self[row])
elif parent.parent() == self.mRootIndex:
# TSS node
tsd = self.tsdFromIdx(parent)
if row < 0 or row >= len(tsd):
return QModelIndex()
return self.createIndex(row, column, tsd[row])
return QModelIndex()
def tsdToIdx(self, tsd:TimeSeriesDate)->QModelIndex:
"""
Returns an QModelIndex pointing on a TimeSeriesDate of interest
:param tsd: TimeSeriesDate
:return: QModelIndex
"""
row = self.mTSDs.index(tsd)
return self.index(row, 0)
def tsdFromIdx(self, index: QModelIndex) -> TimeSeriesDate:
"""
Returns the TimeSeriesDate related to an QModelIndex `index`.
:param index: QModelIndex
:return: TreeNode
"""
if index.row() == -1 and index.column() == -1:
return None
elif not index.isValid():
return None
else:
node = index.internalPointer()
if isinstance(node, TimeSeriesDate):
return node
elif isinstance(node, TimeSeriesSource):
return None
def data(self, index, role):
"""
:param index: QModelIndex
:param role: Qt.ItemRole
:return: object
"""
assert isinstance(index, QModelIndex)
if not index.isValid():
return None
node = index.internalPointer()
tsd = None
tss = None
if isinstance(node, TimeSeriesSource):
tss = node
elif isinstance(node, TimeSeriesDate):
tsd = node
if role == Qt.UserRole:
return node
cName = self.mColumnNames[index.column()]
if role in [Qt.DisplayRole]:
if cName == self.cnDate:
if cName == self.cnImages:
return tss.uri()
if cName == self.cnNB:
return tss.nb
if cName == self.cnNL:
return tss.nl
if cName == self.cnNS:
return tss.ns
if cName == self.cnCRS:
return tss.crs().description()
if cName == self.cnSensor:
return tsd.sensor().name()
if role == Qt.DecorationRole and index.column() == 0:
return None
ext = tss.spatialExtent()
if isinstance(self.mCurrentSpatialExtent, SpatialExtent) and isinstance(ext, SpatialExtent):
ext = ext.toCrs(self.mCurrentSpatialExtent.crs())
b = isinstance(ext, SpatialExtent) and ext.intersects(self.mCurrentSpatialExtent)
if b:
return QIcon(r':/timeseriesviewer/icons/mapview.svg')
else:
return QIcon(r':/timeseriesviewer/icons/mapviewHidden.svg')
else:
print(ext)
return None
if role == Qt.BackgroundColorRole and tsd in self.mVisibleDate:
if isinstance(node, TimeSeriesDate):
if role in [Qt.DisplayRole]:
if cName == self.cnSensor:
return tsd.sensor().name()
if cName == self.cnImages:
return len(tsd)
if cName == self.cnDate:
return str(tsd.date())
return Qt.Checked if tsd.isVisible() else Qt.Unchecked
if role == Qt.BackgroundColorRole and tsd in self.mVisibleDate:
return None
def setData(self, index: QModelIndex, value: typing.Any, role: int):
if not index.isValid():
return False
result = False
bVisibilityChanged = False
node = index.internalPointer()
if isinstance(node, TimeSeriesDate):
if role == Qt.CheckStateRole and index.column() == 0:
node.setVisibility(value == Qt.Checked)
result = True
bVisibilityChanged = True
if result == True:
self.dataChanged.emit(index, index, [role])
if bVisibilityChanged:
self.sigVisibilityChanged.emit()
return result

Benjamin Jakimow
committed
def findDate(self, date)->TimeSeriesDate:
"""
Returns a TimeSeriesDate closest to that in date

Benjamin Jakimow
committed
:param date: numpy.datetime64 | str | TimeSeriesDate
:return: TimeSeriesDate
"""
if isinstance(date, str):
date = np.datetime64(date)
if isinstance(date, TimeSeriesDate):
date = date.date()
assert isinstance(date, np.datetime64)
if len(self) == 0:
return None
dtAbs = np.abs(date - np.asarray([tsd.date() for tsd in self.mTSDs]))
i = np.argmin(dtAbs)
return self.mTSDs[i]
def flags(self, index):
assert isinstance(index, QModelIndex)
if not index.isValid():
return Qt.NoItemFlags
#cName = self.mColumnNames.index(index.column())
flags = Qt.ItemIsEnabled | Qt.ItemIsSelectable
if isinstance(index.internalPointer(), TimeSeriesDate) and index.column() == 0:
flags = flags | Qt.ItemIsUserCheckable
return flags
def getSpatialPropertiesFromDataset(ds):
assert isinstance(ds, gdal.Dataset)
nb = ds.RasterCount
nl = ds.RasterYSize
ns = ds.RasterXSize
proj = ds.GetGeoTransform()
px_x = float(abs(proj[1]))
px_y = float(abs(proj[5]))
crs = QgsCoordinateReferenceSystem(ds.GetProjection())
return nb, nl, ns, crs, px_x, px_y

Benjamin Jakimow
committed
def extractWavelengthsFromGDALMetaData(ds:gdal.Dataset)->(list, str):

Benjamin Jakimow
committed
Reads the wavelength info from standard metadata strings
:param ds: gdal.Dataset

Benjamin Jakimow
committed
:return: (list, str)

Benjamin Jakimow
committed
regWLkey = re.compile('^(center )?wavelength[_ ]*$', re.I)
regWLUkey = re.compile('^wavelength[_ ]*units?$', re.I)
regNumeric = re.compile(r"([-+]?\d*\.\d+|[-+]?\d+)", re.I)

Benjamin Jakimow
committed
def findKey(d:dict, regex)->str:
for key in d.keys():
if regex.search(key):
return key

Benjamin Jakimow
committed
# 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()

Benjamin Jakimow
committed
keyWLU = findKey(md, regWLUkey)
keyWL = findKey(md, regWLkey)

Benjamin Jakimow
committed
if isinstance(keyWL, str) and isinstance(keyWLU, str):
wl.append(float(md[keyWL]))
wlu.append(LUT_WAVELENGTH_UNITS[md[keyWLU].lower()])

Benjamin Jakimow
committed
if len(wlu) == len(wl) and len(wl) == ds.RasterCount:
return wl, wlu[0]

Benjamin Jakimow
committed
# 2. try data set level
for domain in ds.GetMetadataDomainList():
md = ds.GetMetadata_Dict(domain)

Benjamin Jakimow
committed
keyWLU = findKey(md, regWLUkey)
keyWL = findKey(md, regWLkey)
if isinstance(keyWL, str) and isinstance(keyWLU, str):

Benjamin Jakimow
committed
wlu = LUT_WAVELENGTH_UNITS[md[keyWLU].lower()]
matches = regNumeric.findall(md[keyWL])
wl = [float(n) for n in matches]
if len(wl) == ds.RasterCount:
return wl, wlu

Benjamin Jakimow
committed
return None, None

Benjamin Jakimow
committed
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),

Benjamin Jakimow
committed
]
return wl, wlu
return None, None

Benjamin Jakimow
committed
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
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):
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