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 62 : GDALRasterViewshedAlgorithm::GDALRasterViewshedAlgorithm(bool standaloneStep)
37 : : GDALRasterPipelineNonNativelyStreamingAlgorithm(NAME, DESCRIPTION,
38 62 : HELP_URL, standaloneStep)
39 : {
40 124 : AddArg("position", 'p', _("Observer position"), &m_observerPos)
41 124 : .AddAlias("pos")
42 124 : .SetMetaVar("<X,Y> or <X,Y,H>")
43 62 : .SetMinCount(2)
44 62 : .SetMaxCount(3)
45 62 : .SetRepeatedArgAllowed(false);
46 62 : AddArg("height", 'z', _("Observer height"), &m_opts.observer.z);
47 :
48 : auto &sdFilenameArg =
49 : AddArg("sd-filename", 0, _("Filename of standard-deviation raster"),
50 62 : &m_sdFilename, GDAL_OF_RASTER);
51 62 : 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 124 : &m_opts.targetHeight)
57 62 : .SetDefault(m_opts.targetHeight);
58 : AddArg("mode", 0, _("Sets what information the output contains."),
59 124 : &m_outputMode)
60 62 : .SetChoices("normal", "DEM", "ground", "cumulative")
61 62 : .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 124 : &m_opts.maxDistance)
67 62 : .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 124 : &m_opts.minDistance)
73 62 : .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 124 : &m_opts.startAngle)
80 62 : .SetMinValueIncluded(0)
81 62 : .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 124 : &m_opts.endAngle)
87 62 : .SetMinValueIncluded(0)
88 62 : .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 124 : &m_opts.highPitch)
97 62 : .SetMaxValueIncluded(90)
98 62 : .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 124 : &m_opts.lowPitch)
105 62 : .SetMaxValueExcluded(90)
106 62 : .SetMinValueIncluded(-90);
107 :
108 : AddArg("curvature-coefficient", 0,
109 : _("Coefficient to consider the effect of the curvature and "
110 : "refraction."),
111 124 : &m_opts.curveCoeff)
112 62 : .SetMinValueIncluded(0);
113 :
114 62 : AddBandArg(&m_band).SetDefault(m_band);
115 : AddArg("visible-value", 0, _("Pixel value to set for visible areas"),
116 124 : &m_opts.visibleVal)
117 62 : .SetDefault(m_opts.visibleVal)
118 62 : .SetMinValueIncluded(0)
119 62 : .SetMaxValueIncluded(255);
120 : AddArg("invisible-value", 0, _("Pixel value to set for invisible areas"),
121 124 : &m_opts.invisibleVal)
122 62 : .SetDefault(m_opts.invisibleVal)
123 62 : .SetMinValueIncluded(0)
124 62 : .SetMaxValueIncluded(255);
125 : AddArg("maybe-visible-value", 0,
126 : _("Pixel value to set for potentially visible areas"),
127 124 : &m_opts.maybeVisibleVal)
128 62 : .SetDefault(m_opts.maybeVisibleVal)
129 62 : .SetMinValueIncluded(0)
130 62 : .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 124 : &m_opts.outOfRangeVal)
135 62 : .SetDefault(m_opts.outOfRangeVal)
136 62 : .SetMinValueIncluded(0)
137 62 : .SetMaxValueIncluded(255);
138 : AddArg("dst-nodata", 0,
139 : _("The value to be set for the cells in the output raster that have "
140 : "no data."),
141 124 : &m_opts.nodataVal)
142 62 : .SetMinValueIncluded(0)
143 62 : .SetMaxValueIncluded(255);
144 : AddArg("observer-spacing", 0, _("Cell Spacing between observers"),
145 124 : &m_opts.observerSpacing)
146 62 : .SetDefault(m_opts.observerSpacing)
147 62 : .SetMinValueIncluded(1);
148 :
149 62 : m_numThreadsStr = std::to_string(m_numThreads);
150 62 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
151 62 : }
152 :
153 : /************************************************************************/
154 : /* GDALRasterViewshedAlgorithm::RunStep() */
155 : /************************************************************************/
156 :
157 17 : bool GDALRasterViewshedAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
158 : {
159 17 : auto pfnProgress = ctxt.m_pfnProgress;
160 17 : auto pProgressData = ctxt.m_pProgressData;
161 17 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
162 17 : CPLAssert(poSrcDS);
163 17 : CPLAssert(!m_outputDataset.GetDatasetRef());
164 :
165 17 : GDALRasterBandH sdBand = nullptr;
166 17 : if (auto sdDataset = m_sdFilename.GetDatasetRef())
167 : {
168 2 : if (sdDataset->GetRasterCount() == 0)
169 : {
170 1 : ReportError(
171 : CE_Failure, CPLE_AppDefined,
172 : "The standard deviation dataset must have one raster band");
173 1 : return false;
174 : }
175 1 : sdBand = GDALRasterBand::FromHandle(sdDataset->GetRasterBand(1));
176 : }
177 :
178 16 : if (GetArg("height")->IsExplicitlySet())
179 : {
180 3 : if (m_observerPos.size() == 3)
181 : {
182 0 : ReportError(CE_Failure, CPLE_AppDefined,
183 : "Height can't be specified in both 'position' and "
184 : "'height' arguments");
185 0 : return false;
186 : }
187 : }
188 :
189 16 : if (m_observerPos.size())
190 : {
191 13 : m_opts.observer.x = m_observerPos[0];
192 13 : m_opts.observer.y = m_observerPos[1];
193 13 : if (m_observerPos.size() == 3)
194 13 : m_opts.observer.z = m_observerPos[2];
195 : else
196 0 : m_opts.observer.z = 2;
197 : }
198 :
199 16 : if (!GetArg("curvature-coefficient")->IsExplicitlySet())
200 : {
201 15 : m_opts.curveCoeff = gdal::viewshed::adjustCurveCoeff(
202 : m_opts.curveCoeff, GDALDataset::ToHandle(poSrcDS));
203 : }
204 :
205 16 : if (m_outputMode == "normal")
206 11 : m_opts.outputMode = gdal::viewshed::OutputMode::Normal;
207 5 : else if (m_outputMode == "DEM")
208 1 : m_opts.outputMode = gdal::viewshed::OutputMode::DEM;
209 4 : else if (m_outputMode == "ground")
210 1 : m_opts.outputMode = gdal::viewshed::OutputMode::Ground;
211 3 : else if (m_outputMode == "cumulative")
212 3 : m_opts.outputMode = gdal::viewshed::OutputMode::Cumulative;
213 :
214 16 : m_opts.numJobs = static_cast<uint8_t>(std::clamp(m_numThreads, 0, 255));
215 :
216 : m_opts.outputFilename =
217 16 : CPLGenerateTempFilenameSafe(
218 48 : CPLGetBasenameSafe(poSrcDS->GetDescription()).c_str()) +
219 16 : ".tif";
220 16 : m_opts.outputFormat = "GTiff";
221 :
222 16 : if (m_opts.outputMode == gdal::viewshed::OutputMode::Cumulative)
223 : {
224 : static const std::vector<std::string> badArgs{
225 : "visible-value", "invisible-value", "max-distance",
226 : "min-distance", "start-angle", "end-angle",
227 12 : "low-pitch", "high-pitch", "position"};
228 :
229 30 : for (const auto &arg : badArgs)
230 27 : if (GetArg(arg)->IsExplicitlySet())
231 : {
232 : 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 : }
237 :
238 3 : auto poSrcDriver = poSrcDS->GetDriver();
239 5 : if (EQUAL(poSrcDS->GetDescription(), "") || !poSrcDriver ||
240 2 : EQUAL(poSrcDriver->GetDescription(), "MEM"))
241 : {
242 1 : ReportError(
243 : CE_Failure, CPLE_AppDefined,
244 : "In cumulative mode, the input dataset must be opened by name");
245 1 : return false;
246 : }
247 4 : gdal::viewshed::Cumulative oViewshed(m_opts);
248 4 : const bool bSuccess = oViewshed.run(
249 2 : m_inputDataset[0].GetName().c_str(),
250 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
251 2 : if (bSuccess)
252 : {
253 2 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(
254 : GDALDataset::Open(m_opts.outputFilename.c_str(),
255 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
256 : nullptr, nullptr, nullptr)));
257 : }
258 : }
259 : else
260 : {
261 : static const std::vector<std::string> badArgs{
262 15 : "observer-spacing", GDAL_ARG_NAME_NUM_THREADS};
263 39 : for (const auto &arg : badArgs)
264 26 : if (GetArg(arg)->IsExplicitlySet())
265 : {
266 : 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 : }
271 14 : static const std::vector<std::string> goodArgs{"position"};
272 26 : for (const auto &arg : goodArgs)
273 13 : if (!GetArg(arg)->IsExplicitlySet())
274 : {
275 : 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 : }
280 :
281 26 : gdal::viewshed::Viewshed oViewshed(m_opts);
282 13 : const bool bSuccess = oViewshed.run(
283 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(m_band)), sdBand,
284 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
285 13 : if (bSuccess)
286 : {
287 13 : m_outputDataset.Set(oViewshed.output());
288 : }
289 : }
290 :
291 15 : auto poOutDS = m_outputDataset.GetDatasetRef();
292 15 : if (poOutDS && poOutDS->GetDescription()[0])
293 : {
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 15 : VSIUnlink(poOutDS->GetDescription());
301 15 : poOutDS->MarkSuppressOnClose();
302 : }
303 :
304 15 : return poOutDS != nullptr;
305 : }
306 :
307 : GDALRasterViewshedAlgorithmStandalone::
308 : ~GDALRasterViewshedAlgorithmStandalone() = default;
309 :
310 : //! @endcond
|
