LCOV - code coverage report
Current view: top level - alg/viewshed - viewshed.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 100 178 56.2 %
Date: 2025-07-09 17:50:03 Functions: 7 8 87.5 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Viewshed Generation
       4             :  * Purpose:  Core algorithm implementation for viewshed generation.
       5             :  * Author:   Tamas Szekeres, szekerest@gmail.com
       6             :  *
       7             :  * (c) 2024 info@hobu.co
       8             :  *
       9             :  ******************************************************************************
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include <algorithm>
      15             : #include <array>
      16             : 
      17             : #include "gdal_alg.h"
      18             : #include "gdal_priv_templates.hpp"
      19             : 
      20             : #include "progress.h"
      21             : #include "util.h"
      22             : #include "viewshed.h"
      23             : #include "viewshed_executor.h"
      24             : 
      25             : /************************************************************************/
      26             : /*                        GDALViewshedGenerate()                        */
      27             : /************************************************************************/
      28             : 
      29             : /**
      30             :  * Create viewshed from raster DEM.
      31             :  *
      32             :  * This algorithm will generate a viewshed raster from an input DEM raster
      33             :  * by using a modified algorithm of "Generating Viewsheds without Using
      34             :  * Sightlines" published at
      35             :  * https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
      36             :  * This approach provides a relatively fast calculation, since the output raster
      37             :  * is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
      38             :  * be used as an example of how to use this function. The output raster will be
      39             :  * of type Byte or Float64.
      40             :  *
      41             :  * \note The algorithm as implemented currently will only output meaningful
      42             :  * results if the georeferencing is in a projected coordinate reference system.
      43             :  *
      44             :  * @param hBand The band to read the DEM data from. Only the part of the raster
      45             :  * within the specified maxdistance around the observer point is processed.
      46             :  *
      47             :  * @param pszDriverName Driver name (GTiff if set to NULL)
      48             :  *
      49             :  * @param pszTargetRasterName The name of the target raster to be generated.
      50             :  * Must not be NULL
      51             :  *
      52             :  * @param papszCreationOptions creation options.
      53             :  *
      54             :  * @param dfObserverX observer X value (in SRS units)
      55             :  *
      56             :  * @param dfObserverY observer Y value (in SRS units)
      57             :  *
      58             :  * @param dfObserverHeight The height of the observer above the DEM surface.
      59             :  *
      60             :  * @param dfTargetHeight The height of the target above the DEM surface.
      61             :  * (default 0)
      62             :  *
      63             :  * @param dfVisibleVal pixel value for visibility (default 255)
      64             :  *
      65             :  * @param dfInvisibleVal pixel value for invisibility (default 0)
      66             :  *
      67             :  * @param dfOutOfRangeVal The value to be set for the cells that fall outside of
      68             :  * the range specified by dfMaxDistance.
      69             :  *
      70             :  * @param dfNoDataVal The value to be set for the cells that have no data.
      71             :  *                    If set to a negative value, nodata is not set.
      72             :  *                    Note: currently, no special processing of input cells at a
      73             :  * nodata value is done (which may result in erroneous results).
      74             :  *
      75             :  * @param dfCurvCoeff Coefficient to consider the effect of the curvature and
      76             :  * refraction. The height of the DEM is corrected according to the following
      77             :  * formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
      78             :  * the effect of the atmospheric refraction we can use 0.85714.
      79             :  *
      80             :  * @param eMode The mode of the viewshed calculation.
      81             :  * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
      82             :  * GVM_Min = 4.
      83             :  *
      84             :  * @param dfMaxDistance maximum distance range to compute viewshed.
      85             :  *                      It is also used to clamp the extent of the output
      86             :  * raster. If set to 0, then unlimited range is assumed, that is to say the
      87             :  *                      computation is performed on the extent of the whole
      88             :  * raster.
      89             :  *
      90             :  * @param pfnProgress A GDALProgressFunc that may be used to report progress
      91             :  * to the user, or to interrupt the algorithm.  May be NULL if not required.
      92             :  *
      93             :  * @param pProgressArg The callback data for the pfnProgress function.
      94             :  *
      95             :  * @param heightMode Type of information contained in output raster. Possible
      96             :  * values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
      97             :  *                   GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
      98             :  *
      99             :  *                   GVOT_NORMAL returns a raster of type Byte containing
     100             :  * visible locations.
     101             :  *
     102             :  *                   GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
     103             :  * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
     104             :  * containing the minimum target height for target to be visible from the DEM
     105             :  * surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
     106             :  * and dfInvisibleVal will be ignored.
     107             :  *
     108             :  *
     109             :  * @param papszExtraOptions Future extra options. Must be set to NULL currently.
     110             :  *
     111             :  * @return not NULL output dataset on success (to be closed with GDALClose()) or
     112             :  * NULL if an error occurs.
     113             :  *
     114             :  * @since GDAL 3.1
     115             :  */
     116           0 : GDALDatasetH GDALViewshedGenerate(
     117             :     GDALRasterBandH hBand, const char *pszDriverName,
     118             :     const char *pszTargetRasterName, CSLConstList papszCreationOptions,
     119             :     double dfObserverX, double dfObserverY, double dfObserverHeight,
     120             :     double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
     121             :     double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
     122             :     GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
     123             :     void *pProgressArg, GDALViewshedOutputType heightMode,
     124             :     [[maybe_unused]] CSLConstList papszExtraOptions)
     125             : {
     126             :     using namespace gdal;
     127             : 
     128           0 :     viewshed::Options oOpts;
     129           0 :     oOpts.outputFormat = pszDriverName;
     130           0 :     oOpts.outputFilename = pszTargetRasterName;
     131           0 :     oOpts.creationOpts = papszCreationOptions;
     132           0 :     oOpts.observer.x = dfObserverX;
     133           0 :     oOpts.observer.y = dfObserverY;
     134           0 :     oOpts.observer.z = dfObserverHeight;
     135           0 :     oOpts.targetHeight = dfTargetHeight;
     136           0 :     oOpts.curveCoeff = dfCurvCoeff;
     137           0 :     oOpts.maxDistance = dfMaxDistance;
     138           0 :     oOpts.nodataVal = dfNoDataVal;
     139             : 
     140           0 :     switch (eMode)
     141             :     {
     142           0 :         case GVM_Edge:
     143           0 :             oOpts.cellMode = viewshed::CellMode::Edge;
     144           0 :             break;
     145           0 :         case GVM_Diagonal:
     146           0 :             oOpts.cellMode = viewshed::CellMode::Diagonal;
     147           0 :             break;
     148           0 :         case GVM_Min:
     149           0 :             oOpts.cellMode = viewshed::CellMode::Min;
     150           0 :             break;
     151           0 :         case GVM_Max:
     152           0 :             oOpts.cellMode = viewshed::CellMode::Max;
     153           0 :             break;
     154             :     }
     155             : 
     156           0 :     switch (heightMode)
     157             :     {
     158           0 :         case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
     159           0 :             oOpts.outputMode = viewshed::OutputMode::DEM;
     160           0 :             break;
     161           0 :         case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
     162           0 :             oOpts.outputMode = viewshed::OutputMode::Ground;
     163           0 :             break;
     164           0 :         case GVOT_NORMAL:
     165           0 :             oOpts.outputMode = viewshed::OutputMode::Normal;
     166           0 :             break;
     167             :     }
     168             : 
     169           0 :     if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))
     170             :     {
     171           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     172             :                  "dfVisibleVal out of range. Must be [0, 255].");
     173           0 :         return nullptr;
     174             :     }
     175           0 :     if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))
     176             :     {
     177           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     178             :                  "dfInvisibleVal out of range. Must be [0, 255].");
     179           0 :         return nullptr;
     180             :     }
     181           0 :     if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))
     182             :     {
     183           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     184             :                  "dfOutOfRangeVal out of range. Must be [0, 255].");
     185           0 :         return nullptr;
     186             :     }
     187           0 :     oOpts.visibleVal = dfVisibleVal;
     188           0 :     oOpts.invisibleVal = dfInvisibleVal;
     189           0 :     oOpts.outOfRangeVal = dfOutOfRangeVal;
     190             : 
     191           0 :     gdal::viewshed::Viewshed v(oOpts);
     192             : 
     193           0 :     if (!pfnProgress)
     194           0 :         pfnProgress = GDALDummyProgress;
     195           0 :     v.run(hBand, pfnProgress, pProgressArg);
     196             : 
     197           0 :     return GDALDataset::FromHandle(v.output().release());
     198             : }
     199             : 
     200             : namespace gdal
     201             : {
     202             : namespace viewshed
     203             : {
     204             : 
     205             : namespace
     206             : {
     207             : 
     208          39 : bool getTransforms(GDALRasterBand &band, GDALGeoTransform &fwdTransform,
     209             :                    GDALGeoTransform &revTransform)
     210             : {
     211          39 :     band.GetDataset()->GetGeoTransform(fwdTransform);
     212          39 :     if (!fwdTransform.GetInverse(revTransform))
     213             :     {
     214           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
     215           0 :         return false;
     216             :     }
     217          39 :     return true;
     218             : }
     219             : 
     220             : /// Shrink the extent of a window to just cover the slice defined by rays from
     221             : /// (nX, nY) and [startAngle, endAngle]
     222             : ///
     223             : /// @param oOutExtent  Window to modify
     224             : /// @param nX  X coordinate of ray endpoint.
     225             : /// @param nY  Y coordinate of ray endpoint.
     226             : /// @param startAngle  Start angle of slice (standard mathmatics notion, in radians)
     227             : /// @param endAngle  End angle of slice (standard mathmatics notion, in radians)
     228          52 : void shrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
     229             :                            double startAngle, double endAngle)
     230             : {
     231             :     /// NOTE: This probably doesn't work when the observer is outside the raster and
     232             :     ///   needs to be enhanced for that case.
     233             : 
     234          52 :     if (startAngle == endAngle)
     235          37 :         return;
     236             : 
     237          15 :     Window win = oOutExtent;
     238             : 
     239             :     // Set the X boundaries for the angles
     240          15 :     int startAngleX = hIntersect(startAngle, nX, nY, win);
     241          15 :     int stopAngleX = hIntersect(endAngle, nX, nY, win);
     242             : 
     243          15 :     int xmax = nX;
     244          15 :     if (!rayBetween(startAngle, endAngle, 0))
     245             :     {
     246           7 :         xmax = std::max(xmax, startAngleX);
     247           7 :         xmax = std::max(xmax, stopAngleX);
     248             :         // Add one to xmax since we want one past the stop. [start, stop)
     249           7 :         oOutExtent.xStop = std::min(oOutExtent.xStop, xmax + 1);
     250             :     }
     251             : 
     252          15 :     int xmin = nX;
     253          15 :     if (!rayBetween(startAngle, endAngle, M_PI))
     254             :     {
     255           8 :         xmin = std::min(xmin, startAngleX);
     256           8 :         xmin = std::min(xmin, stopAngleX);
     257           8 :         oOutExtent.xStart = std::max(oOutExtent.xStart, xmin);
     258             :     }
     259             : 
     260             :     // Set the Y boundaries for the angles
     261          15 :     int startAngleY = vIntersect(startAngle, nX, nY, win);
     262          15 :     int stopAngleY = vIntersect(endAngle, nX, nY, win);
     263             : 
     264          15 :     int ymin = nY;
     265          15 :     if (!rayBetween(startAngle, endAngle, M_PI / 2))
     266             :     {
     267           7 :         ymin = std::min(ymin, startAngleY);
     268           7 :         ymin = std::min(ymin, stopAngleY);
     269           7 :         oOutExtent.yStart = std::max(oOutExtent.yStart, ymin);
     270             :     }
     271          15 :     int ymax = nY;
     272          15 :     if (!rayBetween(startAngle, endAngle, 3 * M_PI / 2))
     273             :     {
     274           8 :         ymax = std::max(ymax, startAngleY);
     275           8 :         ymax = std::max(ymax, stopAngleY);
     276             :         // Add one to ymax since we want one past the stop. [start, stop)
     277           8 :         oOutExtent.yStop = std::min(oOutExtent.yStop, ymax + 1);
     278             :     }
     279             : }
     280             : 
     281             : }  // unnamed namespace
     282             : 
     283          39 : Viewshed::Viewshed(const Options &opts) : oOpts{opts}
     284             : {
     285          39 : }
     286             : 
     287             : Viewshed::~Viewshed() = default;
     288             : 
     289             : /// Calculate the extent of the output raster in terms of the input raster and
     290             : /// save the input raster extent.
     291             : ///
     292             : /// @return  false on error, true otherwise
     293          39 : bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)
     294             : {
     295             :     // We start with the assumption that the output size matches the input.
     296          39 :     oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
     297          39 :     oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
     298             : 
     299          39 :     if (!oOutExtent.contains(nX, nY))
     300             :     {
     301           8 :         if (oOpts.startAngle != oOpts.endAngle)
     302             :         {
     303           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     304             :                      "Angle masking is not supported with an out-of-raster "
     305             :                      "observer.");
     306           0 :             return false;
     307             :         }
     308           8 :         CPLError(CE_Warning, CPLE_AppDefined,
     309             :                  "NOTE: The observer location falls outside of the DEM area");
     310             :     }
     311             : 
     312          39 :     constexpr double EPSILON = 1e-8;
     313          39 :     if (oOpts.maxDistance > 0)
     314             :     {
     315             :         //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
     316             :         //  Find the distance in the direction of the transformed unit vector in the X and Y
     317             :         //  directions and use those factors to determine the limiting values in the raster space.
     318           3 :         int nXStart = static_cast<int>(
     319           3 :             std::floor(nX - invGT[1] * oOpts.maxDistance + EPSILON));
     320           3 :         int nXStop = static_cast<int>(
     321           3 :             std::ceil(nX + invGT[1] * oOpts.maxDistance - EPSILON) + 1);
     322             :         //ABELL - These seem to be wrong. The transform of 1 is no transform, so not
     323             :         //  sure why we're adding one in the first case. Really, the transformed distance
     324             :         // should add EPSILON. Not sure what the change should be for a negative transform,
     325             :         // which is what I think is being handled with the 1/0 addition/subtraction.
     326             :         int nYStart =
     327           3 :             static_cast<int>(std::floor(
     328           3 :                 nY - std::fabs(invGT[5]) * oOpts.maxDistance + EPSILON)) -
     329           3 :             (invGT[5] > 0 ? 1 : 0);
     330           3 :         int nYStop = static_cast<int>(
     331           6 :             std::ceil(nY + std::fabs(invGT[5]) * oOpts.maxDistance - EPSILON) +
     332           3 :             (invGT[5] < 0 ? 1 : 0));
     333             : 
     334             :         // If the limits are invalid, set the window size to zero to trigger the error below.
     335           3 :         if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
     336           3 :             nYStart >= oOutExtent.yStop || nYStop < 0)
     337             :         {
     338           0 :             oOutExtent = Window();
     339             :         }
     340             :         else
     341             :         {
     342           3 :             oOutExtent.xStart = std::max(nXStart, 0);
     343           3 :             oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
     344             : 
     345           3 :             oOutExtent.yStart = std::max(nYStart, 0);
     346           3 :             oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
     347             :         }
     348             :     }
     349             : 
     350          39 :     if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
     351             :     {
     352           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     353             :                  "Invalid target raster size due to transform "
     354             :                  "and/or distance limitation.");
     355           0 :         return false;
     356             :     }
     357             : 
     358          39 :     shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);
     359             : 
     360             :     // normalize horizontal index to [ 0, oOutExtent.xSize() )
     361          39 :     oCurExtent = oOutExtent;
     362          39 :     oCurExtent.shiftX(-oOutExtent.xStart);
     363             : 
     364          39 :     return true;
     365             : }
     366             : 
     367             : /// Compute the viewshed of a raster band.
     368             : ///
     369             : /// @param band  Pointer to the raster band to be processed.
     370             : /// @param pfnProgress  Pointer to the progress function. Can be null.
     371             : /// @param pProgressArg  Argument passed to the progress function
     372             : /// @return  True on success, false otherwise.
     373          39 : bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
     374             :                    void *pProgressArg)
     375             : {
     376          39 :     pSrcBand = static_cast<GDALRasterBand *>(band);
     377             : 
     378          39 :     GDALGeoTransform fwdTransform, invTransform;
     379          39 :     if (!getTransforms(*pSrcBand, fwdTransform, invTransform))
     380           0 :         return false;
     381             : 
     382             :     // calculate observer position
     383             :     double dfX, dfY;
     384          39 :     invTransform.Apply(oOpts.observer.x, oOpts.observer.y, &dfX, &dfY);
     385          39 :     if (!GDALIsValueInRange<int>(dfX))
     386             :     {
     387           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
     388           0 :         return false;
     389             :     }
     390          39 :     if (!GDALIsValueInRange<int>(dfY))
     391             :     {
     392           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
     393           0 :         return false;
     394             :     }
     395          39 :     int nX = static_cast<int>(dfX);
     396          39 :     int nY = static_cast<int>(dfY);
     397             : 
     398          39 :     if (oOpts.startAngle < 0 || oOpts.startAngle >= 360)
     399             :     {
     400           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     401             :                  "Start angle out of range. Must be [0, 360).");
     402           0 :         return false;
     403             :     }
     404          39 :     if (oOpts.endAngle < 0 || oOpts.endAngle >= 360)
     405             :     {
     406           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     407             :                  "End angle out of range. Must be [0, 360).");
     408           0 :         return false;
     409             :     }
     410          39 :     if (oOpts.highPitch > 90)
     411             :     {
     412           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     413             :                  "Invalid highPitch. Cannot be greater than 90.");
     414           0 :         return false;
     415             :     }
     416          39 :     if (oOpts.lowPitch < -90)
     417             :     {
     418           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     419             :                  "Invalid lowPitch. Cannot be less than -90.");
     420           0 :         return false;
     421             :     }
     422          39 :     if (oOpts.highPitch <= oOpts.lowPitch)
     423             :     {
     424           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     425             :                  "Invalid pitch. highPitch must be > lowPitch");
     426           0 :         return false;
     427             :     }
     428             : 
     429             :     // Normalize angle to radians and standard math arrangement.
     430          39 :     oOpts.startAngle = normalizeAngle(oOpts.startAngle);
     431          39 :     oOpts.endAngle = normalizeAngle(oOpts.endAngle);
     432             : 
     433             :     // Must calculate extents in order to make the output dataset.
     434          39 :     if (!calcExtents(nX, nY, invTransform))
     435           0 :         return false;
     436             : 
     437          39 :     poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
     438          39 :     if (!poDstDS)
     439           0 :         return false;
     440             : 
     441             :     // Create the progress reporter.
     442          78 :     Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());
     443             : 
     444             :     // Execute the viewshed algorithm.
     445          39 :     GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
     446          39 :     ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
     447          39 :                               oCurExtent, oOpts, oProgress,
     448          78 :                               /* emitWarningIfNoData = */ true);
     449          39 :     executor.run();
     450          39 :     oProgress.emit(1);
     451          39 :     return static_cast<bool>(poDstDS);
     452             : }
     453             : 
     454             : // Adjust the coefficient of curvature for non-earth SRS.
     455             : /// \param curveCoeff  Current curve coefficient
     456             : /// \param hSrcDS  Source dataset
     457             : /// \return  Adjusted curve coefficient.
     458          14 : double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)
     459             : {
     460             :     const OGRSpatialReference *poSRS =
     461          14 :         GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
     462          14 :     if (poSRS)
     463             :     {
     464             :         OGRErr eSRSerr;
     465          13 :         const double dfSemiMajor = poSRS->GetSemiMajor(&eSRSerr);
     466          13 :         if (eSRSerr != OGRERR_FAILURE &&
     467          13 :             fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) >
     468             :                 0.05 * SRS_WGS84_SEMIMAJOR)
     469             :         {
     470           0 :             curveCoeff = 1.0;
     471           0 :             CPLDebug("gdal_viewshed",
     472             :                      "Using -cc=1.0 as a non-Earth CRS has been detected");
     473             :         }
     474             :     }
     475          14 :     return curveCoeff;
     476             : }
     477             : 
     478             : #if defined(__clang__) || defined(__GNUC__)
     479             : #pragma GCC diagnostic ignored "-Wmissing-declarations"
     480             : #endif
     481             : 
     482          13 : void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
     483             :                                double startAngle, double endAngle)
     484             : {
     485          13 :     shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);
     486          13 : }
     487             : 
     488             : }  // namespace viewshed
     489             : }  // namespace gdal

Generated by: LCOV version 1.14