Fix generation of raster elevation profiles for exactly horizontal/

vertical lines Fixes #51196
qgis · Jan 31, 2023 · 9e4c39e · 9e4c39e
1 parent ca6ea75
commit 9e4c39e
Show file tree

Hide file tree

Showing 2 changed files with 75 additions and 10 deletions.
diff --git a/src/core/raster/qgsrasterlayerprofilegenerator.cpp b/src/core/raster/qgsrasterlayerprofilegenerator.cpp
@@ -178,16 +178,16 @@ bool QgsRasterLayerProfileGenerator::generateProfile( const QgsProfileGeneration
   int subRegionHeight = 0;
   int subRegionLeft = 0;
   int subRegionTop = 0;
-  const QgsRectangle rasterSubRegion = mRasterProvider->xSize() > 0 && mRasterProvider->ySize() > 0 ?
-                                       QgsRasterIterator::subRegion(
-                                         mRasterProvider->extent(),
-                                         mRasterProvider->xSize(),
-                                         mRasterProvider->ySize(),
-                                         transformedCurve->boundingBox(),
-                                         subRegionWidth,
-                                         subRegionHeight,
-                                         subRegionLeft,
-                                         subRegionTop ) : transformedCurve->boundingBox();
+  QgsRectangle rasterSubRegion = mRasterProvider->xSize() > 0 && mRasterProvider->ySize() > 0 ?
+                                 QgsRasterIterator::subRegion(
+                                   mRasterProvider->extent(),
+                                   mRasterProvider->xSize(),
+                                   mRasterProvider->ySize(),
+                                   transformedCurve->boundingBox(),
+                                   subRegionWidth,
+                                   subRegionHeight,
+                                   subRegionLeft,
+                                   subRegionTop ) : transformedCurve->boundingBox();
 
   const bool zeroXYSize = mRasterProvider->xSize() == 0 || mRasterProvider->ySize() == 0;
   if ( zeroXYSize )
@@ -196,6 +196,19 @@ bool QgsRasterLayerProfileGenerator::generateProfile( const QgsProfileGeneration
     const double conversionFactor = curveLengthInPixels / transformedCurve->length();
     subRegionWidth = rasterSubRegion.width() * conversionFactor;
     subRegionHeight = rasterSubRegion.height() * conversionFactor;
+
+    // ensure we fetch at least 1 pixel wide/high blocks, otherwise exactly vertical/horizontal profile lines will result in zero size blocks
+    // see https://github.com/qgis/QGIS/issues/51196
+    if ( subRegionWidth == 0 )
+    {
+      subRegionWidth = 1;
+      rasterSubRegion.setXMaximum( rasterSubRegion.xMinimum() + 1 / conversionFactor );
+    }
+    if ( subRegionHeight == 0 )
+    {
+      subRegionHeight = 1;
+      rasterSubRegion.setYMaximum( rasterSubRegion.yMinimum() + 1 / conversionFactor );
+    }
   }
 
   // iterate over the raster blocks, throwing away any which don't intersect the profile curve

diff --git a/tests/src/python/test_qgsrasterlayerprofilegenerator.py b/tests/src/python/test_qgsrasterlayerprofilegenerator.py
@@ -110,6 +110,58 @@ def testGenerationWithStepSize(self):
         self.assertEqual(r.zRange().lower(), 74)
         self.assertEqual(r.zRange().upper(), 154)
 
+    def testGenerationWithVerticalLine(self):
+        rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
+        self.assertTrue(rl.isValid())
+
+        curve = QgsLineString()
+        curve.fromWkt('LineString (321878.13400000002002344 130592.75222520538954996, 321878.13400000002002344 129982.02943174661777448)')
+        req = QgsProfileRequest(curve)
+        req.setStepDistance(10)
+
+        rl.elevationProperties().setEnabled(True)
+
+        req.setCrs(QgsCoordinateReferenceSystem('EPSG:27700'))
+        generator = rl.createProfileGenerator(req)
+        self.assertTrue(generator.generateProfile())
+
+        r = generator.takeResults()
+        results = r.distanceToHeightMap()
+        self.assertEqual(len(results), 63)
+        first_point = min(results.keys())
+        last_point = max(results.keys())
+        self.assertEqual(results[first_point], 120)
+        self.assertEqual(results[last_point], 86)
+
+        self.assertEqual(r.zRange().lower(), 74)
+        self.assertEqual(r.zRange().upper(), 120)
+
+    def testGenerationWithHorizontallLine(self):
+        rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
+        self.assertTrue(rl.isValid())
+
+        curve = QgsLineString()
+        curve.fromWkt('LineString (321471.82703730149660259 130317.67500000000291038, 322294.53625493601430207 130317.67500000000291038)')
+        req = QgsProfileRequest(curve)
+        req.setStepDistance(10)
+
+        rl.elevationProperties().setEnabled(True)
+
+        req.setCrs(QgsCoordinateReferenceSystem('EPSG:27700'))
+        generator = rl.createProfileGenerator(req)
+        self.assertTrue(generator.generateProfile())
+
+        r = generator.takeResults()
+        results = r.distanceToHeightMap()
+        self.assertEqual(len(results), 84)
+        first_point = min(results.keys())
+        last_point = max(results.keys())
+        self.assertEqual(results[first_point], 122)
+        self.assertEqual(results[last_point], 122)
+
+        self.assertEqual(r.zRange().lower(), 76)
+        self.assertEqual(r.zRange().upper(), 130)
+
     def testSnapping(self):
         rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
         self.assertTrue(rl.isValid())