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