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 :
16 : #include "gdal_alg.h"
17 : #include "gdal_priv_templates.hpp"
18 :
19 : #include "progress.h"
20 : #include "util.h"
21 : #include "viewshed.h"
22 : #include "viewshed_executor.h"
23 :
24 : /************************************************************************/
25 : /* GDALViewshedGenerate() */
26 : /************************************************************************/
27 :
28 : /**
29 : * Create viewshed from raster DEM.
30 : *
31 : * This algorithm will generate a viewshed raster from an input DEM raster
32 : * by using a modified algorithm of "Generating Viewsheds without Using
33 : * Sightlines" published at
34 : * https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
35 : * This approach provides a relatively fast calculation, since the output raster
36 : * is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
37 : * be used as an example of how to use this function. The output raster will be
38 : * of type Byte or Float64.
39 : *
40 : * \note The algorithm as implemented currently will only output meaningful
41 : * results if the georeferencing is in a projected coordinate reference system.
42 : *
43 : * @param hBand The band to read the DEM data from. Only the part of the raster
44 : * within the specified maxdistance around the observer point is processed.
45 : *
46 : * @param pszDriverName Driver name (GTiff if set to NULL)
47 : *
48 : * @param pszTargetRasterName The name of the target raster to be generated.
49 : * Must not be NULL
50 : *
51 : * @param papszCreationOptions creation options.
52 : *
53 : * @param dfObserverX observer X value (in SRS units)
54 : *
55 : * @param dfObserverY observer Y value (in SRS units)
56 : *
57 : * @param dfObserverHeight The height of the observer above the DEM surface.
58 : *
59 : * @param dfTargetHeight The height of the target above the DEM surface.
60 : * (default 0)
61 : *
62 : * @param dfVisibleVal pixel value for visibility (default 255)
63 : *
64 : * @param dfInvisibleVal pixel value for invisibility (default 0)
65 : *
66 : * @param dfOutOfRangeVal The value to be set for the cells that fall outside of
67 : * the range specified by dfMaxDistance.
68 : *
69 : * @param dfNoDataVal The value to be set for the cells that have no data.
70 : * If set to a negative value, nodata is not set.
71 : * Note: currently, no special processing of input cells at a
72 : * nodata value is done (which may result in erroneous results).
73 : *
74 : * @param dfCurvCoeff Coefficient to consider the effect of the curvature and
75 : * refraction. The height of the DEM is corrected according to the following
76 : * formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
77 : * the effect of the atmospheric refraction we can use 0.85714.
78 : *
79 : * @param eMode The mode of the viewshed calculation.
80 : * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
81 : * GVM_Min = 4.
82 : *
83 : * @param dfMaxDistance maximum distance range to compute viewshed.
84 : * It is also used to clamp the extent of the output
85 : * raster. If set to 0, then unlimited range is assumed, that is to say the
86 : * computation is performed on the extent of the whole
87 : * raster.
88 : *
89 : * @param pfnProgress A GDALProgressFunc that may be used to report progress
90 : * to the user, or to interrupt the algorithm. May be NULL if not required.
91 : *
92 : * @param pProgressArg The callback data for the pfnProgress function.
93 : *
94 : * @param heightMode Type of information contained in output raster. Possible
95 : * values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
96 : * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
97 : *
98 : * GVOT_NORMAL returns a raster of type Byte containing
99 : * visible locations.
100 : *
101 : * GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
102 : * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
103 : * containing the minimum target height for target to be visible from the DEM
104 : * surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
105 : * and dfInvisibleVal will be ignored.
106 : *
107 : *
108 : * @param papszExtraOptions Future extra options. Must be set to NULL currently.
109 : *
110 : * @return not NULL output dataset on success (to be closed with GDALClose()) or
111 : * NULL if an error occurs.
112 : *
113 : * @since GDAL 3.1
114 : */
115 1 : GDALDatasetH GDALViewshedGenerate(
116 : GDALRasterBandH hBand, const char *pszDriverName,
117 : const char *pszTargetRasterName, CSLConstList papszCreationOptions,
118 : double dfObserverX, double dfObserverY, double dfObserverHeight,
119 : double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
120 : double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
121 : GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
122 : void *pProgressArg, GDALViewshedOutputType heightMode,
123 : [[maybe_unused]] CSLConstList papszExtraOptions)
124 : {
125 : using namespace gdal;
126 :
127 2 : viewshed::Options oOpts;
128 1 : oOpts.outputFormat = pszDriverName;
129 1 : oOpts.outputFilename = pszTargetRasterName;
130 1 : oOpts.creationOpts = papszCreationOptions;
131 1 : oOpts.observer.x = dfObserverX;
132 1 : oOpts.observer.y = dfObserverY;
133 1 : oOpts.observer.z = dfObserverHeight;
134 1 : oOpts.targetHeight = dfTargetHeight;
135 1 : oOpts.curveCoeff = dfCurvCoeff;
136 1 : oOpts.maxDistance = dfMaxDistance;
137 1 : oOpts.nodataVal = dfNoDataVal;
138 :
139 1 : switch (eMode)
140 : {
141 1 : case GVM_Edge:
142 1 : oOpts.cellMode = viewshed::CellMode::Edge;
143 1 : break;
144 0 : case GVM_Diagonal:
145 0 : oOpts.cellMode = viewshed::CellMode::Diagonal;
146 0 : break;
147 0 : case GVM_Min:
148 0 : oOpts.cellMode = viewshed::CellMode::Min;
149 0 : break;
150 0 : case GVM_Max:
151 0 : oOpts.cellMode = viewshed::CellMode::Max;
152 0 : break;
153 : }
154 :
155 1 : switch (heightMode)
156 : {
157 0 : case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
158 0 : oOpts.outputMode = viewshed::OutputMode::DEM;
159 0 : break;
160 1 : case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
161 1 : oOpts.outputMode = viewshed::OutputMode::Ground;
162 1 : break;
163 0 : case GVOT_NORMAL:
164 0 : oOpts.outputMode = viewshed::OutputMode::Normal;
165 0 : break;
166 : }
167 :
168 1 : if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))
169 : {
170 0 : CPLError(CE_Failure, CPLE_AppDefined,
171 : "dfVisibleVal out of range. Must be [0, 255].");
172 0 : return nullptr;
173 : }
174 1 : if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))
175 : {
176 0 : CPLError(CE_Failure, CPLE_AppDefined,
177 : "dfInvisibleVal out of range. Must be [0, 255].");
178 0 : return nullptr;
179 : }
180 1 : if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))
181 : {
182 0 : CPLError(CE_Failure, CPLE_AppDefined,
183 : "dfOutOfRangeVal out of range. Must be [0, 255].");
184 0 : return nullptr;
185 : }
186 1 : oOpts.visibleVal = dfVisibleVal;
187 1 : oOpts.invisibleVal = dfInvisibleVal;
188 1 : oOpts.outOfRangeVal = dfOutOfRangeVal;
189 :
190 1 : gdal::viewshed::Viewshed v(oOpts);
191 :
192 1 : if (!pfnProgress)
193 1 : pfnProgress = GDALDummyProgress;
194 1 : v.run(hBand, pfnProgress, pProgressArg);
195 :
196 1 : return GDALDataset::FromHandle(v.output().release());
197 : }
198 :
199 : namespace gdal
200 : {
201 : namespace viewshed
202 : {
203 :
204 : namespace
205 : {
206 :
207 37 : bool getTransforms(GDALRasterBand &band, double *pFwdTransform,
208 : double *pRevTransform)
209 : {
210 37 : band.GetDataset()->GetGeoTransform(pFwdTransform);
211 37 : if (!GDALInvGeoTransform(pFwdTransform, pRevTransform))
212 : {
213 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
214 0 : return false;
215 : }
216 37 : return true;
217 : }
218 :
219 : } // unnamed namespace
220 :
221 37 : Viewshed::Viewshed(const Options &opts) : oOpts{opts}
222 : {
223 37 : }
224 :
225 : Viewshed::~Viewshed() = default;
226 :
227 : /// Calculate the extent of the output raster in terms of the input raster and
228 : /// save the input raster extent.
229 : ///
230 : /// @return false on error, true otherwise
231 37 : bool Viewshed::calcExtents(int nX, int nY,
232 : const std::array<double, 6> &adfInvTransform)
233 : {
234 : // We start with the assumption that the output size matches the input.
235 37 : oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
236 37 : oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
237 :
238 37 : if (!oOutExtent.contains(nX, nY))
239 8 : CPLError(CE_Warning, CPLE_AppDefined,
240 : "NOTE: The observer location falls outside of the DEM area");
241 :
242 37 : constexpr double EPSILON = 1e-8;
243 37 : if (oOpts.maxDistance > 0)
244 : {
245 : //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
246 : // Find the distance in the direction of the transformed unit vector in the X and Y
247 : // directions and use those factors to determine the limiting values in the raster space.
248 2 : int nXStart = static_cast<int>(
249 2 : std::floor(nX - adfInvTransform[1] * oOpts.maxDistance + EPSILON));
250 2 : int nXStop = static_cast<int>(
251 2 : std::ceil(nX + adfInvTransform[1] * oOpts.maxDistance - EPSILON) +
252 2 : 1);
253 : int nYStart =
254 2 : static_cast<int>(std::floor(
255 2 : nY - std::fabs(adfInvTransform[5]) * oOpts.maxDistance +
256 2 : EPSILON)) -
257 2 : (adfInvTransform[5] > 0 ? 1 : 0);
258 2 : int nYStop = static_cast<int>(
259 2 : std::ceil(nY + std::fabs(adfInvTransform[5]) * oOpts.maxDistance -
260 2 : EPSILON) +
261 2 : (adfInvTransform[5] < 0 ? 1 : 0));
262 :
263 : // If the limits are invalid, set the window size to zero to trigger the error below.
264 2 : if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
265 2 : nYStart >= oOutExtent.yStop || nYStop < 0)
266 : {
267 0 : oOutExtent = Window();
268 : }
269 : else
270 : {
271 2 : oOutExtent.xStart = std::max(nXStart, 0);
272 2 : oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
273 :
274 2 : oOutExtent.yStart = std::max(nYStart, 0);
275 2 : oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
276 : }
277 : }
278 :
279 37 : if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
280 : {
281 0 : CPLError(CE_Failure, CPLE_AppDefined,
282 : "Invalid target raster size due to transform "
283 : "and/or distance limitation.");
284 0 : return false;
285 : }
286 :
287 : // normalize horizontal index to [ 0, oOutExtent.xSize() )
288 37 : oCurExtent = oOutExtent;
289 37 : oCurExtent.shiftX(-oOutExtent.xStart);
290 :
291 37 : return true;
292 : }
293 :
294 : /// Compute the viewshed of a raster band.
295 : ///
296 : /// @param band Pointer to the raster band to be processed.
297 : /// @param pfnProgress Pointer to the progress function. Can be null.
298 : /// @param pProgressArg Argument passed to the progress function
299 : /// @return True on success, false otherwise.
300 37 : bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
301 : void *pProgressArg)
302 : {
303 37 : pSrcBand = static_cast<GDALRasterBand *>(band);
304 :
305 : std::array<double, 6> adfFwdTransform;
306 : std::array<double, 6> adfInvTransform;
307 37 : if (!getTransforms(*pSrcBand, adfFwdTransform.data(),
308 : adfInvTransform.data()))
309 0 : return false;
310 :
311 : // calculate observer position
312 : double dfX, dfY;
313 37 : GDALApplyGeoTransform(adfInvTransform.data(), oOpts.observer.x,
314 : oOpts.observer.y, &dfX, &dfY);
315 37 : if (!GDALIsValueInRange<int>(dfX))
316 : {
317 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
318 0 : return false;
319 : }
320 37 : if (!GDALIsValueInRange<int>(dfY))
321 : {
322 0 : CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
323 0 : return false;
324 : }
325 37 : int nX = static_cast<int>(dfX);
326 37 : int nY = static_cast<int>(dfY);
327 :
328 : // Must calculate extents in order to make the output dataset.
329 37 : if (!calcExtents(nX, nY, adfInvTransform))
330 0 : return false;
331 :
332 37 : poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
333 37 : if (!poDstDS)
334 2 : return false;
335 :
336 : // Create the progress reporter.
337 70 : Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());
338 :
339 : // Execute the viewshed algorithm.
340 35 : GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
341 35 : ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
342 70 : oCurExtent, oOpts, oProgress);
343 35 : executor.run();
344 35 : oProgress.emit(1);
345 35 : return static_cast<bool>(poDstDS);
346 : }
347 :
348 : } // namespace viewshed
349 : } // namespace gdal
|
