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
|
