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 25 : GDALRasterViewshedAlgorithm::GDALRasterViewshedAlgorithm(bool standaloneStep)
37 : : GDALRasterPipelineNonNativelyStreamingAlgorithm(NAME, DESCRIPTION,
38 25 : HELP_URL, standaloneStep)
39 : {
40 50 : AddArg("position", 'p', _("Observer position"), &m_observerPos)
41 50 : .AddAlias("pos")
42 50 : .SetMetaVar("<X,Y> or <X,Y,H>")
43 25 : .SetMinCount(2)
44 25 : .SetMaxCount(3)
45 25 : .SetRequired()
46 25 : .SetRepeatedArgAllowed(false);
47 : AddArg("target-height", 0,
48 : _("Height of the target above the DEM surface in the height unit of "
49 : "the DEM."),
50 50 : &m_targetHeight)
51 25 : .SetDefault(m_targetHeight);
52 : AddArg("mode", 0, _("Sets what information the output contains."),
53 50 : &m_outputMode)
54 25 : .SetChoices("normal", "DEM", "ground", "cumulative")
55 25 : .SetDefault(m_outputMode);
56 :
57 : AddArg("max-distance", 0,
58 : _("Maximum distance from observer to compute visibility. It is also "
59 : "used to clamp the extent of the output raster."),
60 50 : &m_maxDistance)
61 25 : .SetMinValueIncluded(0);
62 : AddArg("curvature-coefficient", 0,
63 : _("Coefficient to consider the effect of the curvature and "
64 : "refraction."),
65 50 : &m_curveCoefficient)
66 25 : .SetMinValueIncluded(0);
67 :
68 25 : AddBandArg(&m_band).SetDefault(m_band);
69 : AddArg("visible-value", 0, _("Pixel value to set for visible areas"),
70 50 : &m_visibleVal)
71 25 : .SetDefault(m_visibleVal)
72 25 : .SetMinValueIncluded(0)
73 25 : .SetMaxValueIncluded(255);
74 : AddArg("invisible-value", 0, _("Pixel value to set for invisible areas"),
75 50 : &m_invisibleVal)
76 25 : .SetDefault(m_invisibleVal)
77 25 : .SetMinValueIncluded(0)
78 25 : .SetMaxValueIncluded(255);
79 : AddArg("out-of-range-value", 0,
80 : _("Pixel value to set for the cells that fall outside of the range "
81 : "specified by the observer location and the maximum distance"),
82 50 : &m_outOfRangeVal)
83 25 : .SetDefault(m_outOfRangeVal)
84 25 : .SetMinValueIncluded(0)
85 25 : .SetMaxValueIncluded(255);
86 : AddArg("dst-nodata", 0,
87 : _("The value to be set for the cells in the output raster that have "
88 : "no data."),
89 50 : &m_dstNoData)
90 25 : .SetMinValueIncluded(0)
91 25 : .SetMaxValueIncluded(255);
92 : AddArg("observer-spacing", 0, _("Cell Spacing between observers"),
93 50 : &m_observerSpacing)
94 25 : .SetDefault(m_observerSpacing)
95 25 : .SetMinValueIncluded(1);
96 :
97 25 : m_numThreadsStr = std::to_string(m_numThreads);
98 25 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
99 25 : }
100 :
101 : /************************************************************************/
102 : /* GDALRasterViewshedAlgorithm::RunStep() */
103 : /************************************************************************/
104 :
105 15 : bool GDALRasterViewshedAlgorithm::RunStep(GDALProgressFunc pfnProgress,
106 : void *pProgressData)
107 : {
108 15 : auto poSrcDS = m_inputDataset.GetDatasetRef();
109 15 : CPLAssert(poSrcDS);
110 15 : CPLAssert(!m_outputDataset.GetDatasetRef());
111 :
112 30 : gdal::viewshed::Options opts;
113 :
114 15 : opts.observer.x = m_observerPos[0];
115 15 : opts.observer.y = m_observerPos[1];
116 15 : if (m_observerPos.size() == 3)
117 15 : opts.observer.z = m_observerPos[2];
118 : else
119 0 : opts.observer.z = 2;
120 :
121 15 : opts.targetHeight = m_targetHeight;
122 :
123 15 : opts.maxDistance = m_maxDistance;
124 :
125 15 : opts.curveCoeff = m_curveCoefficient;
126 15 : if (!GetArg("curvature-coefficient")->IsExplicitlySet())
127 : {
128 14 : opts.curveCoeff = gdal::viewshed::adjustCurveCoeff(
129 : opts.curveCoeff,
130 : GDALDataset::ToHandle(m_inputDataset.GetDatasetRef()));
131 : }
132 :
133 15 : opts.visibleVal = m_visibleVal;
134 15 : opts.invisibleVal = m_invisibleVal;
135 15 : opts.outOfRangeVal = m_outOfRangeVal;
136 15 : opts.nodataVal = m_dstNoData;
137 :
138 15 : if (m_outputMode == "normal")
139 10 : opts.outputMode = gdal::viewshed::OutputMode::Normal;
140 5 : else if (m_outputMode == "DEM")
141 1 : opts.outputMode = gdal::viewshed::OutputMode::DEM;
142 4 : else if (m_outputMode == "ground")
143 1 : opts.outputMode = gdal::viewshed::OutputMode::Ground;
144 3 : else if (m_outputMode == "cumulative")
145 3 : opts.outputMode = gdal::viewshed::OutputMode::Cumulative;
146 :
147 15 : opts.observerSpacing = m_observerSpacing;
148 15 : opts.numJobs = static_cast<uint8_t>(std::clamp(m_numThreads, 0, 255));
149 :
150 : opts.outputFilename =
151 15 : CPLGenerateTempFilenameSafe(
152 45 : CPLGetBasenameSafe(poSrcDS->GetDescription()).c_str()) +
153 15 : ".tif";
154 15 : opts.outputFormat = "GTiff";
155 :
156 15 : if (opts.outputMode == gdal::viewshed::OutputMode::Cumulative)
157 : {
158 3 : auto poSrcDriver = poSrcDS->GetDriver();
159 5 : if (EQUAL(poSrcDS->GetDescription(), "") || !poSrcDriver ||
160 2 : EQUAL(poSrcDriver->GetDescription(), "MEM"))
161 : {
162 1 : ReportError(
163 : CE_Failure, CPLE_AppDefined,
164 : "In cumulative mode, the input dataset must be opened by name");
165 1 : return false;
166 : }
167 4 : gdal::viewshed::Cumulative oViewshed(opts);
168 4 : const bool bSuccess = oViewshed.run(
169 2 : m_inputDataset.GetName().c_str(),
170 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
171 2 : if (bSuccess)
172 : {
173 2 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(
174 : GDALDataset::Open(opts.outputFilename.c_str(),
175 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
176 : nullptr, nullptr, nullptr)));
177 : }
178 : }
179 : else
180 : {
181 24 : gdal::viewshed::Viewshed oViewshed(opts);
182 12 : const bool bSuccess = oViewshed.run(
183 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(m_band)),
184 : pfnProgress ? pfnProgress : GDALDummyProgress, pProgressData);
185 12 : if (bSuccess)
186 : {
187 12 : m_outputDataset.Set(oViewshed.output());
188 : }
189 : }
190 :
191 14 : auto poOutDS = m_outputDataset.GetDatasetRef();
192 14 : if (poOutDS && poOutDS->GetDescription()[0])
193 : {
194 : // In file systems that allow it (all but Windows...), we want to
195 : // delete the temporary file as soon as soon as possible after
196 : // having open it, so that if someone kills the process there are
197 : // no temp files left over. If that unlink() doesn't succeed
198 : // (on Windows), then the file will eventually be deleted when
199 : // poTmpDS is cleaned due to MarkSuppressOnClose().
200 14 : VSIUnlink(poOutDS->GetDescription());
201 14 : poOutDS->MarkSuppressOnClose();
202 : }
203 :
204 14 : return poOutDS != nullptr;
205 : }
206 :
207 : //! @endcond
|
