#2559

Committed 22 Aug 2025 01:01PM UTC coverage: 34.745% (-0.1%) from 34.841%

Build # #2559

Build Type

push

coveralls-python

Committed by

web-flow

Commit Message

Merge ce1721d94 into 7d44cf35b

Coverage Stats

23 of 104 new or added lines in 2 files covered. (22.12%)

4767 of 13720 relevant lines covered (34.74%)

0.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

21.57

/processing/water_depth_difference_algorithm.py

from pathlib import Path

import numpy as np
from osgeo import gdal
from qgis.core import (
    QgsProcessingAlgorithm,
    QgsProcessingContext,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterFileDestination,
    QgsProcessingUtils,
)

gdal.UseExceptions()

STYLE_DIR = Path(__file__).parent / "styles"


def get_authority_code(raster: gdal.Dataset) -> str:
    """Return authority code (e.g. EPSG:28992) for given raster"""
    srs = raster.GetProjection()
    key = "GEOGCS" if srs.IsGeographic() else "PROJCS"
    return srs.GetAuthorityCode(key)


def get_extent(raster):
    gt = raster.GetGeoTransform()
    ulx = gt[0]
    uly = gt[3]
    lrx = gt[0] + (raster.RasterXSize * gt[1])
    lry = gt[3] + (raster.RasterYSize * gt[5])  # '+' because gt[5] is negative
    return [ulx, uly, lrx, lry]


def get_shared_extent(extent1, extent2):
    ulx = max(extent1[0], extent2[0])
    uly = min(extent1[1], extent2[1])
    lrx = min(extent1[2], extent2[2])
    lry = max(extent1[3], extent2[3])

    return [ulx, uly, lrx, lry]


def water_depth_diff(raster1_fn: Path | str, raster2_fn: Path | str, output: Path | str):
    """
    Calculate difference between two overlapping water depth tiffs. Result is raster2 - raster1.
    Regards nodata as 0 water depth
    Nodatavalue of output is always 0
    If input rasters have different extents, difference raster will be calculated for the overlapping part of the input
    rasters
    """
    output = Path(output)

    raster1 = gdal.Open(raster1_fn)
    raster2 = gdal.Open(raster2_fn)

    # raise error if pixel sizes differ
    if not raster1.GetGeoTransform()[1] == raster2.GetGeoTransform()[1] and \
            raster1.GetGeoTransform()[5] == raster2.GetGeoTransform()[5]:
        raise ValueError("Input rasters have different pixel sizes")

    # raise error if pixel skews differ
    if not raster1.GetGeoTransform()[2] == raster2.GetGeoTransform()[2] and \
            raster1.GetGeoTransform()[4] == raster2.GetGeoTransform()[4]:
        raise ValueError("Input rasters have different pixel skew")

    # raise error if projections differ
    # compare on EPSG code to prevent errors from insignificant differences between the projections
    authority_code_raster1 = get_authority_code(raster1)
    authority_code_raster2 = get_authority_code(raster2)
    if not authority_code_raster1 == authority_code_raster2:
        raise ValueError(f"Input rasters have different CRS ({authority_code_raster1} vs {authority_code_raster2}")

    # Clip rasters if they do not have the same extent
    raster1extent = get_extent(raster1)
    raster2extent = get_extent(raster2)

    if raster1extent == raster2extent:
        raster1_array = raster1.ReadAsArray()
        raster2_array = raster2.ReadAsArray()
        gt = raster1.GetGeoTransform()

    else:
        shared_extent = get_shared_extent(raster1extent, raster2extent)

        raster1_shared_extent = gdal.Translate('', raster1, format='MEM', projWin=shared_extent)
        raster1_array = raster1_shared_extent.ReadAsArray()
        raster1_shared_extent = None

        raster2_shared_extent = gdal.Translate('', raster2, format='MEM', projWin=shared_extent)
        raster2_array = raster2_shared_extent.ReadAsArray()
        raster2_shared_extent = None

        gt = list(raster1.GetGeoTransform())  # (upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)
        gt[0] = shared_extent[0]
        gt[3] = shared_extent[1]
        gt = tuple(gt)

    # Replace all nodata pixels by 0
    ndv1 = raster1.GetRasterBand(1).GetNoDataValue()
    raster1_array[np.where(raster1_array == ndv1)] = 0
    ndv2 = raster2.GetRasterBand(1).GetNoDataValue()
    raster2_array[np.where(raster2_array == ndv2)] = 0

    # Calculate difference
    result = np.subtract(raster2_array, raster1_array)

    # Write to raster with 0 as nodatavalue
    height = result.shape[0]
    width = result.shape[1]
    wkt = raster1.GetProjection()

    if output.exists():
        output.unlink()

    dst_drv = gdal.GetDriverByName('GTiff')
    dst_ds = dst_drv.Create(
        output,
        xsize=width,
        ysize=height,
        bands=1,
        eType=gdal.GDT_Float32,
        options=["COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9", "TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512"],
    )
    dst_ds.SetGeoTransform(gt)
    dst_ds.SetProjection(wkt)
    dst_ds.GetRasterBand(1).SetNoDataValue(0)
    dst_ds.GetRasterBand(1).WriteArray(result)
    dst_ds = None


