diff --git a/timeseriesviewer/instruments.txt b/timeseriesviewer/instruments.txt index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..1f6a23f6dd986ab1f0f9c9656f945b8d58efc8a4 100644 --- a/timeseriesviewer/instruments.txt +++ b/timeseriesviewer/instruments.txt @@ -0,0 +1,3 @@ +#define new instruments by: +#<name>;<spat. resolution x>;<spat. resolution y>;<#bands>;[<comma separated center wavelength with unit>] +Landsat legacy bands; 30; 30;6; \ No newline at end of file diff --git a/timeseriesviewer/main.py b/timeseriesviewer/main.py index dabe33f2de3bd1de08163c5caffa0e537114f311..abe58628baa04472a1545bf8a5b3961f950015f0 100644 --- a/timeseriesviewer/main.py +++ b/timeseriesviewer/main.py @@ -812,11 +812,52 @@ def getBoundingBoxPolygon(points, srs=None): +def transformGeometry(geom, crsSrc, crsDst, trans=None): + if trans is None: + assert isinstance(crsSrc, QgsCoordinateReferenceSystem) + assert isinstance(crsDst, QgsCoordinateReferenceSystem) + return transformGeometry(geom, None, None, trans=QgsCoordinateTransform(crsSrc, crsDst)) + else: + assert isinstance(trans, QgsCoordinateTransform) + return trans.transform(geom) + + + + def getDS(ds): if type(ds) is not gdal.Dataset: ds = gdal.Open(ds) return ds +regAcqDate = re.compile(r'acquisition[ ]*(time|date|day)', re.I) +def getImageDate2(lyr): + assert isinstance(lyr, QgsRasterLayer) + mdLines = str(lyr.metadata()).splitlines() + date = None + #find date in metadata + for line in [l for l in mdLines if regAcqDate.search(l)]: + date = parseAcquisitionDate(line) + if date: + return date + #find date in filename + dn, fn = os.path.split(str(lyr.dataProvider().dataSourceUri())) + date = parseAcquisitionDate(fn) + if date: return date + + #find date in file directory path + date = parseAcquisitionDate(date) + + return date + +def getBandNames(lyr): + assert isinstance(lyr, QgsRasterLayer) + dp = lyr.dataProvider() + assert isinstance(dp, QgsRasterDataProvider) + if str(dp.name()) == 'gdal': + s = "" + else: + return lyr + def getImageDate(ds): if type(ds) is str: ds = gdal.Open(ds) @@ -837,8 +878,6 @@ def getImageDate(ds): raise Exception('Can not identify acquisition date of {}'.format(path)) - - class TimeSeriesDatum(object): """ Collects all data sets related to one sensor @@ -848,55 +887,17 @@ class TimeSeriesDatum(object): self.pathImg = pathImg self.pathMsk = None - self.lyrImg = QgsRasterLayer(pathImg) - - dsImg = gdal.Open(pathImg) - assert dsImg + self.lyrImg = QgsRasterLayer(pathImg, os.path.basename(pathImg), False) + self.crs = self.lyrImg.dataProvider().crs() + self.sensor = SensorInstrument.fromRasterLayer(self.lyrImg) - date = getImageDate(dsImg) - assert date is not None - self.date = date.astype(str) + self.date = getImageDate2(self.lyrImg) + assert self.date is not None, 'Unable to find acquisition date of {}'.format(pathImg) - self.ns = dsImg.RasterXSize - self.nl = dsImg.RasterYSize - self.nb = dsImg.RasterCount - - self.srs_wkt = dsImg.GetProjection() - self.gt = list(dsImg.GetGeoTransform()) - - refBand = dsImg.GetRasterBand(1) - self.etype = refBand.DataType - - self.nodata = refBand.GetNoDataValue() - - - - self.bandnames = list() - for b in range(self.nb): - name = dsImg.GetRasterBand(b+1).GetDescription() - if name is None or name == '': - name = 'Band {}'.format(b+1) - self.bandnames.append(name) - - self.wavelength = None - domains = dsImg.GetMetadataDomainList() - if domains: - for domain in domains: - md = dsImg.GetMetadata_Dict(domain) - - if 'wavelength' in md.keys(): - try: - wl = md['wavelength'] - wl = re.split('[;,{} ]+', wl) - wl = [w for w in wl if len(w) > 0] - wl = [float(w) for w in wl] - assert len(wl) == self.nb - self.wavelength = wl - break - except: - pass - - #self.sensor = SensorInstrument(self.nb, self.gt[1], self.gt[5], self.bandnames, self.wavelength) + self.ns = self.lyrImg.width() + self.nl = self.lyrImg.height() + self.nb = self.lyrImg.bandCount() + self.srs_wkt = str(self.crs.toWkt()) if pathMsk: @@ -908,31 +909,21 @@ class TimeSeriesDatum(object): def getDate(self): return np.datetime64(self.date) + def getSpatialReference(self): + return self.crs srs = osr.SpatialReference() srs.ImportFromWkt(self.srs_wkt) return srs - def getBoundingBox(self, srs=None): - ext = list() - - - for px in [0,self.ns]: - for py in [0, self.nl]: - ext.append(px2Coordinate(self.gt, px, py)) + def getBoundingBox(self, srs=None): - - if srs is not None: - assert type(srs) is osr.SpatialReference - my_srs = self.getSpatialReference() - if not my_srs.IsSame(srs): - #todo: consider srs differences - trans = osr.CoordinateTransformation(my_srs, srs) - ext = trans.TransformPoints(ext) - ext = [(e[0], e[1]) for e in ext] - - return ext + bbox = self.lyrImg.extent() + if srs: + assert isinstance(srs, QgsCoordinateReferenceSystem) + bbox = transformGeometry(bbox, self.crs, srs) + return bbox def setMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True): diff --git a/timeseriesviewer/sandbox.py b/timeseriesviewer/sandbox.py index c678dd63a25bc6b1536ae67c12a14f9f984907e0..e15b8a535c485c43721381a95b692b55f919f8e8 100644 --- a/timeseriesviewer/sandbox.py +++ b/timeseriesviewer/sandbox.py @@ -17,8 +17,55 @@ class HiddenCanvas(QgsMapCanvas): +class LayerLoaderR(QRunnable): + + layerReady = pyqtSignal(QgsRasterLayer) + finished = pyqtSignal(list) + def __init__(self): + super(LayerLoader, self).__init__() + + def loadLayers(self, paths): + lyrs = [] + for path in paths: + print('Load '+path) + lyr = QgsRasterLayer(path) + if lyr: + self.layerReady.emit(lyr) + lyrs.append(lyr) + + self.finished.emit(lyrs) + +class LayerLoader(QObject): + + layerReady = pyqtSignal(QgsRasterLayer) + finished = pyqtSignal(list) + def __init__(self): + super(LayerLoader, self).__init__() + + def loadLayers(self, paths): + lyrs = [] + for path in paths: + print('Load '+path) + lyr = QgsRasterLayer(path) + if lyr: + self.layerReady.emit(lyr) + lyrs.append(lyr) + + self.finished.emit(lyrs) +def getLoadedLayer(layer): + print('LAYER READY') + print(layer) + +def getLoadedLayers(layers): + for lyr in layers: + getLoadedLayer(lyr) + + +def loadingDone(): + print('DONE') + if __name__ == '__main__': import site, sys #add site-packages to sys.path as done by enmapboxplugin.py @@ -41,10 +88,33 @@ if __name__ == '__main__': qgsApp.setPrefixPath(PATH_QGS, True) qgsApp.initQgis() - #run tests - if True: test_gui() - if False: test_component() + from timeseriesviewer import file_search + paths = file_search(r'C:\Users\geo_beja\Repositories\QGIS_Plugins\SenseCarbonTSViewer\example\Images', + '*.bsq') + #run + import numpy as np + pool = QThreadPool() + pool.maxThreadCount(3) + + pool.start() + + t0 = np.datetime64('now') + #paths = [] + objThread = QThread() + loader = LayerLoader() + #loader.loadLayers(paths) + loader.moveToThread(objThread) + loader.finished.connect(objThread.quit) + #loader.layerReady.connect(getLoadedLayer) + objThread.started.connect(lambda:loader.loadLayers(paths)) + + objThread.finished.connect(loadingDone) + objThread.start() + + while not objThread.isFinished(): + + print('DT: {}'.format(np.datetime64('now') - t0)) #close QGIS qgsApp.exec_()