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