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, double *pFwdTransform,
209 : double *pRevTransform)
210 : {
211 39 : band.GetDataset()->GetGeoTransform(pFwdTransform);
212 39 : if (!GDALInvGeoTransform(pFwdTransform, pRevTransform))
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,
294 : const std::array<double, 6> &adfInvTransform)
295 : {
296 : // We start with the assumption that the output size matches the input.
297 39 : oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
298 39 : oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
299 :
300 39 : if (!oOutExtent.contains(nX, nY))
301 : {
302 8 : if (oOpts.startAngle != oOpts.endAngle)
303 : {
304 0 : CPLError(CE_Failure, CPLE_AppDefined,
305 : "Angle masking is not supported with an out-of-raster "
306 : "observer.");
307 0 : return false;
308 : }
309 8 : CPLError(CE_Warning, CPLE_AppDefined,
310 : "NOTE: The observer location falls outside of the DEM area");
311 : }
312 :
313 39 : constexpr double EPSILON = 1e-8;
314 39 : if (oOpts.maxDistance > 0)
315 : {
316 : //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
317 : // Find the distance in the direction of the transformed unit vector in the X and Y
318 : // directions and use those factors to determine the limiting values in the raster space.
319 3 : int nXStart = static_cast<int>(
320 3 : std::floor(nX - adfInvTransform[1] * oOpts.maxDistance + EPSILON));
321 3 : int nXStop = static_cast<int>(
322 3 : std::ceil(nX + adfInvTransform[1] * oOpts.maxDistance - EPSILON) +
323 3 : 1);
324 : //ABELL - These seem to be wrong. The transform of 1 is no transform, so not
325 : // sure why we're adding one in the first case. Really, the transformed distance
326 : // should add EPSILON. Not sure what the change should be for a negative transform,
327 : // which is what I think is being handled with the 1/0 addition/subtraction.
328 : int nYStart =
329 3 : static_cast<int>(std::floor(
330 3 : nY - std::fabs(adfInvTransform[5]) * oOpts.maxDistance +
331 3 : EPSILON)) -
332 3 : (adfInvTransform[5] > 0 ? 1 : 0);
333 3 : int nYStop = static_cast<int>(
334 3 : std::ceil(nY + std::fabs(adfInvTransform[5]) * oOpts.maxDistance -
335 3 : EPSILON) +
336 3 : (adfInvTransform[5] < 0 ? 1 : 0));
337 :
338 : // If the limits are invalid, set the window size to zero to trigger the error below.
339 3 : if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
340 3 : nYStart >= oOutExtent.yStop || nYStop < 0)
341 : {
342 0 : oOutExtent = Window();
343 : }
344 : else
345 : {
346 3 : oOutExtent.xStart = std::max(nXStart, 0);
347 3 : oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
348 :
349 3 : oOutExtent.yStart = std::max(nYStart, 0);
350 3 : oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
351 : }
352 : }
353 :
354 39 : if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
355 : {
356 0 : CPLError(CE_Failure, CPLE_AppDefined,
357 : "Invalid target raster size due to transform "
358 : "and/or distance limitation.");
359 0 : return false;
360 : }
361 :
362 39 : shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);
363 :
364 : // normalize horizontal index to [ 0, oOutExtent.xSize() )
365 39 : oCurExtent = oOutExtent;
366 39 : oCurExtent.shiftX(-oOutExtent.xStart);
367 :
368 39 : return true;
369 : }
370 :
371 : /// Compute the viewshed of a raster band.
372 : ///
373 : /// @param band Pointer to the raster band to be processed.
374 : /// @param pfnProgress Pointer to the progress function. Can be null.
375 : /// @param pProgressArg Argument passed to the progress function
376 : /// @return True on success, false otherwise.
377 39 : bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
378 : void *pProgressArg)
379 : {
380 39 : pSrcBand = static_cast<GDALRasterBand *>(band);
381 :
382 : std::array<double, 6> adfFwdTransform;
383 : std::array<double, 6> adfInvTransform;
384 39 : if (!getTransforms(*pSrcBand, adfFwdTransform.data(),
385 : adfInvTransform.data()))
386 0 : return false;
387 :
388 : // calculate observer position
389 : double dfX, dfY;
390 39 : GDALApplyGeoTransform(adfInvTransform.data(), oOpts.observer.x,
391 : oOpts.observer.y, &dfX, &dfY);
392 39 : if (!GDALIsValueInRange<int>(dfX))
393 : {
394 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
395 0 : return false;
396 : }
397 39 : if (!GDALIsValueInRange<int>(dfY))
398 : {
399 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
400 0 : return false;
401 : }
402 39 : int nX = static_cast<int>(dfX);
403 39 : int nY = static_cast<int>(dfY);
404 :
405 39 : if (oOpts.startAngle < 0 || oOpts.startAngle >= 360)
406 : {
407 0 : CPLError(CE_Failure, CPLE_AppDefined,
408 : "Start angle out of range. Must be [0, 360).");
409 0 : return false;
410 : }
411 39 : if (oOpts.endAngle < 0 || oOpts.endAngle >= 360)
412 : {
413 0 : CPLError(CE_Failure, CPLE_AppDefined,
414 : "End angle out of range. Must be [0, 360).");
415 0 : return false;
416 : }
417 39 : if (oOpts.highPitch > 90)
418 : {
419 0 : CPLError(CE_Failure, CPLE_AppDefined,
420 : "Invalid highPitch. Cannot be greater than 90.");
421 0 : return false;
422 : }
423 39 : if (oOpts.lowPitch < -90)
424 : {
425 0 : CPLError(CE_Failure, CPLE_AppDefined,
426 : "Invalid lowPitch. Cannot be less than -90.");
427 0 : return false;
428 : }
429 39 : if (oOpts.highPitch <= oOpts.lowPitch)
430 : {
431 0 : CPLError(CE_Failure, CPLE_AppDefined,
432 : "Invalid pitch. highPitch must be > lowPitch");
433 0 : return false;
434 : }
435 :
436 : // Normalize angle to radians and standard math arrangement.
437 39 : oOpts.startAngle = normalizeAngle(oOpts.startAngle);
438 39 : oOpts.endAngle = normalizeAngle(oOpts.endAngle);
439 :
440 : // Must calculate extents in order to make the output dataset.
441 39 : if (!calcExtents(nX, nY, adfInvTransform))
442 0 : return false;
443 :
444 39 : poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
445 39 : if (!poDstDS)
446 0 : return false;
447 :
448 : // Create the progress reporter.
449 78 : Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());
450 :
451 : // Execute the viewshed algorithm.
452 39 : GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
453 39 : ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
454 39 : oCurExtent, oOpts, oProgress,
455 78 : /* emitWarningIfNoData = */ true);
456 39 : executor.run();
457 39 : oProgress.emit(1);
458 39 : return static_cast<bool>(poDstDS);
459 : }
460 :
461 : // Adjust the coefficient of curvature for non-earth SRS.
462 : /// \param curveCoeff Current curve coefficient
463 : /// \param hSrcDS Source dataset
464 : /// \return Adjusted curve coefficient.
465 14 : double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)
466 : {
467 : const OGRSpatialReference *poSRS =
468 14 : GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
469 14 : if (poSRS)
470 : {
471 : OGRErr eSRSerr;
472 13 : const double dfSemiMajor = poSRS->GetSemiMajor(&eSRSerr);
473 13 : if (eSRSerr != OGRERR_FAILURE &&
474 13 : fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) >
475 : 0.05 * SRS_WGS84_SEMIMAJOR)
476 : {
477 0 : curveCoeff = 1.0;
478 0 : CPLDebug("gdal_viewshed",
479 : "Using -cc=1.0 as a non-Earth CRS has been detected");
480 : }
481 : }
482 14 : return curveCoeff;
483 : }
484 :
485 : #if defined(__clang__) || defined(__GNUC__)
486 : #pragma GCC diagnostic ignored "-Wmissing-declarations"
487 : #endif
488 :
489 13 : void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
490 : double startAngle, double endAngle)
491 : {
492 13 : shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);
493 13 : }
494 :
495 : } // namespace viewshed
496 : } // namespace gdal
|
