search

EnMAP-Box / enmap-box

EnMAP-Box source code repository. See https://enmap-box.readthedocs.io for documentation
GNU General Public License v3.0
39 stars 16 forks source link

[API] support raster data reading with on-the-fly reprojection #1005

Closed janzandr closed 2 weeks ago

janzandr commented 2 weeks ago

It is proposed to support reading raster data in blocks that don't match the native raster CRS. We already have algorithms that allow the processing of raster layers with different CRS, but this involves explicite creation of a warped raster (usually a VRT). Proposed solution will use the QgsRasterPipe class, e.g.

from enmapboxprocessing.utils import Utils
from enmapboxtestdata import enmap_berlin
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem, QgsRasterPipe, QgsRasterProjector, QgsRasterBlock, QgsRectangle

# setup raster pipe for a layer with UTM projection
layer = QgsRasterLayer(enmap_berlin)
pipe = QgsRasterPipe()
pipe.set(layer.dataProvider().clone())
sourceCrs = QgsCoordinateReferenceSystem(layer.crs())
destinationCrs = QgsCoordinateReferenceSystem('EPSG:4326')
projector = QgsRasterProjector()
projector.setCrs(sourceCrs, destinationCrs)
pipe.insert(1, projector)

# read data in WGS84 projection
extent4326 = QgsRectangle(13.24539916899924741, 52.41274014201967901, 13.34677689752254182, 52.52184188795427389)
block: QgsRasterBlock = projector.block(1, extent4326, 10, 10)
array = Utils().qgsRasterBlockToNumpyArray(block)
janzandr commented 2 weeks ago

Let's just add a crs argument to the RasterReader constructor.

reader = RasterReader(enmap, crs=QgsCoordinateReferenceSystem('EPSG:4326'))
extent = QgsRectangle(13.24539916899924741, 52.41274014201967901, 13.34677689752254182, 52.52184188795427389)
array = reader.arrayFromBoundingBoxAndSize(extent, 220, 400, [1])
janzandr commented 2 weeks ago

The raster pipe can also be changed afterwards via RasterReader.setRasterPipeCrs(crs).