class WaterDepthDiffAlgorithm(QgsProcessingAlgorithm):
    INPUT_RASTER1 = "INPUT_RASTER1"
    INPUT_RASTER2 = "INPUT_RASTER2"
    OUTPUT_RASTER = "OUTPUT_RASTER"

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT_RASTER1, "First Raster"
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT_RASTER2, "Second Raster"
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.OUTPUT_RASTER,
                "Output raster",
                fileFilter="GeoTIFF (*.tif)"
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        raster1 = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER1, context)
        raster2 = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER2, context)
        output = self.parameterAsFileOutput(parameters, self.OUTPUT_RASTER, context)

        water_depth_diff(raster1.source(), raster2.source(), output)
        layer = QgsProcessingUtils.mapLayerFromString(output, context)
        self.output_layer_id = layer.id()
        context.temporaryLayerStore().addMapLayer(layer)
        layer_details = QgsProcessingContext.LayerDetails(
            "Water depth difference", context.project(), "Water depth difference"
        )
        context.addLayerToLoadOnCompletion(layer.id(), layer_details)

        return {self.OUTPUT_RASTER: output}

    def postProcessAlgorithm(self, context, feedback):
        output_layer = context.getMapLayer(self.output_layer_id)
        output_layer.loadNamedStyle(STYLE_DIR / "water_depth_difference.qml")
        context.project().addMapLayer(output_layer)

    def name(self):
        return "water_depth_diff"

    def displayName(self):
        return "Water depth difference"

    def group(self):
        return self.tr("Post-process results")

    def groupId(self):
        return "postprocessing"

    def shortHelpString(self):
        return "Calculate difference between two overlapping water depth rasters. Result is raster 2 - raster 1."

