/src/gdal/apps/gdalalg_raster_viewshed.cpp
Line | Count | Source |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: GDAL |
4 | | * Purpose: gdal "raster viewshed" subcommand |
5 | | * Author: Even Rouault <even dot rouault at spatialys.com> |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com> |
9 | | * |
10 | | * SPDX-License-Identifier: MIT |
11 | | ****************************************************************************/ |
12 | | |
13 | | #include "gdalalg_raster_viewshed.h" |
14 | | |
15 | | #include "cpl_conv.h" |
16 | | #include "cpl_vsi_virtual.h" |
17 | | |
18 | | #include "commonutils.h" |
19 | | #include "gdal_priv.h" |
20 | | |
21 | | #include "viewshed/cumulative.h" |
22 | | #include "viewshed/viewshed.h" |
23 | | |
24 | | #include <algorithm> |
25 | | |
26 | | //! @cond Doxygen_Suppress |
27 | | |
28 | | #ifndef _ |
29 | 0 | #define _(x) (x) |
30 | | #endif |
31 | | |
32 | | /************************************************************************/ |
33 | | /* GDALRasterViewshedAlgorithm::GDALRasterViewshedAlgorithm() */ |
34 | | /************************************************************************/ |
35 | | |
36 | | GDALRasterViewshedAlgorithm::GDALRasterViewshedAlgorithm(bool standaloneStep) |
37 | 0 | : GDALRasterPipelineNonNativelyStreamingAlgorithm(NAME, DESCRIPTION, |
38 | 0 | HELP_URL, standaloneStep) |
39 | 0 | { |
40 | 0 | AddArg("position", 'p', _("Observer position"), &m_observerPos) |
41 | 0 | .AddAlias("pos") |
42 | 0 | .SetMetaVar("<X,Y> or <X,Y,H>") |
43 | 0 | .SetMinCount(2) |
44 | 0 | .SetMaxCount(3) |
45 | 0 | .SetRepeatedArgAllowed(false); |
46 | 0 | AddArg("height", 'z', _("Observer height"), &m_opts.observer.z); |
47 | |
|
48 | 0 | auto &sdFilenameArg = |
49 | 0 | AddArg("sd-filename", 0, _("Filename of standard-deviation raster"), |
50 | 0 | &m_sdFilename, GDAL_OF_RASTER); |
51 | 0 | SetAutoCompleteFunctionForFilename(sdFilenameArg, GDAL_OF_RASTER); |
52 | |
|
53 | 0 | AddArg("target-height", 0, |
54 | 0 | _("Height of the target above the DEM surface in the height unit of " |
55 | 0 | "the DEM."), |
56 | 0 | &m_opts.targetHeight) |
57 | 0 | .SetDefault(m_opts.targetHeight); |
58 | 0 | AddArg("mode", 0, _("Sets what information the output contains."), |
59 | 0 | &m_outputMode) |
60 | 0 | .SetChoices("normal", "DEM", "ground", "cumulative") |
61 | 0 | .SetDefault(m_outputMode); |
62 | |
|
63 | 0 | AddArg("max-distance", 0, |
64 | 0 | _("Maximum distance from observer to compute visibility. It is also " |
65 | 0 | "used to clamp the extent of the output raster."), |
66 | 0 | &m_opts.maxDistance) |
67 | 0 | .SetMinValueIncluded(0); |
68 | 0 | AddArg("min-distance", 0, |
69 | 0 | _("Mask all cells less than this distance from the observer. Must " |
70 | 0 | "be less " |
71 | 0 | "than 'max-distance'."), |
72 | 0 | &m_opts.minDistance) |
73 | 0 | .SetMinValueIncluded(0); |
74 | |
|
75 | 0 | AddArg("start-angle", 0, |
76 | 0 | _("Mask all cells outside of the arc ('start-angle', 'end-angle'). " |
77 | 0 | "Clockwise degrees " |
78 | 0 | "from north. Also used to clamp the extent of the output raster."), |
79 | 0 | &m_opts.startAngle) |
80 | 0 | .SetMinValueIncluded(0) |
81 | 0 | .SetMaxValueExcluded(360); |
82 | 0 | AddArg("end-angle", 0, |
83 | 0 | _("Mask all cells outside of the arc ('start-angle', 'end-angle'). " |
84 | 0 | "Clockwise degrees " |
85 | 0 | "from north. Also used to clamp the extent of the output raster."), |
86 | 0 | &m_opts.endAngle) |
87 | 0 | .SetMinValueIncluded(0) |
88 | 0 | .SetMaxValueExcluded(360); |
89 | |
|
90 | 0 | AddArg("high-pitch", 0, |
91 | 0 | _("Mark all cells out-of-range where the observable height would be " |
92 | 0 | "higher than the " |
93 | 0 | "'high-pitch' angle from the observer. Degrees from horizontal - " |
94 | 0 | "positive is up. " |
95 | 0 | "Must be greater than 'low-pitch'."), |
96 | 0 | &m_opts.highPitch) |
97 | 0 | .SetMaxValueIncluded(90) |
98 | 0 | .SetMinValueExcluded(-90); |
99 | 0 | AddArg("low-pitch", 0, |
100 | 0 | _("Bound observable height to be no lower than the 'low-pitch' " |
101 | 0 | "angle from the observer. " |
102 | 0 | "Degrees from horizontal - positive is up. Must be less than " |
103 | 0 | "'high-pitch'."), |
104 | 0 | &m_opts.lowPitch) |
105 | 0 | .SetMaxValueExcluded(90) |
106 | 0 | .SetMinValueIncluded(-90); |
107 | |
|
108 | 0 | AddArg("curvature-coefficient", 0, |
109 | 0 | _("Coefficient to consider the effect of the curvature and " |
110 | 0 | "refraction."), |
111 | 0 | &m_opts.curveCoeff) |
112 | 0 | .SetMinValueIncluded(0); |
113 | |
|
114 | 0 | AddBandArg(&m_band).SetDefault(m_band); |
115 | 0 | AddArg("visible-value", 0, _("Pixel value to set for visible areas"), |
116 | 0 | &m_opts.visibleVal) |
117 | 0 | .SetDefault(m_opts.visibleVal) |
118 | 0 | .SetMinValueIncluded(0) |
119 | 0 | .SetMaxValueIncluded(255); |
120 | 0 | AddArg("invisible-value", 0, _("Pixel value to set for invisible areas"), |
121 | 0 | &m_opts.invisibleVal) |
122 | 0 | .SetDefault(m_opts.invisibleVal) |
123 | 0 | .SetMinValueIncluded(0) |
124 | 0 | .SetMaxValueIncluded(255); |
125 | 0 | AddArg("maybe-visible-value", 0, |
126 | 0 | _("Pixel value to set for potentially visible areas"), |
127 | 0 | &m_opts.maybeVisibleVal) |
128 | 0 | .SetDefault(m_opts.maybeVisibleVal) |
129 | 0 | .SetMinValueIncluded(0) |
130 | 0 | .SetMaxValueIncluded(255); |
131 | 0 | AddArg("out-of-range-value", 0, |
132 | 0 | _("Pixel value to set for the cells that fall outside of the range " |
133 | 0 | "specified by the observer location and the maximum distance"), |
134 | 0 | &m_opts.outOfRangeVal) |
135 | 0 | .SetDefault(m_opts.outOfRangeVal) |
136 | 0 | .SetMinValueIncluded(0) |
137 | 0 | .SetMaxValueIncluded(255); |
138 | 0 | AddArg("dst-nodata", 0, |
139 | 0 | _("The value to be set for the cells in the output raster that have " |
140 | 0 | "no data."), |
141 | 0 | &m_opts.nodataVal) |
142 | 0 | .SetMinValueIncluded(0) |
143 | 0 | .SetMaxValueIncluded(255); |
144 | 0 | AddArg("observer-spacing", 0, _("Cell Spacing between observers"), |
145 | 0 | &m_opts.observerSpacing) |
146 | 0 | .SetDefault(m_opts.observerSpacing) |
147 | 0 | .SetMinValueIncluded(1); |
148 | |
|
149 | 0 | m_numThreadsStr = std::to_string(m_numThreads); |
150 | 0 | AddNumThreadsArg(&m_numThreads, &m_numThreadsStr); |
151 | 0 | } |
152 | | |
153 | | /************************************************************************/ |
154 | | /* GDALRasterViewshedAlgorithm::RunStep() */ |
155 | | /************************************************************************/ |
156 | | |
157 | | bool GDALRasterViewshedAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt) |
158 | 0 | { |
159 | 0 | auto pfnProgress = ctxt.m_pfnProgress; |
160 | 0 | auto pProgressData = ctxt.m_pProgressData; |
161 | 0 | auto poSrcDS = m_inputDataset[0].GetDatasetRef(); |
162 | 0 | CPLAssert(poSrcDS); |
163 | 0 | CPLAssert(!m_outputDataset.GetDatasetRef()); |
164 | | |
165 | 0 | GDALRasterBandH sdBand = nullptr; |
166 | 0 | if (auto sdDataset = m_sdFilename.GetDatasetRef()) |
167 | 0 | { |
168 | 0 | if (sdDataset->GetRasterCount() == 0) |
169 | 0 | { |
170 | 0 | ReportError( |
171 | 0 | CE_Failure, CPLE_AppDefined, |
172 | 0 | "The standard deviation dataset must have one raster band"); |
173 | 0 | return false; |
174 | 0 | } |
175 | 0 | sdBand = GDALRasterBand::FromHandle(sdDataset->GetRasterBand(1)); |
176 | 0 | } |
177 | | |
178 | 0 | if (GetArg("height")->IsExplicitlySet()) |
179 | 0 | { |
180 | 0 | if (m_observerPos.size() == 3) |
181 | 0 | { |
182 | 0 | ReportError(CE_Failure, CPLE_AppDefined, |
183 | 0 | "Height can't be specified in both 'position' and " |
184 | 0 | "'height' arguments"); |
185 | 0 | return false; |
186 | 0 | } |
187 | 0 | } |
188 | | |
189 | 0 | if (m_observerPos.size()) |
190 | 0 | { |
191 | 0 | m_opts.observer.x = m_observerPos[0]; |
192 | 0 | m_opts.observer.y = m_observerPos[1]; |
193 | 0 | if (m_observerPos.size() == 3) |
194 | 0 | m_opts.observer.z = m_observerPos[2]; |
195 | 0 | else |
196 | 0 | m_opts.observer.z = 2; |
197 | 0 | } |
198 | |
|
199 | 0 | if (!GetArg("curvature-coefficient")->IsExplicitlySet()) |
200 | 0 | { |
201 | 0 | m_opts.curveCoeff = gdal::viewshed::adjustCurveCoeff( |
202 | 0 | m_opts.curveCoeff, GDALDataset::ToHandle(poSrcDS)); |
203 | 0 | } |
204 | |
|
205 | 0 | if (m_outputMode == "normal") |
206 | 0 | m_opts.outputMode = gdal::viewshed::OutputMode::Normal; |
207 | 0 | else if (m_outputMode == "DEM") |
208 | 0 | m_opts.outputMode = gdal::viewshed::OutputMode::DEM; |
209 | 0 | else if (m_outputMode == "ground") |
210 | 0 | m_opts.outputMode = gdal::viewshed::OutputMode::Ground; |
211 | 0 | else if (m_outputMode == "cumulative") |
212 | 0 | m_opts.outputMode = gdal::viewshed::OutputMode::Cumulative; |
213 | |
|
214 | 0 | m_opts.numJobs = static_cast<uint8_t>(std::clamp(m_numThreads, 0, 255)); |
215 | |
|
216 | 0 | m_opts.outputFilename = |
217 | 0 | CPLGenerateTempFilenameSafe( |
218 | 0 | CPLGetBasenameSafe(poSrcDS->GetDescription()).c_str()) + |
219 | 0 | ".tif"; |
220 | 0 | m_opts.outputFormat = "GTiff"; |
221 | |
|
222 | 0 | if (m_opts.outputMode == gdal::viewshed::OutputMode::Cumulative) |
223 | 0 | { |
224 | 0 | static const std::vector<std::string> badArgs{ |
225 | 0 | "visible-value", "invisible-value", "max-distance", |
226 | 0 | "min-distance", "start-angle", "end-angle", |
227 | 0 | "low-pitch", "high-pitch", "position"}; |
228 | |
|
229 | 0 | for (const auto &arg : badArgs) |
230 | 0 | if (GetArg(arg)->IsExplicitlySet()) |
231 | 0 | { |
232 | 0 | std::string err = |
233 | 0 | "Option '" + arg + "' can't be used in cumulative mode."; |
234 | 0 | ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str()); |
235 | 0 | return false; |
236 | 0 | } |
237 | | |
238 | 0 | auto poSrcDriver = poSrcDS->GetDriver(); |
239 | 0 | if (EQUAL(poSrcDS->GetDescription(), "") || !poSrcDriver || |
240 | 0 | EQUAL(poSrcDriver->GetDescription(), "MEM")) |
241 | 0 | { |
242 | 0 | ReportError( |
243 | 0 | CE_Failure, CPLE_AppDefined, |
244 | 0 | "In cumulative mode, the input dataset must be opened by name"); |
245 | 0 | return false; |
246 | 0 | } |
247 | 0 | gdal::viewshed::Cumulative oViewshed(m_opts); |
248 | 0 | const bool bSuccess = oViewshed.run( |
249 | 0 | m_inputDataset[0].GetName().c_str(), |
250 | 0 | pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData); |
251 | 0 | if (bSuccess) |
252 | 0 | { |
253 | 0 | m_outputDataset.Set(std::unique_ptr<GDALDataset>( |
254 | 0 | GDALDataset::Open(m_opts.outputFilename.c_str(), |
255 | 0 | GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, |
256 | 0 | nullptr, nullptr, nullptr))); |
257 | 0 | } |
258 | 0 | } |
259 | 0 | else |
260 | 0 | { |
261 | 0 | static const std::vector<std::string> badArgs{ |
262 | 0 | "observer-spacing", GDAL_ARG_NAME_NUM_THREADS}; |
263 | 0 | for (const auto &arg : badArgs) |
264 | 0 | if (GetArg(arg)->IsExplicitlySet()) |
265 | 0 | { |
266 | 0 | std::string err = |
267 | 0 | "Option '" + arg + "' can't be used in standard mode."; |
268 | 0 | ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str()); |
269 | 0 | return false; |
270 | 0 | } |
271 | 0 | static const std::vector<std::string> goodArgs{"position"}; |
272 | 0 | for (const auto &arg : goodArgs) |
273 | 0 | if (!GetArg(arg)->IsExplicitlySet()) |
274 | 0 | { |
275 | 0 | std::string err = |
276 | 0 | "Option '" + arg + "' must be specified in standard mode."; |
277 | 0 | ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str()); |
278 | 0 | return false; |
279 | 0 | } |
280 | | |
281 | 0 | gdal::viewshed::Viewshed oViewshed(m_opts); |
282 | 0 | const bool bSuccess = oViewshed.run( |
283 | 0 | GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(m_band)), sdBand, |
284 | 0 | pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData); |
285 | 0 | if (bSuccess) |
286 | 0 | { |
287 | 0 | m_outputDataset.Set(oViewshed.output()); |
288 | 0 | } |
289 | 0 | } |
290 | | |
291 | 0 | auto poOutDS = m_outputDataset.GetDatasetRef(); |
292 | 0 | if (poOutDS && poOutDS->GetDescription()[0]) |
293 | 0 | { |
294 | | // In file systems that allow it (all but Windows...), we want to |
295 | | // delete the temporary file as soon as soon as possible after |
296 | | // having open it, so that if someone kills the process there are |
297 | | // no temp files left over. If that unlink() doesn't succeed |
298 | | // (on Windows), then the file will eventually be deleted when |
299 | | // poTmpDS is cleaned due to MarkSuppressOnClose(). |
300 | 0 | VSIUnlink(poOutDS->GetDescription()); |
301 | 0 | poOutDS->MarkSuppressOnClose(); |
302 | 0 | } |
303 | |
|
304 | 0 | return poOutDS != nullptr; |
305 | 0 | } |
306 | | |
307 | | GDALRasterViewshedAlgorithmStandalone:: |
308 | 0 | ~GDALRasterViewshedAlgorithmStandalone() = default; |
309 | | |
310 | | //! @endcond |