Verification.

Author: Andrew Bell

alg/viewshed/viewshed.cpp

@@ -373,6 +373,13 @@ bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)
 bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
                    void *pProgressArg)
 {
+	DatasetPtr sdDataset;
+	if (oOpts.sdFilename.size())
+	{
+		GDALDatasetH ds = GDALOpen(oOpts.sdFilename.c_str(), GA_ReadOnly);
+		sdDataset.reset(GDALDataset::FromHandle(ds));
+		pSdBand = sdDataset->GetRasterBand(1);
+	}
     pSrcBand = static_cast<GDALRasterBand *>(band);
 
     GDALGeoTransform fwdTransform, invTransform;
@@ -443,10 +450,20 @@ bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
 
     // Execute the viewshed algorithm.
     GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
-    ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
+    if (pSdBand)
+	{
+    	ViewshedExecutor executor(*pSrcBand, *pSdBand, *pDstBand, nX, nY,
+								  oOutExtent, oCurExtent, oOpts, oProgress,
+	                              /* emitWarningIfNoData = */ true);
+    	executor.run();
+	}
+	else
+	{
+    	ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
                               oCurExtent, oOpts, oProgress,
                               /* emitWarningIfNoData = */ true);
-    executor.run();
+    	executor.run();
+	}
     oProgress.emit(1);
     return static_cast<bool>(poDstDS);
 }

alg/viewshed/viewshed.h

@@ -69,6 +69,7 @@ class Viewshed
     Window oCurExtent{};
     DatasetPtr poDstDS{};
     GDALRasterBand *pSrcBand = nullptr;
+    GDALRasterBand *pSdBand = nullptr;
 
     DatasetPtr execute(int nX, int nY, const std::string &outFilename);
     void setOutput(double &dfResult, double &dfCellVal, double dfZ);

alg/viewshed/viewshed_executor.cpp

@@ -11,6 +11,9 @@
  * SPDX-License-Identifier: MIT
  ****************************************************************************/
 
+//ABELL
+#include <iostream>
+
 #include <algorithm>
 #include <atomic>
 #include <cassert>
@@ -295,7 +298,6 @@ bool ViewshedExecutor::readLine(int nLine, Lines &lines)
 
     if (sdMode())
     {
-        lines.input = lines.cur;
         double nodata = m_sdBand.GetNoDataValue();
         int sdStatus = m_sdBand.RasterIO(
             GF_Read, oOutExtent.xStart, nLine, oOutExtent.xSize(), 1,
@@ -497,6 +499,9 @@ bool ViewshedExecutor::processFirstLine(Lines &lines)
         m_dfZObserver += lines.cur[m_nX];
 
     LineLimits ll = adjustHeight(nYOffset, lines);
+	if (sdMode())
+		lines.input = lines.cur;
+
     if (oCurExtent.containsX(m_nX))
     {
         if (ll.leftMin != ll.rightMin)
@@ -513,9 +518,13 @@ bool ViewshedExecutor::processFirstLine(Lines &lines)
         else
         {
             CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
+/**
             pQueue->SubmitJob([&]() { processFirstLineLeft(ll, lines); });
             pQueue->SubmitJob([&]() { processFirstLineRight(ll, lines); });
             pQueue->WaitCompletion();
+**/
+processFirstLineLeft(ll, lines);
+processFirstLineRight(ll, lines);
         }
         sdCalc = false;
     };
@@ -871,6 +880,11 @@ void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
         else
             dfZ = oZcalc(nXOffset, nYOffset, lines.cur[iPixel + 1],
                          lines.prev[iPixel], lines.prev[iPixel + 1]);
+if (nLine == 2424 && iPixel == 530)
+{
+	std::cerr << "Val/SD val/dfZ = " << lines.cur[iPixel] << "/" <<
+		lines.sd[iPixel] << "/" << dfZ << "!\n";
+}
         setOutput(lines, iPixel, dfZ);
     }
 
@@ -961,6 +975,8 @@ bool ViewshedExecutor::processLine(int nLine, Lines &lines)
 
     // Adjust height of the read line.
     LineLimits ll = adjustHeight(nYOffset, lines);
+	if (sdMode())
+		lines.input = lines.cur;
 
     // Handle the initial position on the line.
     if (oCurExtent.containsX(m_nX))
@@ -991,6 +1007,8 @@ bool ViewshedExecutor::processLine(int nLine, Lines &lines)
     if (sdMode())
     {
         lines.cur = std::move(lines.input);
+if (nLine == 2424)
+	std::cerr << "Cur[530] = " << lines.cur[530] << "!\n";	
         process(true);
     }
 

alg/viewshed/viewshed_types.h

@@ -79,6 +79,7 @@ struct Options
     double curveCoeff{.85714};  //!< coefficient for atmospheric refraction
     OutputMode outputMode{OutputMode::Normal};  //!< Output information.
         //!< Normal, Height from DEM or Height from ground
+	std::string sdFilename{};           //!< sd raster filename
     std::string outputFormat{};         //!< output raster format
     std::string outputFilename{};       //!< output raster filename
     CPLStringList creationOpts{};       //!< options for output raster creation

apps/gdalalg_raster_viewshed.cpp

@@ -45,6 +45,9 @@ GDALRasterViewshedAlgorithm::GDALRasterViewshedAlgorithm(bool standaloneStep)
         .SetRepeatedArgAllowed(false);
     AddArg("height", 'z', _("Observer height"), &m_opts.observer.z);
 
+    AddArg("sd-filename", 0, _("Filename of standard-deviation raster"),
+           &m_opts.sdFilename);
+
     AddArg("target-height", 0,
            _("Height of the target above the DEM surface in the height unit of "
              "the DEM."),