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("dst-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 98 : .SetMinValueIncluded(0)
141 98 : .SetMaxValueIncluded(255);
142 : AddArg("observer-spacing", 0, _("Cell Spacing between observers"),
143 196 : &m_opts.observerSpacing)
144 98 : .SetDefault(m_opts.observerSpacing)
145 98 : .SetMinValueIncluded(1);
146 :
147 98 : m_numThreadsStr = std::to_string(m_numThreads);
148 98 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
149 98 : }
150 :
151 : /************************************************************************/
152 : /* GDALRasterViewshedAlgorithm::RunStep() */
153 : /************************************************************************/
154 :
155 18 : bool GDALRasterViewshedAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
156 : {
157 18 : auto pfnProgress = ctxt.m_pfnProgress;
158 18 : auto pProgressData = ctxt.m_pProgressData;
159 18 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
160 18 : CPLAssert(poSrcDS);
161 18 : CPLAssert(!m_outputDataset.GetDatasetRef());
162 :
163 18 : GDALRasterBandH sdBand = nullptr;
164 18 : if (auto sdDataset = m_sdFilename.GetDatasetRef())
165 : {
166 2 : if (sdDataset->GetRasterCount() == 0)
167 : {
168 1 : ReportError(
169 : CE_Failure, CPLE_AppDefined,
170 : "The standard deviation dataset must have one raster band");
171 1 : return false;
172 : }
173 1 : sdBand = GDALRasterBand::FromHandle(sdDataset->GetRasterBand(1));
174 : }
175 :
176 17 : if (GetArg("height")->IsExplicitlySet())
177 : {
178 3 : if (m_observerPos.size() == 3)
179 : {
180 0 : ReportError(CE_Failure, CPLE_AppDefined,
181 : "Height can't be specified in both 'position' and "
182 : "'height' arguments");
183 0 : return false;
184 : }
185 : }
186 :
187 17 : if (m_observerPos.size())
188 : {
189 14 : m_opts.observer.x = m_observerPos[0];
190 14 : m_opts.observer.y = m_observerPos[1];
191 14 : if (m_observerPos.size() == 3)
192 14 : m_opts.observer.z = m_observerPos[2];
193 : else
194 0 : m_opts.observer.z = 2;
195 : }
196 :
197 17 : if (!GetArg("curvature-coefficient")->IsExplicitlySet())
198 : {
199 16 : m_opts.curveCoeff = gdal::viewshed::adjustCurveCoeff(
200 : m_opts.curveCoeff, GDALDataset::ToHandle(poSrcDS));
201 : }
202 :
203 17 : if (m_outputMode == "normal")
204 11 : m_opts.outputMode = gdal::viewshed::OutputMode::Normal;
205 6 : else if (m_outputMode == "DEM")
206 1 : m_opts.outputMode = gdal::viewshed::OutputMode::DEM;
207 5 : else if (m_outputMode == "ground")
208 2 : m_opts.outputMode = gdal::viewshed::OutputMode::Ground;
209 3 : else if (m_outputMode == "cumulative")
210 3 : m_opts.outputMode = gdal::viewshed::OutputMode::Cumulative;
211 :
212 17 : m_opts.numJobs = static_cast<uint8_t>(std::clamp(m_numThreads, 0, 255));
213 :
214 : m_opts.outputFilename =
215 17 : CPLGenerateTempFilenameSafe(
216 51 : CPLGetBasenameSafe(poSrcDS->GetDescription()).c_str()) +
217 17 : ".tif";
218 17 : m_opts.outputFormat = "GTiff";
219 :
220 17 : if (m_opts.outputMode == gdal::viewshed::OutputMode::Cumulative)
221 : {
222 : static const std::vector<std::string> badArgs{
223 : "visible-value", "invisible-value", "max-distance",
224 : "min-distance", "start-angle", "end-angle",
225 12 : "low-pitch", "high-pitch", "position"};
226 :
227 30 : for (const auto &arg : badArgs)
228 27 : if (GetArg(arg)->IsExplicitlySet())
229 : {
230 : std::string err =
231 0 : "Option '" + arg + "' can't be used in cumulative mode.";
232 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str());
233 0 : return false;
234 : }
235 :
236 3 : auto poSrcDriver = poSrcDS->GetDriver();
237 5 : if (EQUAL(poSrcDS->GetDescription(), "") || !poSrcDriver ||
238 2 : EQUAL(poSrcDriver->GetDescription(), "MEM"))
239 : {
240 1 : ReportError(
241 : CE_Failure, CPLE_AppDefined,
242 : "In cumulative mode, the input dataset must be opened by name");
243 1 : return false;
244 : }
245 4 : gdal::viewshed::Cumulative oViewshed(m_opts);
246 4 : const bool bSuccess = oViewshed.run(
247 2 : m_inputDataset[0].GetName().c_str(),
248 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
249 2 : if (bSuccess)
250 : {
251 2 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(
252 : GDALDataset::Open(m_opts.outputFilename.c_str(),
253 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
254 : nullptr, nullptr, nullptr)));
255 : }
256 : }
257 : else
258 : {
259 : static const std::vector<std::string> badArgs{
260 16 : "observer-spacing", GDAL_ARG_NAME_NUM_THREADS};
261 42 : for (const auto &arg : badArgs)
262 28 : if (GetArg(arg)->IsExplicitlySet())
263 : {
264 : std::string err =
265 0 : "Option '" + arg + "' can't be used in standard mode.";
266 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str());
267 0 : return false;
268 : }
269 15 : static const std::vector<std::string> goodArgs{"position"};
270 28 : for (const auto &arg : goodArgs)
271 14 : if (!GetArg(arg)->IsExplicitlySet())
272 : {
273 : std::string err =
274 0 : "Option '" + arg + "' must be specified in standard mode.";
275 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s", err.c_str());
276 0 : return false;
277 : }
278 :
279 28 : gdal::viewshed::Viewshed oViewshed(m_opts);
280 14 : const bool bSuccess = oViewshed.run(
281 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(m_band)), sdBand,
282 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
283 14 : if (bSuccess)
284 : {
285 14 : m_outputDataset.Set(oViewshed.output());
286 : }
287 : }
288 :
289 16 : auto poOutDS = m_outputDataset.GetDatasetRef();
290 16 : if (poOutDS && poOutDS->GetDescription()[0])
291 : {
292 : // In file systems that allow it (all but Windows...), we want to
293 : // delete the temporary file as soon as soon as possible after
294 : // having open it, so that if someone kills the process there are
295 : // no temp files left over. If that unlink() doesn't succeed
296 : // (on Windows), then the file will eventually be deleted when
297 : // poTmpDS is cleaned due to MarkSuppressOnClose().
298 16 : VSIUnlink(poOutDS->GetDescription());
299 16 : poOutDS->MarkSuppressOnClose();
300 : }
301 :
302 16 : return poOutDS != nullptr;
303 : }
304 :
305 : GDALRasterViewshedAlgorithmStandalone::
306 : ~GDALRasterViewshedAlgorithmStandalone() = default;
307 :
308 : //! @endcond
|
