From 04913913b2e852917ad0760a03a4de2359cd883c Mon Sep 17 00:00:00 2001 From: Andreas Janz Date: Mon, 28 Oct 2024 14:52:33 +0100 Subject: [PATCH 1/3] resolved #1005 --- enmapboxprocessing/rasterreader.py | 30 ++++++++++++++++--- .../enmapboxprocessing/test_rasterreader.py | 12 +++++++- 2 files changed, 37 insertions(+), 5 deletions(-) diff --git a/enmapboxprocessing/rasterreader.py b/enmapboxprocessing/rasterreader.py index f396a4b5..9e045d1a 100644 --- a/enmapboxprocessing/rasterreader.py +++ b/enmapboxprocessing/rasterreader.py @@ -18,7 +18,7 @@ from qgis.core import (QgsRasterLayer, QgsRasterDataProvider, QgsCoordinateReferenceSystem, QgsRectangle, QgsRasterRange, QgsPoint, QgsRasterBlockFeedback, QgsRasterBlock, QgsPointXY, QgsProcessingFeedback, QgsRasterBandStats, Qgis, QgsGeometry, QgsVectorLayer, QgsWkbTypes, - QgsFeature) + QgsFeature, QgsRasterPipe, QgsRasterProjector) @typechecked @@ -26,8 +26,11 @@ class RasterReader(object): Nanometers = 'Nanometers' Micrometers = 'Micrometers' - def __init__(self, source: RasterSource, openWithGdal: bool = True): - + def __init__( + self, source: RasterSource, openWithGdal: bool = True, + crs: QgsCoordinateReferenceSystem = None + ): + self.maskReader: Optional[RasterReader] = None if isinstance(source, QgsRasterLayer): self.layer = source self.provider: QgsRasterDataProvider = self.layer.dataProvider() @@ -53,6 +56,8 @@ def __init__(self, source: RasterSource, openWithGdal: bool = True): self.gdalDataset = gdalDataset + self.setRasterPipeCrs(crs) + # prepare STAC metadata if exists(self.layer.source() + '.stac.json'): self.stacMetadata = Utils().jsonLoad(self.layer.source() + '.stac.json') @@ -71,6 +76,9 @@ def __init__(self, source: RasterSource, openWithGdal: bool = True): else: self.terraMetadata = None + def setExternalMask(self, layer: QgsRasterLayer): + self.maskReader = RasterReader(layer) + def bandCount(self) -> int: """Return iterator over all band numbers.""" return self.provider.bandCount() @@ -223,6 +231,19 @@ def walkGrid( continue # empty blocks may occure, but can just skip over yield RasterBlockInfo(blockExtent, xOffset, yOffset, width, height) + def setRasterPipeCrs(self, crs: QgsCoordinateReferenceSystem = None): + if crs is None: + projector = self.provider + elif isinstance(crs, QgsCoordinateReferenceSystem): + pipe = QgsRasterPipe() + pipe.set(self.provider.clone()) + projector = QgsRasterProjector() + projector.setCrs(self.provider.crs(), crs) + pipe.insert(1, projector) + else: + raise ValueError() + self.projector = projector + def arrayFromBlock( self, block: RasterBlockInfo, bandList: List[int] = None, overlap: int = None, feedback: QgsRasterBlockFeedback = None @@ -253,7 +274,7 @@ def arrayFromBoundingBoxAndSize( arrays = list() for bandNo in bandList: assert 0 < bandNo <= self.bandCount(), f'bandNo is {bandNo}' - block: QgsRasterBlock = self.provider.block(bandNo, boundingBox, width, height, feedback) + block: QgsRasterBlock = self.projector.block(bandNo, boundingBox, width, height, feedback) array = Utils.qgsRasterBlockToNumpyArray(block=block) arrays.append(array) return arrays @@ -304,6 +325,7 @@ def maskArray( maskNoDataValue=True ) -> Array3d: """Return mask for given data. No data values evaluate to False, all other to True.""" + if bandList is None: bandList = range(1, self.provider.bandCount() + 1) assert len(bandList) == len(array) diff --git a/tests/enmap-box/enmapboxprocessing/test_rasterreader.py b/tests/enmap-box/enmapboxprocessing/test_rasterreader.py index b57536d5..c28f669c 100644 --- a/tests/enmap-box/enmapboxprocessing/test_rasterreader.py +++ b/tests/enmap-box/enmapboxprocessing/test_rasterreader.py @@ -8,7 +8,7 @@ from enmapboxtestdata import enmap, r_terra_timeseries_days, r_terra_timeseries_seconds, netCDF_timeseries_days from enmapboxtestdata import fraction_polygon_l3 from qgis.PyQt.QtCore import QDateTime, QSizeF, QPoint -from qgis.core import QgsRasterRange, QgsRasterLayer, Qgis, QgsRectangle +from qgis.core import QgsRasterRange, QgsRasterLayer, Qgis, QgsRectangle, QgsCoordinateReferenceSystem class TestRasterReader(TestCase): @@ -699,3 +699,13 @@ def test_nonGdalSource(self): self.assertFalse(reader.isSpectralRasterLayer()) self.assertIsNone(reader.wavelength(1)) self.assertIsNone(reader.wavelengthUnits(1)) + + def test_arrayFromBoundingBoxAndSize_withRasterPipe(self): + reader = RasterReader(enmap, True, QgsCoordinateReferenceSystem('EPSG:4326')) + extent4326 = QgsRectangle( + 13.24539916899924741, 52.41274014201967901, 13.34677689752254182, 52.52184188795427389 + ) + array = reader.arrayFromBoundingBoxAndSize(extent4326, 10, 10, [1]) + print(array) + + # self.assertTrue(np.all(np.equal(reader.noDataValue(1), lead))) From 7c285578621b58df8a335deec72180b0244ad03a Mon Sep 17 00:00:00 2001 From: Andreas Janz Date: Mon, 28 Oct 2024 15:30:20 +0100 Subject: [PATCH 2/3] fixed url --- enmapbox/coreapps/_classic/hubflow/html.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/enmapbox/coreapps/_classic/hubflow/html.py b/enmapbox/coreapps/_classic/hubflow/html.py index b578d9f5..c2dd0d4b 100644 --- a/enmapbox/coreapps/_classic/hubflow/html.py +++ b/enmapbox/coreapps/_classic/hubflow/html.py @@ -16,7 +16,7 @@ #--- LICENSE ------------------------------------------------------------------ -# Copyright Philippe Lagadec - see http://www.decalage.info/contact for contact info +# Copyright Philippe Lagadec # # This module provides a few classes to easily generate HTML tables and lists. # From 2ee4436687726cdaa27382b3a0458b7f51864e4f Mon Sep 17 00:00:00 2001 From: Andreas Janz Date: Mon, 28 Oct 2024 15:30:24 +0100 Subject: [PATCH 3/3] Update rasterreader.py --- enmapboxprocessing/rasterreader.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/enmapboxprocessing/rasterreader.py b/enmapboxprocessing/rasterreader.py index 9e045d1a..3aeb247f 100644 --- a/enmapboxprocessing/rasterreader.py +++ b/enmapboxprocessing/rasterreader.py @@ -234,6 +234,7 @@ def walkGrid( def setRasterPipeCrs(self, crs: QgsCoordinateReferenceSystem = None): if crs is None: projector = self.provider + pipe = None elif isinstance(crs, QgsCoordinateReferenceSystem): pipe = QgsRasterPipe() pipe.set(self.provider.clone()) @@ -242,6 +243,7 @@ def setRasterPipeCrs(self, crs: QgsCoordinateReferenceSystem = None): pipe.insert(1, projector) else: raise ValueError() + self.pipe = pipe self.projector = projector def arrayFromBlock(