Merge pull request #1006 from EnMAP-Box/1005-api-support-raster-data-…

…reading-with-on-the-fly-reprojection resolved #1005
EnMAP-Box · Oct 28, 2024 · 338a686 · 338a686
2 parents eb044ae + 2ee4436
commit 338a686
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 6 deletions.
diff --git 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.
 #

diff --git a/enmapboxprocessing/rasterreader.py b/enmapboxprocessing/rasterreader.py
@@ -18,16 +18,19 @@
 from qgis.core import (QgsRasterLayer, QgsRasterDataProvider, QgsCoordinateReferenceSystem, QgsRectangle,
                        QgsRasterRange, QgsPoint, QgsRasterBlockFeedback, QgsRasterBlock, QgsPointXY,
                        QgsProcessingFeedback, QgsRasterBandStats, Qgis, QgsGeometry, QgsVectorLayer, QgsWkbTypes,
-                       QgsFeature)
+                       QgsFeature, QgsRasterPipe, QgsRasterProjector)
 
 
 @typechecked
 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,21 @@ 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
+            pipe = None
+        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.pipe = pipe
+        self.projector = projector
+
     def arrayFromBlock(
             self, block: RasterBlockInfo, bandList: List[int] = None, overlap: int = None,
             feedback: QgsRasterBlockFeedback = None
@@ -253,7 +276,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 +327,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
@@ -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)))