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 : * @param papszExtraOptions Extra options to control the viewshed analysis.
109 : * This is a NULL-terminated list of strings in "KEY=VALUE" format, or NULL for no options.
110 : * The following keys are supported:
111 : * <ul>
112 : * <li><b>START_ANGLE</b>: Mask all cells outside of the arc ('start-angle', 'end-angle'). Clockwise degrees from north. Also used to clamp the extent of the output raster.</li>
113 : * <li><b>END_ANGLE</b>: Mask all cells outside of the arc ('start-angle', 'end-angle'). Clockwise degrees from north. Also used to clamp the extent of the output raster.</li>
114 : * <li><b>LOW_PITCH</b>: Bound observable height to be no lower than the 'low-pitch' angle from the observer. Degrees from horizontal - positive is up. Must be less than 'high-pitch'.</li>
115 : * <li><b>HIGH_PITCH</b>: Mark all cells out-of-range where the observable height would be higher than the 'high-pitch' angle from the observer. Degrees from horizontal - positive is up. Must be greater than 'low-pitch'.</li>
116 : * </ul>
117 : * If NULL, a 360-degree viewshed is calculated.
118 : *
119 : * @return not NULL output dataset on success (to be closed with GDALClose()) or
120 : * NULL if an error occurs.
121 : *
122 : * @since GDAL 3.1
123 : */
124 0 : GDALDatasetH GDALViewshedGenerate(
125 : GDALRasterBandH hBand, const char *pszDriverName,
126 : const char *pszTargetRasterName, CSLConstList papszCreationOptions,
127 : double dfObserverX, double dfObserverY, double dfObserverHeight,
128 : double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
129 : double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
130 : GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
131 : void *pProgressArg, GDALViewshedOutputType heightMode,
132 : [[maybe_unused]] CSLConstList papszExtraOptions)
133 : {
134 : using namespace gdal;
135 :
136 0 : viewshed::Options oOpts;
137 0 : oOpts.outputFormat = pszDriverName;
138 0 : oOpts.outputFilename = pszTargetRasterName;
139 0 : oOpts.creationOpts = papszCreationOptions;
140 0 : oOpts.observer.x = dfObserverX;
141 0 : oOpts.observer.y = dfObserverY;
142 0 : oOpts.observer.z = dfObserverHeight;
143 0 : oOpts.targetHeight = dfTargetHeight;
144 0 : oOpts.curveCoeff = dfCurvCoeff;
145 0 : oOpts.maxDistance = dfMaxDistance;
146 0 : oOpts.nodataVal = dfNoDataVal;
147 :
148 0 : switch (eMode)
149 : {
150 0 : case GVM_Edge:
151 0 : oOpts.cellMode = viewshed::CellMode::Edge;
152 0 : break;
153 0 : case GVM_Diagonal:
154 0 : oOpts.cellMode = viewshed::CellMode::Diagonal;
155 0 : break;
156 0 : case GVM_Min:
157 0 : oOpts.cellMode = viewshed::CellMode::Min;
158 0 : break;
159 0 : case GVM_Max:
160 0 : oOpts.cellMode = viewshed::CellMode::Max;
161 0 : break;
162 : }
163 :
164 0 : switch (heightMode)
165 : {
166 0 : case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
167 0 : oOpts.outputMode = viewshed::OutputMode::DEM;
168 0 : break;
169 0 : case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
170 0 : oOpts.outputMode = viewshed::OutputMode::Ground;
171 0 : break;
172 0 : case GVOT_NORMAL:
173 0 : oOpts.outputMode = viewshed::OutputMode::Normal;
174 0 : break;
175 : }
176 :
177 0 : if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))
178 : {
179 0 : CPLError(CE_Failure, CPLE_AppDefined,
180 : "dfVisibleVal out of range. Must be [0, 255].");
181 0 : return nullptr;
182 : }
183 0 : if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))
184 : {
185 0 : CPLError(CE_Failure, CPLE_AppDefined,
186 : "dfInvisibleVal out of range. Must be [0, 255].");
187 0 : return nullptr;
188 : }
189 0 : if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))
190 : {
191 0 : CPLError(CE_Failure, CPLE_AppDefined,
192 : "dfOutOfRangeVal out of range. Must be [0, 255].");
193 0 : return nullptr;
194 : }
195 0 : oOpts.visibleVal = dfVisibleVal;
196 0 : oOpts.invisibleVal = dfInvisibleVal;
197 0 : oOpts.outOfRangeVal = dfOutOfRangeVal;
198 :
199 : const char *pszStartAngle =
200 0 : CSLFetchNameValue(papszExtraOptions, "START_ANGLE");
201 0 : if (pszStartAngle)
202 0 : oOpts.startAngle = CPLAtof(pszStartAngle);
203 :
204 0 : const char *pszEndAngle = CSLFetchNameValue(papszExtraOptions, "END_ANGLE");
205 0 : if (pszEndAngle)
206 0 : oOpts.endAngle = CPLAtof(pszEndAngle);
207 :
208 0 : const char *pszLowPitch = CSLFetchNameValue(papszExtraOptions, "LOW_PITCH");
209 0 : if (pszLowPitch)
210 0 : oOpts.lowPitch = CPLAtof(pszLowPitch);
211 :
212 : const char *pszHighPitch =
213 0 : CSLFetchNameValue(papszExtraOptions, "HIGH_PITCH");
214 0 : if (pszHighPitch)
215 0 : oOpts.highPitch = CPLAtof(pszHighPitch);
216 :
217 0 : gdal::viewshed::Viewshed v(oOpts);
218 :
219 0 : if (!pfnProgress)
220 0 : pfnProgress = GDALDummyProgress;
221 0 : v.run(hBand, pfnProgress, pProgressArg);
222 :
223 0 : return GDALDataset::FromHandle(v.output().release());
224 : }
225 :
226 : namespace gdal
227 : {
228 : namespace viewshed
229 : {
230 :
231 : namespace
232 : {
233 :
234 39 : bool getTransforms(GDALRasterBand &band, GDALGeoTransform &fwdTransform,
235 : GDALGeoTransform &revTransform)
236 : {
237 39 : band.GetDataset()->GetGeoTransform(fwdTransform);
238 39 : if (!fwdTransform.GetInverse(revTransform))
239 : {
240 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
241 0 : return false;
242 : }
243 39 : return true;
244 : }
245 :
246 : /// Shrink the extent of a window to just cover the slice defined by rays from
247 : /// (nX, nY) and [startAngle, endAngle]
248 : ///
249 : /// @param oOutExtent Window to modify
250 : /// @param nX X coordinate of ray endpoint.
251 : /// @param nY Y coordinate of ray endpoint.
252 : /// @param startAngle Start angle of slice (standard mathematics notion, in radians)
253 : /// @param endAngle End angle of slice (standard mathematics notion, in radians)
254 52 : void shrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
255 : double startAngle, double endAngle)
256 : {
257 : /// NOTE: This probably doesn't work when the observer is outside the raster and
258 : /// needs to be enhanced for that case.
259 :
260 52 : if (startAngle == endAngle)
261 37 : return;
262 :
263 15 : Window win = oOutExtent;
264 :
265 : // Set the X boundaries for the angles
266 15 : int startAngleX = hIntersect(startAngle, nX, nY, win);
267 15 : int stopAngleX = hIntersect(endAngle, nX, nY, win);
268 :
269 15 : int xmax = nX;
270 15 : if (!rayBetween(startAngle, endAngle, 0))
271 : {
272 7 : xmax = std::max(xmax, startAngleX);
273 7 : xmax = std::max(xmax, stopAngleX);
274 : // Add one to xmax since we want one past the stop. [start, stop)
275 7 : oOutExtent.xStop = std::min(oOutExtent.xStop, xmax + 1);
276 : }
277 :
278 15 : int xmin = nX;
279 15 : if (!rayBetween(startAngle, endAngle, M_PI))
280 : {
281 8 : xmin = std::min(xmin, startAngleX);
282 8 : xmin = std::min(xmin, stopAngleX);
283 8 : oOutExtent.xStart = std::max(oOutExtent.xStart, xmin);
284 : }
285 :
286 : // Set the Y boundaries for the angles
287 15 : int startAngleY = vIntersect(startAngle, nX, nY, win);
288 15 : int stopAngleY = vIntersect(endAngle, nX, nY, win);
289 :
290 15 : int ymin = nY;
291 15 : if (!rayBetween(startAngle, endAngle, M_PI / 2))
292 : {
293 7 : ymin = std::min(ymin, startAngleY);
294 7 : ymin = std::min(ymin, stopAngleY);
295 7 : oOutExtent.yStart = std::max(oOutExtent.yStart, ymin);
296 : }
297 15 : int ymax = nY;
298 15 : if (!rayBetween(startAngle, endAngle, 3 * M_PI / 2))
299 : {
300 8 : ymax = std::max(ymax, startAngleY);
301 8 : ymax = std::max(ymax, stopAngleY);
302 : // Add one to ymax since we want one past the stop. [start, stop)
303 8 : oOutExtent.yStop = std::min(oOutExtent.yStop, ymax + 1);
304 : }
305 : }
306 :
307 : } // unnamed namespace
308 :
309 39 : Viewshed::Viewshed(const Options &opts) : oOpts{opts}
310 : {
311 39 : }
312 :
313 : Viewshed::~Viewshed() = default;
314 :
315 : /// Calculate the extent of the output raster in terms of the input raster and
316 : /// save the input raster extent.
317 : ///
318 : /// @return false on error, true otherwise
319 39 : bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)
320 : {
321 : // We start with the assumption that the output size matches the input.
322 39 : oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
323 39 : oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
324 :
325 39 : if (!oOutExtent.contains(nX, nY))
326 : {
327 8 : if (oOpts.startAngle != oOpts.endAngle)
328 : {
329 0 : CPLError(CE_Failure, CPLE_AppDefined,
330 : "Angle masking is not supported with an out-of-raster "
331 : "observer.");
332 0 : return false;
333 : }
334 8 : CPLError(CE_Warning, CPLE_AppDefined,
335 : "NOTE: The observer location falls outside of the DEM area");
336 : }
337 :
338 39 : constexpr double EPSILON = 1e-8;
339 39 : if (oOpts.maxDistance > 0)
340 : {
341 : //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
342 : // Find the distance in the direction of the transformed unit vector in the X and Y
343 : // directions and use those factors to determine the limiting values in the raster space.
344 3 : int nXStart = static_cast<int>(
345 3 : std::floor(nX - invGT[1] * oOpts.maxDistance + EPSILON));
346 3 : int nXStop = static_cast<int>(
347 3 : std::ceil(nX + invGT[1] * oOpts.maxDistance - EPSILON) + 1);
348 : //ABELL - These seem to be wrong. The transform of 1 is no transform, so not
349 : // sure why we're adding one in the first case. Really, the transformed distance
350 : // should add EPSILON. Not sure what the change should be for a negative transform,
351 : // which is what I think is being handled with the 1/0 addition/subtraction.
352 : int nYStart =
353 3 : static_cast<int>(std::floor(
354 3 : nY - std::fabs(invGT[5]) * oOpts.maxDistance + EPSILON)) -
355 3 : (invGT[5] > 0 ? 1 : 0);
356 3 : int nYStop = static_cast<int>(
357 6 : std::ceil(nY + std::fabs(invGT[5]) * oOpts.maxDistance - EPSILON) +
358 3 : (invGT[5] < 0 ? 1 : 0));
359 :
360 : // If the limits are invalid, set the window size to zero to trigger the error below.
361 3 : if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
362 3 : nYStart >= oOutExtent.yStop || nYStop < 0)
363 : {
364 0 : oOutExtent = Window();
365 : }
366 : else
367 : {
368 3 : oOutExtent.xStart = std::max(nXStart, 0);
369 3 : oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
370 :
371 3 : oOutExtent.yStart = std::max(nYStart, 0);
372 3 : oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
373 : }
374 : }
375 :
376 39 : if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
377 : {
378 0 : CPLError(CE_Failure, CPLE_AppDefined,
379 : "Invalid target raster size due to transform "
380 : "and/or distance limitation.");
381 0 : return false;
382 : }
383 :
384 39 : shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);
385 :
386 : // normalize horizontal index to [ 0, oOutExtent.xSize() )
387 39 : oCurExtent = oOutExtent;
388 39 : oCurExtent.shiftX(-oOutExtent.xStart);
389 :
390 39 : return true;
391 : }
392 :
393 : /// Compute the viewshed of a raster band.
394 : ///
395 : /// @param band Pointer to the raster band to be processed.
396 : /// @param pfnProgress Pointer to the progress function. Can be null.
397 : /// @param pProgressArg Argument passed to the progress function
398 : /// @return True on success, false otherwise.
399 39 : bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
400 : void *pProgressArg)
401 : {
402 39 : pSrcBand = static_cast<GDALRasterBand *>(band);
403 :
404 39 : GDALGeoTransform fwdTransform, invTransform;
405 39 : if (!getTransforms(*pSrcBand, fwdTransform, invTransform))
406 0 : return false;
407 :
408 : // calculate observer position
409 : double dfX, dfY;
410 39 : invTransform.Apply(oOpts.observer.x, oOpts.observer.y, &dfX, &dfY);
411 39 : if (!GDALIsValueInRange<int>(dfX))
412 : {
413 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
414 0 : return false;
415 : }
416 39 : if (!GDALIsValueInRange<int>(dfY))
417 : {
418 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
419 0 : return false;
420 : }
421 39 : int nX = static_cast<int>(dfX);
422 39 : int nY = static_cast<int>(dfY);
423 :
424 39 : if (oOpts.startAngle < 0 || oOpts.startAngle >= 360)
425 : {
426 0 : CPLError(CE_Failure, CPLE_AppDefined,
427 : "Start angle out of range. Must be [0, 360).");
428 0 : return false;
429 : }
430 39 : if (oOpts.endAngle < 0 || oOpts.endAngle >= 360)
431 : {
432 0 : CPLError(CE_Failure, CPLE_AppDefined,
433 : "End angle out of range. Must be [0, 360).");
434 0 : return false;
435 : }
436 39 : if (oOpts.highPitch > 90)
437 : {
438 0 : CPLError(CE_Failure, CPLE_AppDefined,
439 : "Invalid highPitch. Cannot be greater than 90.");
440 0 : return false;
441 : }
442 39 : if (oOpts.lowPitch < -90)
443 : {
444 0 : CPLError(CE_Failure, CPLE_AppDefined,
445 : "Invalid lowPitch. Cannot be less than -90.");
446 0 : return false;
447 : }
448 39 : if (oOpts.highPitch <= oOpts.lowPitch)
449 : {
450 0 : CPLError(CE_Failure, CPLE_AppDefined,
451 : "Invalid pitch. highPitch must be > lowPitch");
452 0 : return false;
453 : }
454 :
455 : // Normalize angle to radians and standard math arrangement.
456 39 : oOpts.startAngle = normalizeAngle(oOpts.startAngle);
457 39 : oOpts.endAngle = normalizeAngle(oOpts.endAngle);
458 :
459 : // Must calculate extents in order to make the output dataset.
460 39 : if (!calcExtents(nX, nY, invTransform))
461 0 : return false;
462 :
463 39 : poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
464 39 : if (!poDstDS)
465 0 : return false;
466 :
467 : // Create the progress reporter.
468 78 : Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());
469 :
470 : // Execute the viewshed algorithm.
471 39 : GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
472 39 : ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
473 39 : oCurExtent, oOpts, oProgress,
474 78 : /* emitWarningIfNoData = */ true);
475 39 : executor.run();
476 39 : oProgress.emit(1);
477 39 : return static_cast<bool>(poDstDS);
478 : }
479 :
480 : // Adjust the coefficient of curvature for non-earth SRS.
481 : /// \param curveCoeff Current curve coefficient
482 : /// \param hSrcDS Source dataset
483 : /// \return Adjusted curve coefficient.
484 14 : double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)
485 : {
486 : const OGRSpatialReference *poSRS =
487 14 : GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
488 14 : if (poSRS)
489 : {
490 : OGRErr eSRSerr;
491 13 : const double dfSemiMajor = poSRS->GetSemiMajor(&eSRSerr);
492 13 : if (eSRSerr != OGRERR_FAILURE &&
493 13 : fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) >
494 : 0.05 * SRS_WGS84_SEMIMAJOR)
495 : {
496 0 : curveCoeff = 1.0;
497 0 : CPLDebug("gdal_viewshed",
498 : "Using -cc=1.0 as a non-Earth CRS has been detected");
499 : }
500 : }
501 14 : return curveCoeff;
502 : }
503 :
504 : #if defined(__clang__) || defined(__GNUC__)
505 : #pragma GCC diagnostic ignored "-Wmissing-declarations"
506 : #endif
507 :
508 13 : void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
509 : double startAngle, double endAngle)
510 : {
511 13 : shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);
512 13 : }
513 :
514 : } // namespace viewshed
515 : } // namespace gdal
|