1	from pathlib import Path	1✔
2
3	import numpy as np	1✔
4	from osgeo import gdal	1✔
5	from qgis.core import (	1✔
6	QgsProcessingAlgorithm,
7	QgsProcessingContext,
8	QgsProcessingParameterRasterLayer,
9	QgsProcessingParameterFileDestination,
10	QgsProcessingUtils,
11	)
12
13	gdal.UseExceptions()	1✔
14
15	STYLE_DIR = Path(__file__).parent / "styles"	1✔
16
17
18	def get_authority_code(raster: gdal.Dataset) -> str:	1✔
19	"""Return authority code (e.g. EPSG:28992) for given raster"""
NEW 20	srs = raster.GetProjection()	×
NEW 21	key = "GEOGCS" if srs.IsGeographic() else "PROJCS"	×
NEW 22	return srs.GetAuthorityCode(key)	×
23
24
25	def get_extent(raster):	1✔
NEW 26	gt = raster.GetGeoTransform()	×
NEW 27	ulx = gt[0]	×
NEW 28	uly = gt[3]	×
NEW 29	lrx = gt[0] + (raster.RasterXSize * gt[1])	×
NEW 30	lry = gt[3] + (raster.RasterYSize * gt[5]) # '+' because gt[5] is negative	×
NEW 31	return [ulx, uly, lrx, lry]	×
32
33
34	def get_shared_extent(extent1, extent2):	1✔
NEW 35	ulx = max(extent1[0], extent2[0])	×
NEW 36	uly = min(extent1[1], extent2[1])	×
NEW 37	lrx = min(extent1[2], extent2[2])	×
NEW 38	lry = max(extent1[3], extent2[3])	×
39
NEW 40	return [ulx, uly, lrx, lry]	×
41
42
43	def water_depth_diff(raster1_fn: Path \| str, raster2_fn: Path \| str, output: Path \| str):	1✔
44	"""
45	Calculate difference between two overlapping water depth tiffs. Result is raster2 - raster1.
46	Regards nodata as 0 water depth
47	Nodatavalue of output is always 0
48	If input rasters have different extents, difference raster will be calculated for the overlapping part of the input
49	rasters
50	"""
NEW 51	output = Path(output)	×
52
NEW 53	raster1 = gdal.Open(raster1_fn)	×
NEW 54	raster2 = gdal.Open(raster2_fn)	×
55
56	# raise error if pixel sizes differ
NEW 57	if not raster1.GetGeoTransform()[1] == raster2.GetGeoTransform()[1] and \	×
58	raster1.GetGeoTransform()[5] == raster2.GetGeoTransform()[5]:
NEW 59	raise ValueError("Input rasters have different pixel sizes")	×
60
61	# raise error if pixel skews differ
NEW 62	if not raster1.GetGeoTransform()[2] == raster2.GetGeoTransform()[2] and \	×
63	raster1.GetGeoTransform()[4] == raster2.GetGeoTransform()[4]:
NEW 64	raise ValueError("Input rasters have different pixel skew")	×
65
66	# raise error if projections differ
67	# compare on EPSG code to prevent errors from insignificant differences between the projections
NEW 68	authority_code_raster1 = get_authority_code(raster1)	×
NEW 69	authority_code_raster2 = get_authority_code(raster2)	×
NEW 70	if not authority_code_raster1 == authority_code_raster2:	×
NEW 71	raise ValueError(f"Input rasters have different CRS ({authority_code_raster1} vs {authority_code_raster2}")	×
72
73	# Clip rasters if they do not have the same extent
NEW 74	raster1extent = get_extent(raster1)	×
NEW 75	raster2extent = get_extent(raster2)	×
76
NEW 77	if raster1extent == raster2extent:	×
NEW 78	raster1_array = raster1.ReadAsArray()	×
NEW 79	raster2_array = raster2.ReadAsArray()	×
NEW 80	gt = raster1.GetGeoTransform()	×
81
82	else:
NEW 83	shared_extent = get_shared_extent(raster1extent, raster2extent)	×
84
NEW 85	raster1_shared_extent = gdal.Translate('', raster1, format='MEM', projWin=shared_extent)	×
NEW 86	raster1_array = raster1_shared_extent.ReadAsArray()	×
NEW 87	raster1_shared_extent = None	×
88
NEW 89	raster2_shared_extent = gdal.Translate('', raster2, format='MEM', projWin=shared_extent)	×
NEW 90	raster2_array = raster2_shared_extent.ReadAsArray()	×
NEW 91	raster2_shared_extent = None	×
92
NEW 93	gt = list(raster1.GetGeoTransform()) # (upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)	×
NEW 94	gt[0] = shared_extent[0]	×
NEW 95	gt[3] = shared_extent[1]	×
NEW 96	gt = tuple(gt)	×
97
98	# Replace all nodata pixels by 0
NEW 99	ndv1 = raster1.GetRasterBand(1).GetNoDataValue()	×
NEW 100	raster1_array[np.where(raster1_array == ndv1)] = 0	×
NEW 101	ndv2 = raster2.GetRasterBand(1).GetNoDataValue()	×
NEW 102	raster2_array[np.where(raster2_array == ndv2)] = 0	×
103
104	# Calculate difference
NEW 105	result = np.subtract(raster2_array, raster1_array)	×
106
107	# Write to raster with 0 as nodatavalue
NEW 108	height = result.shape[0]	×
NEW 109	width = result.shape[1]	×
NEW 110	wkt = raster1.GetProjection()	×
111
NEW 112	if output.exists():	×
NEW 113	output.unlink()	×
114
NEW 115	dst_drv = gdal.GetDriverByName('GTiff')	×
NEW 116	dst_ds = dst_drv.Create(	×
117	output,
118	xsize=width,
119	ysize=height,
120	bands=1,
121	eType=gdal.GDT_Float32,
122	options=["COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9", "TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512"],
123	)
NEW 124	dst_ds.SetGeoTransform(gt)	×
NEW 125	dst_ds.SetProjection(wkt)	×
NEW 126	dst_ds.GetRasterBand(1).SetNoDataValue(0)	×
NEW 127	dst_ds.GetRasterBand(1).WriteArray(result)	×
NEW 128	dst_ds = None	×
129
130
131	class WaterDepthDiffAlgorithm(QgsProcessingAlgorithm):	1✔
132	INPUT_RASTER1 = "INPUT_RASTER1"	1✔
133	INPUT_RASTER2 = "INPUT_RASTER2"	1✔
134	OUTPUT_RASTER = "OUTPUT_RASTER"	1✔
135
136	def initAlgorithm(self, config=None):	1✔
NEW 137	self.addParameter(	×
138	QgsProcessingParameterRasterLayer(
139	self.INPUT_RASTER1, "First Raster"
140	)
141	)
142
NEW 143	self.addParameter(	×
144	QgsProcessingParameterRasterLayer(
145	self.INPUT_RASTER2, "Second Raster"
146	)
147	)
148
NEW 149	self.addParameter(	×
150	QgsProcessingParameterFileDestination(
151	self.OUTPUT_RASTER,
152	"Output raster",
153	fileFilter="GeoTIFF (*.tif)"
154	)
155	)
156
157	def processAlgorithm(self, parameters, context, feedback):	1✔
NEW 158	raster1 = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER1, context)	×
NEW 159	raster2 = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER2, context)	×
NEW 160	output = self.parameterAsFileOutput(parameters, self.OUTPUT_RASTER, context)	×
161
NEW 162	water_depth_diff(raster1.source(), raster2.source(), output)	×
NEW 163	layer = QgsProcessingUtils.mapLayerFromString(output, context)	×
NEW 164	self.output_layer_id = layer.id()	×
NEW 165	context.temporaryLayerStore().addMapLayer(layer)	×
NEW 166	layer_details = QgsProcessingContext.LayerDetails(	×
167	"Water depth difference", context.project(), "Water depth difference"
168	)
NEW 169	context.addLayerToLoadOnCompletion(layer.id(), layer_details)	×
170
NEW 171	return {self.OUTPUT_RASTER: output}	×
172
173	def postProcessAlgorithm(self, context, feedback):	1✔
NEW 174	output_layer = context.getMapLayer(self.output_layer_id)	×
NEW 175	output_layer.loadNamedStyle(STYLE_DIR / "water_depth_difference.qml")	×
NEW 176	context.project().addMapLayer(output_layer)	×
177
178	def name(self):	1✔
NEW 179	return "water_depth_diff"	×
180
181	def displayName(self):	1✔
NEW 182	return "Water depth difference"	×
183
184	def group(self):	1✔
NEW 185	return self.tr("Post-process results")	×
186
187	def groupId(self):	1✔
NEW 188	return "postprocessing"	×
189
190	def shortHelpString(self):	1✔
NEW 191	return "Calculate difference between two overlapping water depth rasters. Result is raster 2 - raster 1."	×

nens / ThreeDiToolbox / #2559

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous